• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Nonlinear Responses of Oceanic Temperature to Wind Stress Anomalies in Tropical Pacific and Indian Oceans:A Study Based on Numerical Experiments with an OGCM

    2015-01-05 02:01:57HUALijuan華莉娟andYUYongqiang俞永強
    Journal of Meteorological Research 2015年4期

    HUA Lijuan(華莉娟)and YU Yongqiang(俞永強)

    1 State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics,Institute of Atmospheric Physics,Chinese Academy of Sciences,Beijing 100029

    2 College of Earth Sciences,University of Chinese Academy of Sciences,Beijing 100049

    Nonlinear Responses of Oceanic Temperature to Wind Stress Anomalies in Tropical Pacific and Indian Oceans:A Study Based on Numerical Experiments with an OGCM

    HUA Lijuan1,2(華莉娟)and YU Yongqiang2?(俞永強)

    1 State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics,Institute of Atmospheric Physics,Chinese Academy of Sciences,Beijing 100029

    2 College of Earth Sciences,University of Chinese Academy of Sciences,Beijing 100049

    As a highly nonlinear dynamic system,oceanic general circulation models(OGCMs)usually exhibit nonlinear responses to prescribed wind stress forcing.To explore mechanisms for these nonlinear responses, we designed and conducted three idealized numerical experiments with an OGCM with modified wind stress forcing.In the experiments,the climatological mean wind stress was identical,and the only differences in external forcing were wind stress anomalies.The wind anomalies were set to zero in a control run,and the observed wind stress anomalies with and without reversed signs were superimposed on the mean climatology in two sensitivity experiments.Forced by the prescribed wind stress anomalies in sensitivity runs,the OGCM well reproduced the El Ni?no-Southern Oscillation(ENSO)and the Pacific and Indian Ocean Dipole(IOD) in the Indian Ocean,as wellas the asymmetry between positive and negative phases of these modes.Relative to the control run,the two sensitivity runs exhibited almost identical changes in the mean climate state, although the wind stress anomalies were reversed in these two experiments.Thus,it was concluded that the asymmetry of wind stress anomalies contributes only slightly to the mean state changes and ocean internal dynamics was the main contributor.Further heat budget analysis suggested that nonlinear temperature advection terms,including both mean advection and perturbed advection,favor the ENSO/IOD rectified effect on the mean state.

    ENSO,IOD,numerical simulation,asymmetry,rectification,nonlinearity

    1.Introduction

    It is well known that the probability density function of the observed tropical Pacific sea temperature does not always follow a standard normal distribution(Burgers and Stephenson,1999;An and Jin, 2004;Zhang et al.,2009).Instead,it shows a significantly skewed distribution in most regions.For example,in the equatorial central and eastern Pacific,both the strength and duration of warm temperature anomalies in El Ni?no events are markedly different from those in their counterparts;this difference is commonly referred to as El Ni?no-Southern Oscillation(ENSO)asymmetry(Zebiak and Cane,1987; Okumura and Deser,2010).Further studies indicated that the asymmetry of the temperature anomalies is closely related to nonlinear dynamical heating(NDH) (An and Jin,2004;An et al.,2005;Su et al.,2010). Furthermore,ENSO asymmetry can be caused by high-frequency winds through a nonlinear rectification process,which increases the skewness of sea surface temperature(SST)anomalies in the eastern equato-rial Pacific(Rong et al.,2011).Other researchers have reported that the atmospheric responses to SST anomalies of the same magnitudes but opposite signs are asymmetric,which might,in turn,reinforce the original SST asymmetry(Ohba and Ueda,2009). The interaction between ENSO and the climatological mean state is another important issue,but the specifics of this interaction remain unknown.Through numerical experiments with simplified or fully coupled models,many studies have emphasized the impact of the mean state on ENSO events,especially in modulations of ENSO period and amplitude(An and Jin,2001;Fedorov and Philander,2001).The climatological basic state of the tropical atmosphere can influence ENSO events by modulating westerly winds over the tropical Pacific(Yan and Zhang,2002).However,recent studies have begun to consider the time-mean effect of ENSO events(Sun et al.,2012;Sun et al.,2014;Hua et al.,2015).Sun(2003)found that stronger ENSO variations attempt to cool the western Pacific and warm the eastern Pacific.Sun and Zhang(2006)further noted that the time-mean effect of ENSO would cool the warm pool while warming the thermocline water and broad region in the surface eastern Pacific.Some of these studies have been empirically based(Rodger et al.,2004;Choi et al.,2012).These studies have suggested that decadal variability in the tropical Pacific SST probably results from the residual effect of ENSO.

    The asymmetry of temperature anomalies is not only observed in the tropical Pacific,but also in the Indian Ocean.The Indian Ocean Dipole(IOD)is a zonal mode of SST interannual variation(Saji et al.,1999).Observational(Saji and Yamagata,2003), modeled(Cai et al.,2005),and theoretical(Zhong et al.,2005)studies have indicated that the IOD is a seasonally dependent mode that develops rapidly in summer and reaches a mature phase in autumn. IOD events are significantly asymmetric,which is similar to ENSO events,and the strength of a positive IOD is greater than that of a negative IOD. Meanwhile,strong negative IOD temperature skewness anomalies exist in the eastern IOD in autumn (Hong et al.,2008).Hong et al.(2008)considered that the asymmetry in wind-advection-SST feedback, SST-cloud-radiation feedback,and wind-evaporation-SST feedback all contribute to this negative skewness. In particular,wind-advection-SST feedback is related to NDH.It shows that the NDH favors cold anomalies, rather than restraining warm anomalies,in the eastern IOD heat budget(Hong et al.,2008;Cai and Qiu, 2013).In addition,Zheng et al.(2010)noticed,using coupled model simulations,that SST-thermocline depth feedback is asymmetric.Furthermore,this asymmetry might induce the negative skewness of temperature anomalies in the equatorial southeastern Indian Ocean.Ogata et al.(2013)confirmed the above assumption through observations and an ocean general circulation model(OGCM)study.However, the rectified effect of IOD asymmetry on the climatological mean state has not been considered until now.

    Although many climate models,including standalone OGCMs and coupled general circulation models,have succeeded in reproducing the skewness of temperature anomalies,the mechanism of temperature skewness remains unclear.In this study,we conduct numerical experiments with an OGCM and then evaluate how well the model reproduces the observed asymmetry of temperature anomalies in two ocean basins.Many studies have shown a close relationship between temperature asymmetry and climatological mean state;however,it is difficult to identify a physical mechanism and causality.In our study,the OGCM willbe forced by the climatological mean wind and anomalous wind on an interannualtimescale.The simulated mean state is clearly different between the two experiments.Therefore,we can compare the differences between the two experiments and further explore how the physical mechanism influences the temperature asymmetry and mean state.Subsequently, the numerical experiments would help us comprehend the asymmetry and rectification resulting from the external forcing or ocean dynamics.

    2.Methodology

    2.1 Model

    The OGCM used in this study is a climate sys-tem ocean model,LICOM2.0(Liu et al.,2004a;Liu et al.,2012),which was developed at the State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics,Institute of Atmospheric Physics.This model has been widely applied in many studies(Liu et al.,2004b;Yu et al., 2011;Li et al.,2012).It is a global ocean model that covers the region from 75°S to 88°N,and the North Pole is modeled as an isolated island.In particular,the model’s dynamical framework is based on a latitudelongitude grid system with 1°×1°horizontal resolution;however,the meridional resolution increases to 0.5°between 10°S and 10°N,and there are 30 layers with 15 equal-depth levels in the upper 150 m. The model adopts some new physical parameterization schemes,such as a new turbulent mixing scheme (Canuto and Dubovikov,2005),solar radiation penetration,and an improved isopycnal-mixing scheme (Gent and Williams,1990;Large et al.,1997).It can well reproduce the basic temporal and spatial structures of the observed ENSO and IOD,especially the significant asymmetry between El Ni?no and La Ni?na events,as well as positive and negative IOD events (Hua et al.,2010).By using climatological monthlymean wind stress,LICOM2.0 was first integrated from a motionless state to an equilibrium state through a 500-yr spin up.The end state of year 500 was then used as the initial condition for the numerical experiments that will be described below.

    2.2 Experiment design

    In the two sensitivity runs,the surface wind stresses were climatological winds plus the interannual monthly anomalies,and the surface boundary conditions for temperature and salinity were the same as those for the control run.In the first sensitivity run(hereafter referred to as WSA),which was the same as the sensitivity run described by Sun et al.(2014),the interannual monthly anomalies during 1958-2001(obtained from the European Centre for Medium-Range Weather Forecasts 40-yr reanalysis(ERA40;Uppala et al.,2005))over the tropical Pacific(30°N-30°S,120°E-80°W)and Indian Ocean (30°N-30°S,30°-120°E)were superimposed on the climatological surface wind stress.The second sensitivity run(hereafter referred to as WSA?R)was the same as WSA except that the signs of the interannual monthly anomalies were reversed.Notably,because of the restoring boundary conditions,the simulated surface temperature and salinity were restored to the climatologicalmean value with a seasonalcycle.Then, as a negative feedback,the restoring boundary conditions would have a damping effect on any perturba-tions in SST and salinity that deviated from the climatological mean state.Consequently,the anomalies in temperature and other variables and the skewness of the anomalies could only result from the wind stress forcing or the ocean dynamical process in the numerical experiments.

    According to the experimental design,the simulated temperature asymmetry,as well as the ENSO/ IOD rectification,should originate only from the external forcing and internal nonlinear dynamics in the ocean.Comparison between the two sensitivity runs may help to quantify the relative contributions of the above two factors.The ocean reanalysis data from simple ocean data assimilation(SODA;Carton and Giese,2008)were compared with model simulations. The ocean model of SODA has a horizontal resolution of 0.5°by 0.5°and 40 vertical levels.The monthly dataset used in this study was the temperature data from SODA2.0.2.Climatology fields were defined as the averages for the period 1958-2001.

    2.3 Diagnosis method

    The diagnosis methods used in the paper follow:

    2)The time mean temperature equation can be expressed as

    3.Model verification

    Forced with the observed wind stress and heat flux,the ocean model well captured the basic characteristics of the observed mean state.Figure 1 presents the climatological mean temperature and its annual standard deviation in the equatorial upper ocean(2°N-2°S)as a function of longitude and depth in the two basins,in which both the simulation and reanalyses were used from 1958 to 2001.These plots demonstrate that the model could capture the basic structure of the observed temperature,including the depth and tilt of the equatorial thermocline.In the Pacific,the model succeeds in describing the spatial pattern of the thermocline,which is deeper in the west and shallower in the east,and the large-scale pattern of the simulated temperature coincides well with the reanalysis;however,the simulated temperature near the western boundary is approximately 1℃higher than that obtained from the reanalysis.Meanwhile,an excessive cold tongue and cold bias are shown in the equatorial eastern Pacific in WSA.In the Indian Ocean,the model also well reproduces the spatial structure of the observed thermocline,which is shallower in the western basin and deeper in the easternbasin.Both the observed and simulated tilts of the thermocline are opposite to those in the Pacific because the annual averaged wind stress prevailing in the equatorial Indian Ocean is westerly,but that in the Pacific is easterly.The simulated depth of the mean thermocline is almost 120 m,which agrees well with the reanalysis.However,there still exists some model bias;for instance,the simulated vertical temperature gradient is smaller than that obtained from the reanalysis.In general,the overall pattern sim-ulated by LICOM2.0 resembles the reanalysis in the two basins(implying that the model has enough simulation ability),and a few discrepancies,such as cold bias in the equatorial eastern Pacific,are attributed to the coarse resolution of the model.The simulated standard deviation in the tropical Pacific agrees well with the reanalysis except for the smaller magnitude along the thermocline.In the Indian Ocean,the spatial pattern of the standard deviation is basically similar to the reanalysis,although there are differences in magnitude.

    Fig.1.Temperature(℃;contour)and annual standard deviation(>0.3 with color)in the equatorial upper Pacific and Indian Oceans(2°N-2°S)in(a,b)SODA,(c,d)WSA,and(e,f)the corresponding temperature(contour)and standard deviation(color)differences between WSA and SODA.

    Fig.2.(a,b)Nino3 and IOD index,(c,d)temperature anomalies in subsurface of eastern Pacific(2°N-2°S,100°W; 100 m)and eastern Indian Ocean(2°N-2°S,95°E;60 m),and(e,f)temperature anomalies in subsurface of western Pacific(2°N-2°S,160°W;100 m)and western Indian Ocean(2°N-2°S,45°E;120 m)from SODA data(black line),WSA (red line),and WSA?R(green line).

    One of the principal causes of the simulated sea surface temperature anomalies as a response to the observed wind stress is ENSO in the tropical Pacific and IOD in the Indian Ocean.The observed ENSO/IOD variability was well reproduced by the OGCM,as shown in Fig.2,which compares the time series of the Nino3/IOD index from the reanalysis with thatfrom the WSA/WSA?R experiments.The simulated time series of the Nino3 from the WSA/WSA?R indices have correlations of 0.87/-0.79 with those from the reanalysis and 0.66/-0.50 for IOD index.The time series of the Nino3 from WSA has correlations of -0.92 and-0.72 with those from WSA?R and from the IOD index,respectively.The correlation between the time series of Nino3 and the IOD index is 0.36, 0.29,and 0.27 for the reanalysis,WSA,and WSA?R, respectively.The simulated time series from WSA has a similar temporal evolution to that obtained from the reanalysis,but the simulated amplitude is weaker,especially for the IOD index.This discrepancy is expected,given that the thermal boundary conditions are simplified and that the SST is restored to the observed climatological datasets(Levitus and Boyer,1994;Levitus et al.,1994).Moreover,because of the stronger magnitude associated with positive ENSO/IOD phases than with negative phases, the positive asymmetry can be found in both Nino3 and IOD indices(WSA experiment).More details can be seen in Table 1.Both the Nino3 and IOD indices show the temporal evolution at the surface,but the changes in the temperature anomalies in the subsurface also play an important role in ENSO/IOD dynamics.Therefore,the time series of the subsurface temperature anomalies in the eastern/western Pacific and Indian Ocean are presented for SODA,WSA,and WSA?R in Fig.2.In addition,the skewness and standard deviation of subsurface temperature anomalies can be seen in Tables 2 and 3.

    4.Asymmetry

    4.1 Subsurface

    4.1.1 ENSO/IOD asymmetry

    It is wellknown that ENSO is the most important coupled ocean-atmosphere phenomenon in the tropical Pacific(Zebiak and Cane,1987;Battisti,1988)and the first EOF mode of the equatorial subsurface temperature anomalies is ENSO(Cai et al.,2003;Rodger et al.,2004;Choi et al.,2013).Therefore,temperature asymmetry in the equatorial Pacific is mostly induced by ENSO asymmetry.Similarly,Cai et al.(2008) found,using temperature data from Scripps Institute of Oceanography,that the first EOF mode of the subsurface temperature anomalies is described by an eastwest dipole in the equatorial Indian Ocean.With SODA reanalysis data,we reached the same conclusion (figure omitted).In other words,the IOD mode is the first EOF mode for subsurface temperature anomalies,and the subsurface temperature asymmetry in the equatorial Indian Ocean also results from the IOD asymmetry.As suggested by An et al.(2005),skewness is an effective metric to quantify the asymmetryoftemperature anomalies.Figure 3 illustrates temperature skewness as a function of longitude and depth in the upper equatorial ocean(2°N-2°S).The positive (negative)skewness implies that the strength of warm (cold)anomalies is larger than that of cold(warm) anomalies.The skewness from WSA coincides well with that of SODA in the upper ocean.The maximum positive skewness is at 100 m at around 100°W,and the largest negative skewness is at the same depth at around 160°E.However,the positive skewness occurs in the region(170°E-160°W;0-80 m)that is different from SODA;and a small amount of negative skewness around 140°Win the upper 40 m in WSA is not found in SODA.Although there exist some discrepancies,the model could still represent the most significant skewness in the subsurface.In WSA?R,the location ofthe extreme center is not exactly the same as that in WSA,but there is a common phenomenon at a depth of 40-160 m.Meanwhile,the quantitative difference between the skewness from WSA and WSA?R is not negligible.From WSA minus WSA?R(figure omitted),we found that the largest positive value was 2.5 and was located at 100 m at around 100°W,while the largest negative value was-2 at 100 m at around 160°E.This result means that both the positive skewness in the equatorial eastern Pacific and the negative skewness in the western Pacific are stronger in WSA than in WSA?R.

    Table 1.Skewness and standard deviation of time series of Nino3 and IOD index from SODAdata,WSA, and WSA?R experiments

    Table 2.As in Table 1,but for time series of temperature anomalies in subsurface of eastern Pacific(2°N-2°S,100°W;100 m)and eastern Indian Ocean(2°N-2°S,95°E;60 m)

    Table 3.As in Table 1,but for time series of temperature anomalies in subsurface of western Pacific(2°N-2°S,160°W;100 m)and western Indian Ocean(2°N-2°S,45°E;120 m)

    Fig.3.Skewness of temperature in equatorialupper Pacific and Indian Oceans(2°N-2°S)in(a,b)SODA,(c,d)WSA, and(e,f)WSA?R.

    In the Indian Ocean,the skewness results from WSA and WSA?R and those from SODA are similar in the upper 140 m.There exists negative skewness in the upper 100 m of the equatorial eastern basin and positive skewness at 100-140 m.The largest negative skewness is between 90°-100°E at 20-40 m for both simulations and SODA,but the simulated values are larger.The positive skewness for WSA and WSA?R is also greater than that for SODA.The quantitative difference between the skewness results of WSA and WSA?R in the Indian Ocean(figure omitted)revealed that the largest negative skewness was nearly-7 and was located at 140 m between 80°-90°E,with a negative minimum of around-1 at 20-40 m between 90°and 100°E.Thus,considering the skewness difference between WSA and WSA?R,the negative skewness was slightly larger in the upper eastern Indian Ocean in WSA;more significant positive skewness was exhibited below 100 m in the eastern basin in WSA?R.The quantitative difference between WSA and WSA?R cannot be ignored in the two basins,which means that the asymmetry of the wind stress anomalies would cause changes in the magnitude of the temperature skewness.Specifically,the thermocline responses to the reversed wind stress anomalies were asymmetric. The thermocline became shallower when easterly wind stress anomalies occurred and deeper when westerly anomalies ofthe same magnitudes occurred.The more significant positive skewness in WSA?R is attributed to the stronger oceanic responses to the westerly wind stress anomalies.Furthermore,even if the climatological mean wind stress was applied in the control run,the simulated temperature still exhibited some skewness that was entirely different from that obtained from the reanalysis,WSA,and WSA?R(figure omitted).We found that the large-scale spatial pattern in the control run was obviously different from that in other plots,for example,differences in the location of the significant skewness and in the magnitude of the skewness.Considering the above findings,we conclude that the temperature skewness in the upper 40-160 m in the Pacific,and in the upper 140 m in the Indian Ocean,is not associated with the skewness of the wind stress.Our results resemble Figs.6 and 7c in Ogata et al.(2013);the skewness was zero in numerical experiments with wind stress anomalies as the external forcing(Ogata et al.,2013),but the SST response was strongly negatively skewed.The numerical experimental design in our study was not identical to that of Ogata et al.,but the conclusions are similar. Therefore,the temperature asymmetry in the subsurface to the first order may be attributed to the internal dynamical response to the external forcing.

    4.1.2 Role of NDH in ENSO/IOD asymmetry

    Besides the skewness of the temperature anomalies,another metric to quantify asymmetry is the residual of temperature anomalies associated with positive and negative phases of ENSO/IOD events(Rodgers et al.,2004;Su et al.,2010;Ogata et al.,2013;Sun et al.,2014).We performed heat budget analysis in the positive and negative phases of ENSO/IOD events (figure omitted).The results showed that the zonal and meridional terms of NDH enhanced the positive temperature anomalies in the warm events in the upper equatorial eastern Pacific but weakened the negative anomalies in the cold events.The vertical term of NDH contributed greatly to temperature decrease in the warm events in the equatorial eastern Indian Ocean,but mitigated temperature increase in the cold events.Our findings based on the numerical experiments are consistent with those from reanalysis data by Su et al.(2010).

    4.2 Heat budget in the mixed layer

    Figure 4 illustrates the temperature asymmetry in the mixed layer.In the Pacific,negative skewnessexists in the southeastern basin(3°N-5°S,140°-80°W)in WSA?R,but the opposite occurs in WSA. This difference indicates that the skewness of the wind stress anomalies would play a role in the skewness of the temperature anomalies in such a region.However, there is a common negative skewness in the western Pacific(120°-160°E)and in the southeastern tropical Indian Ocean,which confirms the more significant role of ocean dynamics in such regions.

    Fig.4.As in Fig.3,but for mixed layer condition.

    5.Time-mean effect

    5.1 Subsurface

    5.1.1 ENSO/IOD rectification

    Because ofasymmetry,the residualof the positive and negative phases of ENSO/IOD events is eithermuch larger or smaller than zero.The spatial pattern of the residual of ENSO/IOD events is similar to that of the above temperature asymmetry.Given the fact that warm and cold temperature anomalies cannot completely cancel each other out,the long-term mean ENSO/IOD events should influence the climatological mean state or decadaloscillation,as has also been suggested by some recent studies(Rodger et al.,2004;Yu and Kim,2011;Choiet al.,2012).As shown in Fig.5, although the climatological mean wind stress used in all numericalexperiments was identical,the simulated climatological mean temperature differs between the sensitivity experiments and control run.In the Pacific,due to interannual wind stress forcing,the simulated mean temperature in both WSA and WSA?R is 1℃colder than that of the control run in the equatorial western Pacific,while the warm anomalies show a maximum of around 1℃in the equatorial eastern Pacific.In the Indian Ocean,both WSA and WSA?R show temperature decreases in the upper 80 m in theequatorial eastern basin and temperature increases along the thermocline at 80-160 m in the Indian Ocean.In terms of the subsurface,the large-scale spatial pattern in Fig.5 is similar to that of Fig.3,but the locations of the maximum/minimum values in Fig. 5 are not exactly identical to those in Fig.3.In Fig. 5,the time-mean effect from WSA has a high spatialcorrelation with WSA?R,implying that the mean temperature change in the subsurface is due to the asymmetry of the oceanic response to the anomalous surface forcing.Besides the temperature change,the current and salinity differences between the sensitivity experiments and control run showed identical behaviors(figure omitted).Furthermore,we conducted an additionalexperiment to confirm the above results. The experiment was designed to force the ocean model using idealized interannual wind stress anomalies but without any asymmetry.The time-mean temperature differences between the simulation results from this experiment and those from the control run were consistent with Fig.5.

    Fig.5.Temperature differences between control run and(a,c)WSA and(b,d)WSA?R in equatorial upper Pacific and Indian Oceans(2°N-2°S).

    5.1.2 Role of NDH in ENSO/IOD rectification

    Because the observed wind stress anomalies revealed obvious interannual variations,there existed dramatic interannual variations in the simulated temperatures and currents of the sensitivity experiments; however,these were negligible in the controlrun.Consequently,the climatological mean nonlinear advection terms in the sensitivity runs were different from those in the control run.Liang et al.(2012)and Sun et al.(2014)noted that ENSO rectification of the mean state is induced by the nonlinear heating terms of the heat budget formula.Therefore,we present NDH differences between the sensitivity runs(WSA and WSA?R)and the control run(Fig.6),respectively.The results in the Figs.6a and 6c are similar to those in Figs.6b and 6d for the two basins,indicating that the situation is concurrent with Fig.5. Only a minor difference between WSA and WSA?R can be found,in which the maximum difference is limited to 0.1℃mon?1.In the Pacific,the largest positive NDH occurs at 20-60 m between 100°and 80°W, but the negative NDH mainly exists in the upper 80 m between 120°and 160°E.The warming center is located at 20-60 m between 100°and 80°W in Fig.5, while the cooling center at 120 m is around 160°E. Thus,the mean temperature change in the equatorial eastern Pacific(Fig.5)could be explained by NDH,but the temperature decrease in the equatorial western Pacific could not.In the Indian Ocean,the largest negative NDH is located at 20-80 m between 85°and 98°E,and the largest positive NDH is at 80-120 m between 70°and 80°E.The cooling center is located at 20-80 m between 85°and 98°E in Fig.5, and the warming center is at 80-120 m between 70°

    and 100°E.Therefore,NDH could explain the temperature decrease in the Indian Ocean,as well as in the partial warming region.To sum up,NDH favors the rectified effect of temperature asymmetry on the mean state.Further research(Hua et al.,2015)has shown that the upper ocean temperature and velocity anomalies(T?,u?,v?,w?)respond linearly to wind stress anomalies.Therefore,the significant correlation between the temperature anomaly(T?)and velocity anomalies(u?,v?,w?)(and thus the NDH)retains the same sign even though the wind stress anomalies are reversed.In this way,the impact of the rectification on the mean state remains almost unchanged. We further compared the three components of the NDH with the total NDH(figure omitted)and found that the most important components were the zonal and meridional in the Pacific Ocean,while the key component was the vertical in the Indian Ocean.

    5.1.3 Role of mean advection in ENSO/IOD rectification Given the fact that NDH could not account for the mean state change in some regions,such as the western Pacific,Fig.7 shows the differences in mean advection between the sensitivity runs(WSA and WSA?R)and control run.In the Pacific,the cold advection corresponding with the cold temperature anomalies covers the western region(140°E-160°W)below 100 m.The warm advection is located at 0-20 m between 120°and 80°W,and it is also in line with the warm anomalies near the surface of the Pacific in Fig.5.In the Indian Ocean,the warm advection center is located at 60-120 m between 90°and 100°E,which could well explain the warm anomalies near the eastern boundary of theIndian Ocean in Fig.5.The cold advection below 120 m covers the partial region with cold anomalies in Fig.5.Considering that neither the total NDH nor the mean advection could provide an adequate explanation of the mean state change,the non-linearity of the vertical mixing(R in the temperature formula) should also play a role.The largest temperature differences mainly occurred around the thermocline,where vertical mixing cannot be ignored.Moreover,because any verticalmixing scheme is nonlinear and dependent on the mean state,it is vital to diagnose the contribution of vertical mixing to the heat budget in future studies.Although temperature change is related,to some extent,to vertical mixing,understanding the effects of NDH and mean advection on the mean state is essential.

    5.2 Heat budget in the mixed layer

    Fig.6.Differences in NDH between control run and(a,c)WSA and(b,d)WSA?R in equatorial upper Pacific and Indian Oceans(2°N-2°S).

    The pattern of mean temperature change in the mixed layer from WSA is very similar to that fromWSA?R(Fig.8),illustrating the importance of nonlinear dynamics.As can be seen,the temperature increases in the eastern Pacific(15°S-10°N,120°-80°W) and decreases in the western Pacific(140°-160°E).In the Indian Ocean,there is a warming center located over 0°-10°S,40°-80°E,and there are also some cooling regions such as north of equator,south of 10°S.

    Figure 9 further shows the corresponding differences in the NDH between the sensitivity runs and control run.The NDH contributes much to the cooling region in the western Pacific(140°-160°E)as well as to some warming regions(0°-3°N,120°-80°W).In the Indian Ocean,the NDH shows a negative value in most regions that cover the cooling regions in Fig.8.

    The mean advection differences are shown in Fig. 10.These can explain the temperature increases in most regions of the eastern Pacific(15°S-15°N,120°-80°W)and also the warming center(0°-10°S,40°-80°E)in the Indian Ocean.

    Fig.7.Differences in mean advection between control run and(a,c)WSA and(b,d)WSA R in equatorial upper Pacific and Indian Oceans(2°N-2°S).

    6.Conclusions

    The temperature anomalies associated with ENSO/IOD in the tropical Pacific and Indian Oceans not only show significant asymmetry,but also have a rectification effect on the mean state or decadal climate variations.We conducted several numerical experiments with an OGCM in this study.These experiments were designed(1)to evaluate the model’s ability to reproduce the interannual variations,especially the temperature asymmetry in tropical basins;(2)to help quantify the relative contributions of wind stress forcing and ocean dynamical processes to the asymmetry of temperature anomalies;and(3)to understand the rectification effect of the interannual variations on the mean state and its physical mechanism.

    Fig.8.As in Fig.5,but for mixed layer condition.

    When forced by the observed wind stress with interannual anomalies,the model well reproduced thelarge-scale patterns of the mean states in these two basins.Meanwhile,the generalspatialpattern of temperature asymmetry simulated by the ocean model was similar to the reanalysis in the upper 40-160 m in Pacific and 0-140 m in Indian Ocean.The asymmetry in such subsurface regions mainly results from the nonlinear response to the surface wind stress anomalies through internal dynamic processes,rather than the skewness of winds.The climatological mean states simulated by the two numerical experiments were alike,but different from that of the control run, which indicates that the two sensitivity runs,in which the wind stress anomalies were opposite to each other, exerted similar time-mean effects on the mean state. Further analysis showed that NDH plays a crucial role in the warming of the equatorial eastern Pacific below20 m and in the cooling ofthe upper equatorialeastern Indian Ocean.Mean advection contributes much to the cooling in the equatorial western Pacific below 100 m as well as to the warming at the eastern boundary of the Indian Ocean.

    Fig.9.As in Fig.6,but for most important components,i.e.,(a,b)zonal NDH for Pacific and(c,d)vertical NDH for Indian Ocean,in mixed layer condition.Units are℃month?1.

    There are some interannual variation modes in SST anomalies in the tropical Atlantic(Servain,1991; M′elice and Servain,2003),but previous studies have rarely addressed the temperature asymmetry in this basin.Besides numerical experiments in the Pacific and Indian oceans,we conducted similar numericalexperiments in the tropical Atlantic(30°N-30°S,80°W-20°E)in this study.Comparison of these simulations indicated that the asymmetric responses of the temperature and velocity anomalies to the wind stress anomaly are highly universal in the tropical oceans; temperature anomalies are less asymmetric in the tropical Atlantic than they are in the tropical Pacificand Indian oceans,and the corresponding rectification of the mean state is also much weaker.

    Fig.10.As in Fig.7,but for most important components,i.e.,vertical advection for(a,b)Pacific and(c,d)Indian Oceans,in mixed layer condition.Units are℃month?1.

    REFERENCES

    An,S.I.,and Jin Feifei,2001:Collective role of thermocline and zonal advective feedbacks in the ENSO mode.J.Climate,14,3421-3432.

    An,S.I.,and Jin Feifei,2004:Nonlinearity and asymmetry of ENSO.J.Climate,17,2399-2412.

    An,S.I.,Y.G.Ham,J.S.Kug,et al.,2005:El Ni?no-La Ni?na asymmetry in the coupled model intercomparison project simulations.J.Climate,18,2617-2627. Battisti,D.S.,1988:Dynamics and thermodynamics of a warming event in a coupled tropical atmosphereocean model.J.Atmos.Sci.,45,2889-2919.

    Burgers,G.,and D.B.Stephenson,1999:The“normality”of El Ni?no.Geophys.Res.Lett.,26, 1027-1030,doi:10.1029/1999GL900161.

    Cai Wenju,H.H.Hendon,and G.Meyers,2005:Indian Ocean dipolelike variability in the CSIRO Mark 3 coupled climate model.J.Climate,18,1449-1468.

    Cai Wenju and Qiu Yun,2013:An observation-based assessment of nonlinear feedback processes associated with the Indian Ocean dipole.J.Climate,26, 2880-2890.

    Cai Yi,Wang Zhanggui,Yu Zhouwen,et al.,2003:The EOF analysis of temperature and zonal flow in the equatorial Pacific Ocean and the study of the El Ni?no forecasting.Acta Oceanol.Sinica,25,12-18. (in Chinese)

    Cai Yi,Li Hai,and Zhang Renhe,2008:A study on the relationship between ENSO and tropical Indian Ocean temperature.Acta Meteor.Sinica,66,120-124.(in Chinese)

    Canuto,V.M.,and M.S.Dubovikov,2005:Modeling mesoscale eddies.Ocean Modeling,8,1-30.

    Carton,J.A.,and B.S.Giese,2008:A reanalysis of ocean climate using simple ocean data assimilation (SODA).Mon.Wea.Rev.,136,2999-3017.

    Choi,J.,S.I.An,and S.W.Yeh,2012:Decadal amplitude modulation of two types of ENSO and its relationship with the mean state.Climate Dyn.,38, 2631-2644.

    Choi,J.,S.I.An,S.W.Yeh,et al.,2013:ENSO-like and ENSO-induced tropical Pacific decadal variability in CGCMs.J.Climate,26,1485-1501.

    Fedorov,A.V.,and S.G.Philander,2001:A stability analysis of tropical ocean-atmosphere interactions: Bridging measurements and theory for El Ni?no.J. Climate,14,3086-3101.

    Gent,P.R.,and J.C.McWilliams,1990:Isopycnal mixing in ocean circulation models.J.Phys.Oceanogr.,20,150-155.

    Hua Lijuan,Yu Yongqiang,and Yin Baoshu,2010:Numerical modeling of the asymmetry of Indian Ocean dipole and its mechanism.Chinese J.Atmos.Sci.,34,1046-1058.(in Chinese)

    Hua Lijuan,Yu Yongqiang,and Sun Dezheng,2015:A further study of ENSO rectification:Results from an OGCM with a seasonal cycle.J.Climate,28, 1362-1382.

    Hong,C.C.,T.Li,L.Ho,et al.,2008:Asymmetry of the Indian Ocean dipole.Part I:Observational analysis. J.Climate,21,4834-4848.

    Large,W.G.,G.Danabasoglu,S.C.Doney,et al.,1997: Sensitivity to surface forcing and boundary layer mixing in a Global Ocean Model:Annual-mean climatology.J.Phys.Oceanogr.,27,2418-2447.

    Levitus,S.,and T.P.Boyer,1994:World Ocean Atlas 1994 Volume 4:Temperature.NOAA Atlas NESDIS 4.US Department of Commerce,Washington D.C.,1-117.

    Levitus,S.,R.Burgett,and T.P.Boyer,1994:World Ocean Atlas 1994 Volume 3:Salinity.NOAA Atlas NESDIS 3.US Department of Commerce,Washington D.C.,1-99.

    Li Yangchun,Xu Yongfu,Chu Min,et al.,2012:Influences of climate change on the uptake and storage of anthropogenic CO2in the global ocean.Acta Meteor.Sinica,26,304-317,doi:10.1007/s13351-012-0304-z.

    Liang Jin,Yang Xiuqun,and Sun Dezheng,2012:The effect of ENSO events on the tropical Pacific mean climate:Insights from an analytical model.J.Climate,25,7590-7606.

    Liu Hailong,Zhang Xuehong,Yu Yongqiang,et al., 2004a:Manual for LASG/IAP Climate System Ocean Model.Science Press,Beijing,1-108.(in Chinese)

    Liu Hailong,Zhang Xuehong,Li Wei,et al.,2004b:An eddy-permitting oceanic general circulation model and its preliminary evaluation.Adv.Atmos.Sci.,21,675-690.

    Liu Hailong,Lin Pengfei,Yu Yongqiang,et al.,2012: The baseline evaluation of LASG/IAP Climate System Ocean Model(LICOM)Version 2.Acta Meteor. Sinica,26,318-329,doi:10.1007/s13351-012-0305-y.

    M′elice,J.L.,and J.Servain,2003:The tropical Atlantic meridional SST gradient index and its relationships with the SOI,NAO,and southern Ocean.Climate Dyn.,20,447-464.

    Ogata,T.,Xie Shangping,Lan Jian,et al.,2013:Importance of ocean dynamics for the skewness of the Indian Ocean dipole mode.J.Climate,26,2145-2159.

    Ohba,M.,and H.Ueda,2009:Role of nonlinear atmospheric response to SST on the asymmetric transition process of ENSO.J.Climate,22,177-192.

    Okumura,Y.M.,and C.Deser,2010:Asymmetry in the duration of El Ni?no and La Ni?na.J.Climate,23, 5826-5843.

    Rodgers,K.B.,P.Friederichs,and M.Latif,2004:Tropical Pacific decadal variability and its relation to decadalmodulations of ENSO.J.Climate,17,3761-3774.

    Roeske,F.,2001:An Atlas of Surface Flues Based on the ECMWF Reanalysis-A Climatological Dataset to Force Global Ocean General Circulation Models. Report No.23,Hamburg:Max-Planck-Institut f¨ur Meteorologie,1-31.

    Rong Xinyao,Zhang Renhe,T.Li,et al.,2011:Upscale feedback of high-frequency winds to ENSO.Quart. J.R.Meteor.Soc.,137,894-907.

    Sa ji,N.H.,B.N.Goswami,P.N.Vinayachandran,et al., 1999:A dipole mode in the tropical Indian Ocean. Nature,401,360-363.

    Sa ji,N.H.,and T.Yamagata,2003:Structure of SST and surface wind variability during Indian Ocean dipole mode events:COADS observations.J.Climate,16,2735-2751.

    Servain,J.,1991:Simple climatic indices for the tropical Atlantic Ocean and some applications.J.Geophys. Res.,96,15137-15146.

    Su Jingzhi,Zhang Renhe,T.Li,et al.,2010:Causes of the El Ni?no and La Ni?na amplitude asymmetry in the equatorial eastern Pacific.J.Climate,23, 605-617.

    Sun Dezheng,2003:A possible effect of an increase in the warm-pool SST on the magnitude of El Ni?no warming.J.Climate,16,185-205.

    Sun Dezheng and Zhang Tao,2006:A regulatory effect of ENSO on the time-mean thermal stratification of the equatorial upper ocean.Geophys.Res.Lett.,33,L07710.

    Sun Dezheng,Zhang Tao,Sun Yan,et al.,2014:Rectification of El Ni?no-Southern Oscillation into climate anomalies of decadal and longer time scales:Results from forced ocean GCM experiments.J.Climate,27,2545-2561.

    Sun Yan,Sun Dezheng,Wu Lixin,et al.,2012:The western Pacific warm pool and ENSO asymmetry in CMIP3 models.Adv.Atmos.Sci.,30,940-953.

    Uppala,S.M.,P.W.K°allberg,A.J.Simmons,et al., 2005:The ERA-40 re-analysis.Quart.J.Roy.Meteor.Soc.,131,2961-3012.

    Yan Bangliang and Zhang Renhe,2002:The role of atmosphere climate basic state in the formation of westerly over the equatorial Pacific.Acta Oceanol. Sin.,24,39-50.(in Chinese)

    Yu Jinyi and S.T.Kim,2011:Reversed spatial asymmetries between El Ni?no and La Ni?na and their linkage to decadal ENSO modulation in CMIP3 models.J. Climate,24,5423-5434.

    Yu Yongqiang,Zheng Weipeng,Wang Bin,et al.,2011: Versions g1.0 and g1.1 of the LASG/IAP flexible global ocean-atmosphere-land system model.Adv. Atmos.Sci.,28,99-117.

    Zebiak,S.E.,and M.A.Cane,1987:A model El Ni?no-Southern oscillation.Mon.Wea.Rev.,115,2262-2278.

    Zhang Tao,Sun Dezheng,R.Neale,et al.,2009:An evaluation of ENSO asymmetry in the community climate system models:A view from the subsurface. J.Climate,22,5933-5961.

    Zheng Xiaotong,Xie Shangping,G.A.Vecchi,et al., 2010:Indian Ocean dipole response to global warming:Analysis of ocean-atmospheric feedbacks in a coupled model.J.Climate,23,1240-1253.

    Zhong Aihong,H.H.Hendon,and O.Alves,2005:Indian Ocean variability and its association with ENSO in a global coupled model.J.Climate,18,3634-3649.

    Hua Lijuan and Yu Yongqiang,2015:Nonlinear responses of oceanic temperature to wind stress anomalies in the tropical Pacific and Indian Oceans:A study based on numerical experiments with an OGCM.J.Meteor.Res.,29(4),608-626,

    10.1007/s13351-015-4115-x.

    Supported by the“Strategic Priority Research Program Climate Change:Carbon Budget and Relevant Issues”of the Chinese Academy of Sciences(XDA05110302),National(Key)Basic Research and Development(973)Program of China(2013CB956204), and Jiangsu Collaborative Innovation Center for Climate Change.

    ?yyq@lasg.iap.ac.cn.

    ?The Chinese Meteorological Society and Springer-Verlag Berlin Heidelberg 2015

    November 6,2014;in final form June 8,2015)

    亚洲内射少妇av| 亚洲欧美日韩卡通动漫| 久久久午夜欧美精品| 日产精品乱码卡一卡2卡三| 麻豆久久精品国产亚洲av| 欧美日韩一区二区视频在线观看视频在线 | 久久久久久九九精品二区国产| 国产大屁股一区二区在线视频| 亚洲精品乱码久久久v下载方式| 成熟少妇高潮喷水视频| 亚洲久久久久久中文字幕| 国国产精品蜜臀av免费| 国产伦精品一区二区三区视频9| 日日撸夜夜添| 又粗又爽又猛毛片免费看| 婷婷精品国产亚洲av在线| 亚洲自拍偷在线| 亚洲18禁久久av| 丰满人妻一区二区三区视频av| 中文字幕av成人在线电影| 久久精品久久久久久噜噜老黄 | 在线免费观看的www视频| 久久精品久久久久久噜噜老黄 | 在线播放无遮挡| 久久久久久久久中文| 欧美绝顶高潮抽搐喷水| 三级国产精品欧美在线观看| 亚洲精品一区av在线观看| 99久久中文字幕三级久久日本| 舔av片在线| 午夜福利在线观看吧| 欧美3d第一页| 国国产精品蜜臀av免费| 婷婷色综合大香蕉| 日韩三级伦理在线观看| 日韩欧美免费精品| 免费人成视频x8x8入口观看| 亚洲av免费高清在线观看| 51国产日韩欧美| 最近的中文字幕免费完整| 99热全是精品| 一级av片app| 精品人妻一区二区三区麻豆 | 国产白丝娇喘喷水9色精品| 国产白丝娇喘喷水9色精品| 亚洲国产高清在线一区二区三| 欧美不卡视频在线免费观看| 你懂的网址亚洲精品在线观看 | 精品日产1卡2卡| 久久久a久久爽久久v久久| 人人妻人人看人人澡| 精品不卡国产一区二区三区| 亚洲精品日韩av片在线观看| 一个人观看的视频www高清免费观看| 亚洲av电影不卡..在线观看| 欧美区成人在线视频| 久久久成人免费电影| 一个人观看的视频www高清免费观看| 久久久国产成人精品二区| 青春草视频在线免费观看| 国产毛片a区久久久久| ponron亚洲| 婷婷色综合大香蕉| 我的女老师完整版在线观看| 极品教师在线视频| 国产综合懂色| 成人综合一区亚洲| 人人妻人人看人人澡| 日韩,欧美,国产一区二区三区 | 午夜爱爱视频在线播放| 日本撒尿小便嘘嘘汇集6| 国产精品久久电影中文字幕| 成人综合一区亚洲| 99久久九九国产精品国产免费| 神马国产精品三级电影在线观看| 成人午夜高清在线视频| 午夜久久久久精精品| 在线a可以看的网站| 免费电影在线观看免费观看| 可以在线观看的亚洲视频| 精品久久久久久久人妻蜜臀av| 精品人妻偷拍中文字幕| 观看美女的网站| 午夜福利在线观看免费完整高清在 | 五月伊人婷婷丁香| 成人国产麻豆网| 高清日韩中文字幕在线| 激情 狠狠 欧美| 久久久精品94久久精品| 美女xxoo啪啪120秒动态图| 日本免费a在线| 精品一区二区三区视频在线| 日本三级黄在线观看| 一级毛片aaaaaa免费看小| 日韩欧美在线乱码| 黄色一级大片看看| 日本-黄色视频高清免费观看| 久久精品久久久久久噜噜老黄 | 免费看光身美女| h日本视频在线播放| 日韩中字成人| 给我免费播放毛片高清在线观看| 18+在线观看网站| 欧美日韩在线观看h| 狠狠狠狠99中文字幕| 欧美高清性xxxxhd video| 免费av不卡在线播放| 国产私拍福利视频在线观看| 可以在线观看的亚洲视频| 91在线观看av| 精品一区二区免费观看| 国产淫片久久久久久久久| 成人精品一区二区免费| 亚洲成人中文字幕在线播放| 国内精品宾馆在线| 91狼人影院| 久久精品国产自在天天线| 精品久久久久久久末码| 少妇的逼水好多| 自拍偷自拍亚洲精品老妇| 嫩草影院精品99| 久久久国产成人免费| 国内久久婷婷六月综合欲色啪| 日韩成人伦理影院| 一区二区三区高清视频在线| 亚洲真实伦在线观看| 长腿黑丝高跟| 97热精品久久久久久| 日本一二三区视频观看| 欧美色欧美亚洲另类二区| 免费黄网站久久成人精品| 国产成人a区在线观看| 午夜日韩欧美国产| 伦理电影大哥的女人| 精品一区二区三区人妻视频| 久久国产乱子免费精品| 国产精品人妻久久久影院| 一级毛片久久久久久久久女| 好男人在线观看高清免费视频| 又爽又黄无遮挡网站| 免费无遮挡裸体视频| 免费不卡的大黄色大毛片视频在线观看 | 欧洲精品卡2卡3卡4卡5卡区| 免费观看在线日韩| av在线天堂中文字幕| 成人亚洲精品av一区二区| 亚洲人成网站在线播放欧美日韩| 久久精品影院6| 麻豆乱淫一区二区| 免费av不卡在线播放| 天堂√8在线中文| 免费在线观看成人毛片| 国产精品亚洲一级av第二区| 久久综合国产亚洲精品| 男女那种视频在线观看| 人人妻人人澡欧美一区二区| 精品久久久久久久末码| 日日摸夜夜添夜夜添av毛片| 亚洲一级一片aⅴ在线观看| 简卡轻食公司| 亚洲av五月六月丁香网| 国产成人福利小说| 亚洲美女搞黄在线观看 | 午夜激情福利司机影院| 国产亚洲精品久久久久久毛片| 变态另类丝袜制服| 免费黄网站久久成人精品| av在线蜜桃| 乱人视频在线观看| 天堂网av新在线| 亚洲激情五月婷婷啪啪| 变态另类丝袜制服| 老熟妇乱子伦视频在线观看| 精品乱码久久久久久99久播| a级毛片a级免费在线| 久久久久九九精品影院| a级毛色黄片| 欧美最黄视频在线播放免费| 久久久久免费精品人妻一区二区| 日韩一区二区视频免费看| 成年版毛片免费区| 婷婷精品国产亚洲av在线| 精品乱码久久久久久99久播| 一区二区三区四区激情视频 | 一区二区三区四区激情视频 | 特大巨黑吊av在线直播| 国产成人精品久久久久久| or卡值多少钱| 亚洲第一电影网av| 男女啪啪激烈高潮av片| 高清毛片免费观看视频网站| 久久鲁丝午夜福利片| 97超视频在线观看视频| 女的被弄到高潮叫床怎么办| 97碰自拍视频| 国产亚洲精品久久久com| 在线播放无遮挡| 国产成人91sexporn| 51国产日韩欧美| 国产白丝娇喘喷水9色精品| 麻豆国产97在线/欧美| 成熟少妇高潮喷水视频| 神马国产精品三级电影在线观看| 给我免费播放毛片高清在线观看| 亚洲av五月六月丁香网| 久久综合国产亚洲精品| 人妻丰满熟妇av一区二区三区| 亚洲欧美日韩无卡精品| 国产真实乱freesex| 少妇高潮的动态图| 床上黄色一级片| 女的被弄到高潮叫床怎么办| 人人妻,人人澡人人爽秒播| 免费搜索国产男女视频| 久久久久免费精品人妻一区二区| 十八禁网站免费在线| 成人二区视频| 亚洲精品日韩在线中文字幕 | 在线播放国产精品三级| 国产真实伦视频高清在线观看| 午夜福利视频1000在线观看| 亚洲色图av天堂| 欧美一区二区国产精品久久精品| 国产成人a区在线观看| 精华霜和精华液先用哪个| 国产淫片久久久久久久久| 国产精品国产高清国产av| 有码 亚洲区| 亚洲国产精品久久男人天堂| 国产 一区精品| 男人和女人高潮做爰伦理| 成人高潮视频无遮挡免费网站| 欧美精品国产亚洲| 国产aⅴ精品一区二区三区波| 国内少妇人妻偷人精品xxx网站| 亚洲国产精品成人综合色| 可以在线观看毛片的网站| 一a级毛片在线观看| 最近2019中文字幕mv第一页| 精品不卡国产一区二区三区| 亚洲国产精品成人综合色| 精品久久久久久久人妻蜜臀av| 菩萨蛮人人尽说江南好唐韦庄 | 欧美人与善性xxx| 亚洲一区二区三区色噜噜| 亚洲综合色惰| 国产精品久久视频播放| 在线播放国产精品三级| 久久精品国产亚洲av涩爱 | 久久久午夜欧美精品| 变态另类丝袜制服| 好男人在线观看高清免费视频| 日本一本二区三区精品| 中文亚洲av片在线观看爽| 亚洲国产欧洲综合997久久,| 亚洲精品影视一区二区三区av| 国产视频内射| 国产伦一二天堂av在线观看| 人人妻,人人澡人人爽秒播| 偷拍熟女少妇极品色| 成人鲁丝片一二三区免费| 亚洲成人久久性| 亚洲国产精品成人综合色| 欧美高清成人免费视频www| 日韩欧美三级三区| 国产av一区在线观看免费| 日本与韩国留学比较| 毛片女人毛片| 伦精品一区二区三区| 国产不卡一卡二| 少妇猛男粗大的猛烈进出视频 | 国产精品一区二区三区四区久久| 六月丁香七月| 国产在视频线在精品| 欧美三级亚洲精品| 久久婷婷人人爽人人干人人爱| 麻豆国产97在线/欧美| 热99re8久久精品国产| 在线国产一区二区在线| 97人妻精品一区二区三区麻豆| av视频在线观看入口| 观看免费一级毛片| 久久热精品热| 淫秽高清视频在线观看| 免费看av在线观看网站| 大香蕉久久网| 日韩三级伦理在线观看| 国产麻豆成人av免费视频| 精品少妇黑人巨大在线播放 | 22中文网久久字幕| 最近视频中文字幕2019在线8| 蜜桃亚洲精品一区二区三区| 一级黄色大片毛片| 亚洲精品色激情综合| 美女 人体艺术 gogo| 亚洲性久久影院| 午夜a级毛片| 此物有八面人人有两片| 黄色视频,在线免费观看| 国产69精品久久久久777片| 成人高潮视频无遮挡免费网站| 丰满的人妻完整版| 春色校园在线视频观看| 99久久精品一区二区三区| 亚洲美女黄片视频| 久久久久国内视频| 国产精品乱码一区二三区的特点| 国产精品女同一区二区软件| 波多野结衣巨乳人妻| 长腿黑丝高跟| 99热精品在线国产| 国产黄a三级三级三级人| 真实男女啪啪啪动态图| 97超视频在线观看视频| 又爽又黄a免费视频| 男人舔女人下体高潮全视频| 老司机福利观看| 婷婷精品国产亚洲av在线| 久99久视频精品免费| 亚洲第一区二区三区不卡| 中国美女看黄片| 黄色一级大片看看| 天美传媒精品一区二区| 精品少妇黑人巨大在线播放 | 国产亚洲精品久久久久久毛片| 欧美bdsm另类| 欧美成人a在线观看| 久久久a久久爽久久v久久| 亚洲精品成人久久久久久| 欧美性猛交╳xxx乱大交人| 老司机福利观看| 国产精品久久久久久久电影| 亚洲精品国产成人久久av| 色av中文字幕| 精品熟女少妇av免费看| 淫秽高清视频在线观看| 国产精品久久电影中文字幕| 免费在线观看影片大全网站| 少妇裸体淫交视频免费看高清| 在线播放国产精品三级| 91午夜精品亚洲一区二区三区| 日日干狠狠操夜夜爽| 精品人妻一区二区三区麻豆 | 国产精品女同一区二区软件| 久久久久精品国产欧美久久久| 亚洲精品亚洲一区二区| 一个人观看的视频www高清免费观看| 99精品在免费线老司机午夜| 一级a爱片免费观看的视频| 99热这里只有是精品在线观看| 欧美又色又爽又黄视频| 欧美色欧美亚洲另类二区| 国产三级中文精品| 亚洲成a人片在线一区二区| 久久欧美精品欧美久久欧美| 麻豆精品久久久久久蜜桃| 日韩人妻高清精品专区| 韩国av在线不卡| 国产精品久久久久久亚洲av鲁大| 99热这里只有是精品在线观看| 99国产极品粉嫩在线观看| 少妇人妻精品综合一区二区 | 深夜a级毛片| 成人鲁丝片一二三区免费| 久久人人精品亚洲av| 色视频www国产| av天堂在线播放| 午夜激情欧美在线| 国产精品福利在线免费观看| 亚洲一区高清亚洲精品| 18禁黄网站禁片免费观看直播| 蜜桃亚洲精品一区二区三区| avwww免费| 校园春色视频在线观看| 国产精品av视频在线免费观看| 欧美另类亚洲清纯唯美| 九九久久精品国产亚洲av麻豆| 亚洲国产精品久久男人天堂| 国产私拍福利视频在线观看| 亚洲欧美日韩高清在线视频| 国产精品三级大全| 亚洲成av人片在线播放无| 国产一区二区三区在线臀色熟女| 欧美激情国产日韩精品一区| 日日干狠狠操夜夜爽| 欧美激情久久久久久爽电影| 欧美成人一区二区免费高清观看| 成年女人看的毛片在线观看| 久99久视频精品免费| 男女啪啪激烈高潮av片| 日本黄色片子视频| 国产不卡一卡二| 九九爱精品视频在线观看| 九九久久精品国产亚洲av麻豆| 亚洲国产精品国产精品| 五月玫瑰六月丁香| 少妇的逼水好多| 久久久精品欧美日韩精品| 可以在线观看毛片的网站| 国产欧美日韩精品亚洲av| 亚洲国产精品久久男人天堂| 国产一区二区三区av在线 | 女人十人毛片免费观看3o分钟| a级一级毛片免费在线观看| 久久中文看片网| 亚洲人成网站在线播| 人人妻人人澡欧美一区二区| 午夜老司机福利剧场| 亚洲成av人片在线播放无| 久久婷婷人人爽人人干人人爱| 美女 人体艺术 gogo| 亚洲电影在线观看av| 日本免费a在线| 一级毛片aaaaaa免费看小| 国产高清激情床上av| 夜夜夜夜夜久久久久| 久久久久免费精品人妻一区二区| 热99在线观看视频| 精品无人区乱码1区二区| 99久久久亚洲精品蜜臀av| 在线观看午夜福利视频| 国产大屁股一区二区在线视频| 中国美白少妇内射xxxbb| 真实男女啪啪啪动态图| 久久久久免费精品人妻一区二区| 久久精品夜夜夜夜夜久久蜜豆| 国产精品久久久久久久久免| 村上凉子中文字幕在线| 一级黄片播放器| 国产高清视频在线播放一区| 久久午夜福利片| 又黄又爽又免费观看的视频| 舔av片在线| 天堂av国产一区二区熟女人妻| 国产亚洲91精品色在线| 中文字幕人妻熟人妻熟丝袜美| 美女大奶头视频| 亚洲一区高清亚洲精品| 亚洲精品亚洲一区二区| 亚洲电影在线观看av| 亚洲熟妇中文字幕五十中出| 老师上课跳d突然被开到最大视频| 搡女人真爽免费视频火全软件 | 国产aⅴ精品一区二区三区波| 少妇丰满av| 国产69精品久久久久777片| 91久久精品电影网| 久久天躁狠狠躁夜夜2o2o| 黄色一级大片看看| 免费观看人在逋| 午夜日韩欧美国产| 十八禁网站免费在线| 国产在视频线在精品| www日本黄色视频网| 日本黄色片子视频| 男人舔奶头视频| 国产精品一区www在线观看| 身体一侧抽搐| 日本一二三区视频观看| 在线观看免费视频日本深夜| 日韩av在线大香蕉| 搞女人的毛片| 听说在线观看完整版免费高清| 欧美xxxx性猛交bbbb| 一级毛片aaaaaa免费看小| 久久婷婷人人爽人人干人人爱| 色尼玛亚洲综合影院| 菩萨蛮人人尽说江南好唐韦庄 | 亚洲性久久影院| 中国国产av一级| 国产精品野战在线观看| 国产精品久久久久久精品电影| 成熟少妇高潮喷水视频| 国产精品无大码| 晚上一个人看的免费电影| 久久欧美精品欧美久久欧美| 天天躁日日操中文字幕| videossex国产| 免费一级毛片在线播放高清视频| 麻豆精品久久久久久蜜桃| 两个人的视频大全免费| 卡戴珊不雅视频在线播放| 人人妻人人看人人澡| 国产一区亚洲一区在线观看| 老司机福利观看| 国产av在哪里看| 熟女人妻精品中文字幕| 99久久久亚洲精品蜜臀av| 男人和女人高潮做爰伦理| 国产亚洲精品久久久久久毛片| 18禁在线无遮挡免费观看视频 | 精品乱码久久久久久99久播| 中文在线观看免费www的网站| 久久久成人免费电影| 国产91av在线免费观看| 国产一区二区在线观看日韩| 蜜桃亚洲精品一区二区三区| 午夜影院日韩av| 亚洲一区高清亚洲精品| 午夜福利18| 女人被狂操c到高潮| 深夜a级毛片| 成年女人毛片免费观看观看9| 亚洲美女黄片视频| 亚洲一区高清亚洲精品| 欧美日韩在线观看h| 久久久国产成人精品二区| 国产精品无大码| 成人二区视频| 黄片wwwwww| 成熟少妇高潮喷水视频| 欧美精品国产亚洲| 免费看光身美女| 日本欧美国产在线视频| 少妇人妻一区二区三区视频| 国产又黄又爽又无遮挡在线| 日韩三级伦理在线观看| 日韩大尺度精品在线看网址| 日韩精品有码人妻一区| 国产一区二区在线观看日韩| 成年免费大片在线观看| 在线天堂最新版资源| 国产成人影院久久av| av福利片在线观看| 午夜精品在线福利| 久久精品人妻少妇| 美女高潮的动态| 亚洲自拍偷在线| 校园春色视频在线观看| 99久久中文字幕三级久久日本| 又爽又黄a免费视频| 97超视频在线观看视频| 免费一级毛片在线播放高清视频| 舔av片在线| 久久久久久久久中文| 嫩草影院新地址| 网址你懂的国产日韩在线| 日本一二三区视频观看| 久久久久久九九精品二区国产| 秋霞在线观看毛片| 日产精品乱码卡一卡2卡三| 亚洲精品456在线播放app| 搡女人真爽免费视频火全软件 | 内射极品少妇av片p| 一个人观看的视频www高清免费观看| 三级男女做爰猛烈吃奶摸视频| 午夜福利视频1000在线观看| 日本成人三级电影网站| 日本欧美国产在线视频| 嫩草影视91久久| 91精品国产九色| 久久精品国产亚洲av香蕉五月| 真实男女啪啪啪动态图| 少妇的逼好多水| 在线国产一区二区在线| 一级毛片久久久久久久久女| 国产白丝娇喘喷水9色精品| 在线观看午夜福利视频| 亚洲一区二区三区色噜噜| 日韩 亚洲 欧美在线| 2021天堂中文幕一二区在线观| 免费看日本二区| 亚洲美女视频黄频| 国产高清三级在线| 国产一区二区激情短视频| 久久精品国产亚洲av涩爱 | 精品免费久久久久久久清纯| 精品福利观看| 国内精品宾馆在线| 又粗又爽又猛毛片免费看| 精品无人区乱码1区二区| 美女大奶头视频| 久久久久久久久久久丰满| 国产高清视频在线观看网站| 99热精品在线国产| 九九热线精品视视频播放| 亚洲性久久影院| 偷拍熟女少妇极品色| 国产探花极品一区二区| 校园人妻丝袜中文字幕| 国产白丝娇喘喷水9色精品| 久久久久九九精品影院| 国产精品一区二区性色av| 亚洲人与动物交配视频| 国产精品美女特级片免费视频播放器| 亚洲中文字幕日韩| 一本久久中文字幕| 人妻夜夜爽99麻豆av| 深夜a级毛片| 久久午夜亚洲精品久久| 大型黄色视频在线免费观看| 三级毛片av免费| av中文乱码字幕在线| 成年女人毛片免费观看观看9| 婷婷六月久久综合丁香| 欧美一级a爱片免费观看看| 久久久精品94久久精品| 内射极品少妇av片p| 精品久久久久久久末码| 日本黄大片高清| 少妇人妻精品综合一区二区 | 中文字幕久久专区| 国产精品精品国产色婷婷| 国内精品一区二区在线观看| 2021天堂中文幕一二区在线观| 在线免费观看不下载黄p国产| 国产女主播在线喷水免费视频网站 | 久久这里只有精品中国| 啦啦啦韩国在线观看视频| 听说在线观看完整版免费高清| 亚洲色图av天堂| 毛片女人毛片| 亚洲激情五月婷婷啪啪| 亚洲欧美精品综合久久99| 我的女老师完整版在线观看|