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

    Northern Hemisphere Sudden Stratospheric Warming and Its Downward Impact in Four Chinese CMIP6 Models

    2021-02-26 08:22:14JianRAOSimingLIUandYuanhaoCHEN
    Advances in Atmospheric Sciences 2021年2期

    Jian RAO, Siming LIU, and Yuanhao CHEN

    1Key Laboratory of Meteorological Disaster, Ministry of Education (KLME) / Joint International Research Laboratory of Climate and Environment Change (ILCEC) / Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters (CIC-FEMD), Nanjing University of Information Science and Technology, Nanjing 210044, China

    2Fredy and Nadine Herrmann Institute of Earth Sciences, The Hebrew University of Jerusalem,Edmond J. Safra Campus, Givat Ram Jerusalem 91904, Israel

    3Department of the Geophysical Sciences, University of Chicago, Chicago 60637, USA

    ABSTRACT

    Key words:sudden stratospheric warming,CMIP6,surface impact,model simulation

    1.Introduction

    Sudden stratospheric warming (SSW) is one of the most radical phenomena in the climate system, which mainly occurs in the Northern Hemisphere midwinter(Andrews et al., 1987; Charlton and Polvani, 2007; Butler et al., 2015), with only two SSWs occurring in the Southern Hemisphere (September 2002 and 2020) in the satellite record (Newman and Nash, 2005; Rao et al., 2020d; Shen et al., 2020a, b). When SSWs appear, the Arctic stratosphere warms by tens of degrees within several days and the meridional temperature gradient in the subpolar region is reversed in the stratosphere. According to the sign of zonal-mean zonal winds at 10 hPa and 60°N, SSWs are further classified into major and minor events: major SSWs are additionally accompanied by a reversal of zonal winds from westerlies to easterlies, whereas minor SSWs only show a deceleration of westerlies without a direction reversal of zonal winds (Charlton and Polvani, 2007; Hu et al., 2014). On average, minor SSWs have a weaker strength and degree of stratosphere–troposphere coupling, so SSWs here refer to major events unless otherwise specified, although some minor SSWs are also accompanied by a clear downward propagation signal and exert strong influences on the lower-tropospheric climate (Wang and Chen, 2010; Rao et al., 2020d;Shen et al., 2020a, b).

    Weakening, distortion, and displacement of the stratospheric polar vortex during major SSWs project onto a negative stratospheric annular mode, which propagates downward gradually in the following month(s) after onset of SSWs (Baldwin and Dunkerton, 1999; Baldwin et al., 2003;Cai and Ren, 2007; Sigmond et al., 2013). Therefore, stratospheric disturbance associated with SSWs is strongly coupled with the troposphere and serves as a potential source for tropospheric variability and predictability on the subseasonal time scale (e.g., Karpechko et al., 2018; Rao et al., 2018; Domeisen et al., 2020; Taguchi, 2020). For example, tropospheric and near-surface predictability is enhanced following SSWs (Sigmond et al., 2013; Tripathi et al., 2015, 2016). On average, a cold Eurasian continent–warm North American continent pattern is observed before SSWs at 850 hPa, while the two continents are anomalously cold after SSWs (Cao et al., 2019; Liu et al., 2019).

    Based on the morphology of the stratospheric polar vortex, SSWs are classified into vortex displacement and vortex split events (Charlton and Polvani, 2007; Seviour et al.,2013; Liu et al., 2019; Rao et al., 2019). Early studies reported that the surface impact of displacement and split SSWs is similar (Charlton and Polvani, 2007; Mitchell et al.,2013), while other recent studies emphasize their differences (Nakagawa and Yamazaki, 2006; Seviour et al.,2016). Rao et al. (2020a) attribute the different impacts of displacement and split SSWs to their intensities, because on average the composite of displacement events is weaker than split events, although the methods and sample sizes might also cause a debate on whether impacts of the two SSW types are significantly different (Seviour et al., 2013, 2016;O’Callaghan et al., 2014).

    Due to their important role in the climate system, a successful simulation of SSWs acts as an indispensable requisite when we score a model. For example, some models from phase 5 of the Coupled Model Intercomparison Project(CMIP5) can simulate several aspects of SSWs to different degrees of success based on single-model assessments [e.g.,HadGEM2 in Osprey et al. (2013); CESM1-WACCM in Cao et al. (2019)] and multimodel studies [e.g., 21 models in Charlton-Perez et al. (2013); 13 high-top models in Seviour et al. (2016)]. Especially, Charlton-Perez et al.(2013) suggest low-top models underestimate stratospheric variability on interannual and daily time scales. Osprey et al. (2013) compared the SSW frequency in the high-top and low-top configuration of their model and found a better reproduction of the SSW frequency in the high-top version than the low-top version. A suitable vertical resolution in the stratosphere is also necessary for models to capably simulate the stratospheric processes and stratosphere–troposphere coupling (Charlton-Perez et al., 2013; Osprey et al., 2013),although the high horizontal resolution is less important than the nice vertical resolution for models to simulate SSWs and the quasi-biennial oscillation (QBO) (Rao et al.,2020b, c).

    Using a loose definition of SSWs and three criteria based on the leading mode of the extratropical zonal winds,Rao et al. (2015) found the CMIP5 multimodel ensemble,including four Chinese models (BCC-CSM1-1, BCCCSM1-1-m, FGOALS-g2, and FGOALS-s2), underestimates the frequency of weak stratospheric polar vortex events. Such a bias is not exclusive to CMIP5 models, with it being also found in other models (Charlton et al., 2007;Mitchell et al., 2012; Ayarzagüena et al., 2013). However,we still know little about the performance of CMIP6 models in simulating SSWs. Recently, different CMIP6 modeling groups released their Diagnostic, Evaluation and Characterization of Klima (DECK) experiments (Kclima in Greek means “Climate”). One of the DECK experiments is a historical Atmospheric Model Intercomparison Project (AMIP) simulation from 1979–2014 (the other three are piControl,4×CO2, and 1ptCO2) (Eyring et al., 2016). Our interest is not to include all of the CMIP6 models but to focus on four Chinese CMIP6 models that have been widely used in operations and/or some research projects in China (Zhou et al.,2020).

    By using a strict WMO definition of SSWs and an available DECK experiment, this paper assesses the general performance of the four Chinese CMIP6 models in simulating SSWs, including their frequency, evolution, types, downward propagation, and surface impact. The paper is constructed as follows. Following the introduction, section 2 describes the models, experiments, and datasets. The SSW frequency and its seasonal distribution are compared in section 3. Evolutions and downward propagation of the stratospheric signals associated with SSWs follow in section 4.The tropospheric and near-surface responses to SSWs are shown in section 5. Finally, the main findings are summarized and discussed in section 6. It is expected that our assessments can help the four modeling groups locate the main biases of their models in the stratosphere, and thereby continue to improve these models in the future.

    2.Models, datasets and methods

    2.1.Models, experiments and datasets

    Table 1 lists the four CMIP6 models with daily outputs available for AMIP experiments. The four models are BCC-CSM2-MR, FGOALS-f3-L, FGOALS-g3, and NESM3,with their full name extensions and affiliations listed in the second column of Table 1. BCC-CSM2-MR is developed by the National Climate Center, China Meteorological Administration, and has a moderate horizontal resolution (T106,i.e., 320 × 160 grids, longitude × latitude). FGOALS-f3-L is developed by the State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics (LASG), Institute of Atmospheric Physics, Chinese Academy of Sciences, and has a moderate (although “L”denotes “l(fā)ow resolution”) horizontal resolution (i.e., C96,i.e., 382 × 194 grids, longitude × latitude). FGOALS-g3 is also developed by LASG, but the atmospheric component is different and has a low horizontal resolution (i.e., 180 × 80 grids, longitude × latitude). NESM3 is developed by the Earth System Modeling Center, Nanjing University of Information Science and Technology, and has a low horizontal resolution (i.e., T63, 190 × 95 grids, longitude × latitude). For full details, readers are directed to the model descriptions(Cao et al., 2018; He et al., 2019; Wu et al., 2019; Li et al.,2020). All models have a model top around 1–2 hPa, incorporating the mid-to-lower stratosphere, where SSWs happen. In contrast, BCC-CSM2-MR and NESM3 (46 and 47 levels in total; 18 and 19 levels around 100–10 hPa) have a nicer vertical resolution than FGOALS-f3-L and FGOALS-g3 (32 and 26 levels in total; 8 and 7 levels around 100–10 hPa).

    Table 1. The four Chinese CMIP6 models used in this study. One of the DECK experiments, AMIP, is commonly available for the four models. The size of the AMIP runs for each model is listed in the third column, and all ensemble members are analyzed in the composite.D/S in the fifth column represents the ratio of the vortex displacement and split SSWs. The CMIP6 AMIP experiments start from 1979 and end in 2014.

    Because daily data from AMIP experiments were available for all of the four models at the beginning of this study(October 2019), we use the AMIP outputs. BCC-CSM2-MR and FGOALS-f3-L have three ensemble members, while FGOALS-g3 and NESM3 have five ensemble members (see the third column of Table 1). All the AMIP experiments are forced by the same external forcings, but the initial fields are different. All the ensemble members from the four Chinese CMIP6 models are used in our paper. Considering that the CMIP6 AMIP runs are from 1979–2014, the extracted Japanese 55-year Reanalysis (JRA-55) from 1979–2014(JRA-55; Kobayashi et al., 2015) is used as a baseline for model evaluations. The SSW events from different reanalyses show little difference, especially during the satellite era since 1979 (Rao et al., 2015; Butler et al., 2017), so only the JRA-55 reanalysis is shown. Variables used in our paper include zonal and meridional winds, heights, and temperatures at pressure levels. Because NESM3 does not provide daily heights, we also calculate Ertel’s potential vorticity(PV) as a substitute for height in the stratosphere.

    2.2.Methods

    A vortex-centric diagnostic procedure developed by Seviour et al. (2013) is used to classify the SSW type. This method is established based on the geometry of the stratospheric polar vortex using the geopotential height or PV at 10 hPa. Two-dimensional vortex moments day by day are calculated in the procedure. Two parameters are required, including the centroid and aspect ratio of the stratospheric polar vortex represented by an equivalent ellipse (Mitchell et al.,2011; Seviour et al., 2016). Time series of the centroid and aspect ratio of the stratospheric polar vortex are calculated using a two-dimensional moment equation. The absolute and relative vortex moments (denoted by Mand J) of the modified PV (or height) field are extracted in the Cartesian coordinate. The latitude of the vortex centroid and the aspect ratio of polar vortex during each SSW event are saved using two-dimensional moment diagnostics and geopotential heights (or PVs) on isobaric levels (Matthewman et al., 2009). Note that the results from geopotential height and PV are highly correlated (Seviour et al., 2013, 2016).

    Following Seviour et al. (2013, 2016), an SSW is classified into the vortex split group if the aspect ratio of the vortex is above 2.4 for at least seven days. An SSW is classified into the vortex displacement group if the centroid of the vortex is situated equatorward of 66°N for at least seven days. This threshold-based method has been confirmed to present a similar classification of split and displaced vortices as conventional methods (e.g., Charlton and Polvani,2007; Mitchell et al., 2011). To show the feasibility of the threshold-based method, examples of vortex displacement and split SSWs are provided in Fig.1 from JRA-55 and four CMIP6 models. Obviously, for displacement SSWs, the vortex is far biased from the North Pole, resembling a commalike shape (Figs. 1a–e). In contrast, for split SSWs, the vortex breaks into two comparable pieces in models, as observed in the selected sample from the reanalysis (Figs. 1f–1j). Although the PV (value range: 30–50 PVU, ?PV is drawn for an easy comparison with other models) is diagnosed for the vortex parameters in NESM3, the displacement and split are also clearly present as in other models.

    3.How often do SSWs appear in CMIP6 models?

    Fig.1. Examples of the two types of SSWs for (a, f) JRA-55 on 16 February 1981 and 14 March 1988, (b, g) BCC-CSM2-MR on 11 March 2013 and 11 February 1982, (c, h) FGOALSf3-L on 24 March 2013 and 10 March 2014, (d, i) FGOALSg3 on 9 March 2014 and 28 February 1982, and (e, j) NESM3 on 13 February 1981 and 18 February 1980. The left-hand column shows the height or PV at 10 hPa for vortex displacement SSWs, and the right-hand column shows the height or PV at 10 hPa for the vortex split SSWs. All examples in the four Chinese CMIP6 models are selected from the first AMIP run. Note that daily heights are unavailable for NESM3 and Ertel’s PVs is exclusively shown for this model (the PV sign is reversed for an easy comparison with other models;?PV value ranges: [?50, ?30] PVU).

    In the JRA-55 reanalysis, 23 SSWs appear during 1979–2014 (~0.64 events per year; Table 2). However, the models (excluding NESM3) tend to underestimate the SSW frequency: 30 events in 108 years for BCC-CSM2-MR (i.e.,36 years from 1979–2014 in three AMIP runs; similar for other models), 31 events in 108 years for FGAOLS-f3-L, 47 events in 180 years for FGAOLS-g3, and 192 events in 180 years for NESM3 (see the fifth column of Table 1). The SSW frequency is 0.28, 0.29, 0.26, and 1.1 events per year for the four models, respectively. Namely, three models underestimate the SSW frequency by half, and NESM3 nearly doubles the observed SSW frequency. To get an overview of SSWs in the four CMIP6 models, the month-by-month distributions of SSWs are shown in Fig.2. SSWs mainly occurin midwinter (January and February; unfilled bars in Fig.1)in observations. Obviously, most models simulate a climate drift for SSWs, and more SSWs appear in late winter (February and March), and SSWs in NESM3 are nearly uniformly distributed in most wintertime months except February.

    Table 2. Onset dates of SSW events and the corresponding type of the stratospheric polar vortex (D indicates a vortex displacement and S indicates a vortex split) in the JRA-55 reanalysis (1979–2014). The ratio of the vortex displacement and split SSWs is 1.3 (13/10) in JRA-55 during 1979–2014.

    Seasonal distributions of vortex displacement and vortex split SSWs from November–March are shown in Fig.3.As seen in Fig.3a, vortex displacement SSWs are nearly uniformly distributed in December–February, followed by March. This peak in February is successfully simulated in NESM3, although SSWs occur much more frequently in this model than in JRA-55. Consistent with the distribution of SSWs in Fig.2, all the other three models simulate much fewer displacement SSWs, and SSWs are drifted to late winter (February and/or March).

    Fig.2. Seasonal distribution of the total frequency of SSWs(units: number per year) from November to March for JRA-55(hatched bars) and CMIP6 models (bars in gray shades).

    Fig.3. Seasonal distribution of the frequency (units: events per year) of (a) vortex displacement SSWs and (b) vortex split SSWs in each wintertime month in the JRA-55 reanalysis during 1979–2014 and AMIP runs during 1979–2014 from four Chinese CMIP6 models.

    A stronger seasonality of split SSWs than displacement SSWs is observed for JRA-55, comparing the unfilled bars in Figs. 3a and b. More split SSWs appear in midwinter (January–February) in observations, and far fewer are observed in other wintertime months. Such a seasonality of SSWs observed in JRA-55 is drifted one month later to February–March for most models except BCC-CSM2-MR. Such a climate drift can be tracked to the seasonal evolution of the stratospheric polar vortex, which tends to get strongest in February (January) in models (reanalyses) [Fig.5 in Rao et al. (2015)]. Compared with the three other models,BCC-CSM2-MR is the only one of the four Chinese CMIP6 models that simulates a stratospheric QBO (Rao et al.,2020b, c), which might also affect SSWs.

    In addition to their contrasting seasonal distributions for both types of SSWs in JRA-55, the difference can also be identified for their intensities. We use the warming anomalies in the stratospheric polar cap to denote the intensity of SSWs. The composite strength of SSWs in each month for each type is shown in Fig.4. To reverse the polar night jet that usually reaches climatological maxima in midwinter(Rao et al., 2015), the polar vortex anomalies are expected to be stronger for midwinter SSWs than events in November and March. This expectation is observed in JRA-55 (<15 K in November and March versus > 20 K in midwinter)and simulated in almost all models for both displacement and split SSWs. Although the SSW frequency is not satisfactorily simulated by most models, the contrast in strength between displacement and split is simulated by models to different degrees of success. Specifically, on average, the strength of split SSWs is larger than displacement SSWs in JRA-55, which is simulated in some models (especially in BCC-CSM2-MR and NESM3).

    4.Stratosphere–troposphere coupling during SSWs

    During SSWs, the weak stratospheric polar vortex marks a phase of weak height (pressure) contrast between the midlatitudes and Arctic, which usually corresponds to a negative Northern Annular Mode (NAM) (Baldwin et al.,2003; Ren and Cai, 2007; Rao et al., 2020a). The negative NAM signals associated with SSWs propagate downward gradually with stratospheric anomalies leading tropospheric anomalies and exhibit strong stratosphere–troposphere coupling. The composite temporal evolution of zonal-mean zonal wind anomalies over the extratropical latitude band(55°–75°N) is shown for two types of SSWs from the reanalysis and four models. As shown in the top row of Fig.5,the easterly anomaly (or the westerly deceleration) develops more than a week before the onset date of SSWs, and the response of near-surface easterlies is most significant around day 20 for both displacement and split SSWs in JRA-55. The easterly anomaly in the stratosphere is stronger most of the time during day 0–60 for split SSWs than displacement in JRA-55 (maximum easterly anomaly magnitude: 25 versus 30 m s; Figs. 5a, f and k).

    The stronger stratospheric anomalies for split SSWs are simulated in three models (i.e., BCC-CSM2-MR, FGOALSf3-L, and NESM3), consistent with the temperature anomalies in Fig.4. It is easy to conclude that the split SSWs have a stronger impact on the troposphere than displacement SSWs, but Rao et al. (2020a) argue that the SSW strength is more important than the vortex morphology in inducing a low-level negative NAM response, which is beyond the scope of this study. It is also noticed that the composite intensity of displacement SSWs from FGOALS-g3 is stronger than split SSWs (maximum easterly anomaly magnitude: 20 versus 25 m s; Figs. 5d, i and n). Although the composite difference between displacement and split SSWs is not significant most of the time after the SSW onset, the stronger easterlies for split SSWs are consistently simulated in the other three models (maximum easterly anomaly magnitude for displacement and split SSWs: 30 versus 35 m sfor BCC-CSM2-MR, 15 versus 20 m sfor FGOALS-f3-L,and 10 versus 15 m sfor NESM3) as observed. The nearsurface response due to the downward propagation of easterly anomalies is also simulated.

    Fig.4. Composite area-weighted polar (60°–90°N) temperature anomaly (units: K) at 10 hPa, ±5 days from the onset date of (a) vortex displacement SSWs and (b) vortex split SSWs in each wintertime month for the JRA-55 reanalysis and four Chinese CMIP6 models. The error bar shows uncertainty.

    Fig.5. Composite pressure–time evolution of the zonal mean zonal wind anomalies area-averaged over 55°–75°N (shading;units: m s?1) from day ?20 to day 60 relative to the onset date for (a–e) vortex displacement SSWs and (f–j) vortex split SSWs for (top row) the JRA-55 reanalysis during 1979–2014, and (second–last rows) four Chinese CMIP6 models during 1979–2014. The last column (k–o) shows the difference of vortex split minus displacement SSWs in each dataset. Black contours mark the composite zonal wind anomalies/differences at the 95% confidence level according to the Student’s t-test.

    The evolutions of the polar cap temperature anomalies during SSWs are also compared for both SSW types in JRA-55 and the models (not shown). The temperature anomalies are observed to be larger for split SSWs than displacement SSWs in JRA-55 (maximum positive anomalies: ~8 versus ~10 K) and simulated to be so in BCC-CSM2-MR(~18 versus ~22 K), FGAOLS-f3-L (~8 versus ~10 K), and NESM3 (~6 versus ~7 K). Warm anomalies can propagate downward to the upper troposphere, but easterly anomalies can reach the near surface following onset of SSWs.

    The stronger downward propagation and strength for split SSWs can be tracked to the stronger wave activities in the (upper) troposphere. As the eddy heat flux is proportional to the vertical component of EP flux, evolutions of eddy heat flux anomalies at 100 hPa in the 45°–75°N latitude band are shown in Fig.6 for displacement and split SSWs, respectively. As seen from JRA-55 (green curves),although the total eddy heat flux anomalies for split SSWs are comparable to displacement SSWs from day ?20 to day 0 (comparable maximum values: ~20 K m s; Figs. 6a and b), the eddy heat flux anomalies by the wavenumber-1 for the vortex displacement SSWs are larger than those for the vortex split SSWs (maximum values: ~20 K versus ~15 K m s,significantly different at the 95% confidence level). This means that within three weeks before onset of the vortex displacement SSWs, the upward propagation of planetary waves (especially wavenumber-1) strengthens (Figs. 6c and d). Similarly, the increase in the eddy heat flux by wavenumber-2 is larger for vortex split SSWs than for vortex displacement SSWs (maximum values: ~20 versus ~15 K m s;Figs. 6e and f). After the onset of SSWs, the negative eddy heat flux anomalies develop. Namely, upward propagation of waves is prohibited after onset of SSWs due to the development of easterlies in the stratosphere. Compared with displacement SSWs, the longer time of eddy heat flux anomalies above zero (rather than the peak maxima) before day 0 in JRA-55 means more extra accumulation of energy in the stratospheric Arctic for split SSWs. This is successfully simulated by BCC-CSM2-MR, FGAOLS-f3-L, and NESM3.

    Fig.6. Temporal evolution of eddy heat flux anomalies at 100 hPa area-averaged in the 45°–75°N latitude band(units: K m s?1) from day ?40 to day 40 with respect to the onset date of (left-hand column) vortex displacement SSWs and (right-hand column) vortex split SSWs from JRA-55 (green) and CMIP6 models (red, blue, orange, and purple). The top row shows the eddy heat by total waves, the middle row shows the wavenumber-1 components, and the bottom row shows the wavenumber-2 components. The thickened part of the dashed line denotes the composite at the 95% confidence level according to the t-test.

    5.Impact of SSWs on the lower troposphere

    5.1.Lower-tropospheric temperature response to SSWs

    Fig.7. (a) Composite temperature anomaly distribution (shading, units: K) at 850 hPa in (a) the JRA-55 reanalysis (1979–2014) and(b–e) four Chinese CMIP6 models (1979–2014) during day ?25 to ?15 (first column), day ?15 to ?5 (second column), day ?5 to 5(middle column), day 5 to 15 (fourth column), and day 15 to 25 (fifth column) relative to the onset date of vortex displacement SSWs. Black contours indicate that the composite temperature anomalies are significant at the 95% confidence level according to the t-test. The latitude range is 20°–90°N.

    Previous studies have confirmed that the continental cold surge is modulated by extreme stratospheric events such as SSWs (Yu et al., 2018; Zhang et al., 2019). Composite temperature anomalies at 850 hPa for displacement SSWs are shown in Fig.7 for JRA-55 and the models. In JRA-55, continental cold anomalies develop in North Eurasia before onset of SSWs from day ?25 to day ?5 (cold anomaly center magnitude: 3 K), and North America is covered with warm anomalies (warm anomaly center magnitude: 4 K). After the onset of displacement SSWs, cold anomalies decay in North Eurasia (cold anomaly center magnitude: 1.5 K), and cold anomalies develop in North America (cold anomaly center magnitude: 1.5 K; Fig.7a). The four Chinese CMIP6 models simulate the evolution of the lower-tropospheric temperature response with different degrees of success (Figs. 7b–e). Specifically, BCC-CSM2-MR fails to reproduce the continental-scale cold anomalies in Eurasia before the onset of displacement SSWs, and North America is covered by large cold anomalies (albeit insignificant for most parts of the anomalies) in this model(Fig.7b). In contrast, the other models generally simulate the cold Eurasian pattern before onset of displacement SSWs (Figs. 7c–e), although the temperature anomalies in FGOALS-f3-L and FGOALS-g3 are not as significant as in JRA-55 and NESM3 due to their different sample sizes. The cold North American pattern is also simulated by all of the four models.

    Similarly, Fig.8 presents composite temperature anomalies at 850 hPa for split SSWs in JRA-55 and the models.A significant cold Eurasia (cold center anomaly magnitude:4 K) and warm North America (warm center anomaly magnitude: 2–4 K) pattern is observed before onset of split SSWs in JRA55, whereas both continents are anomalously cold (3 versus 1.5 K) after onset of split SSWs (Fig.8a).The temperature pattern before day 0 is not well simulated by BCC-CSM2-MR, but the uniform cold pattern in most parts of the Eurasian and North American continents (less significant than the observations) after day 0 (i.e., following onset of split SSWs) is successfully simulated (Fig.8b). As for displacement SSWs, the other three models simulate the cold signals in most parts of Eurasia before and around the onset date of split SSWs, albeit with a low significance level (Figs. 8b–e). However, the cold anomalies after split SSWs in most parts of North America are underestimated.

    Fig.8. (a) Composite temperature anomaly distribution (shading, units: K) at 850 hPa in (a) the JRA-55 reanalysis (1979–2014) and(b–e) four Chinese CMIP6 models (1979–2014) during day ?25 to ?15 (first column), day ?15 to ?5 (second column), day ?5 to 5(middle column), day 5 to 15 (fourth column), and day 15 to 25 (fifth column) relative to the onset date of vortex split SSWs. Black contours indicate that the composite temperature anomalies are significant at the 95% confidence level according to the t-test. The latitude range is 20°–90°N.

    Comparing displacement and split SSWs, the cold patterns after the onset date are different. For displacement SSWs, cold anomalies mainly develop in North America,and the cold anomalies in Eurasia decay, more resembling a wavenumber-1 pattern. For split SSWs, most parts of both the Eurasian and North American continents are covered with cold anomalies, more resembling a wavenumber-2 pattern. Liu et al. (2019) evaluated two high-top models(CESM1-WACCM and CESM2-WACCM) from CMIP5/6,and the simulated low-tropospheric response was similar to that in other CMIP5/6 models (Seviour et al., 2016; Cao et al., 2019; Liu et al., 2019). In short, the surface response to SSWs is more consistently captured by models than the SSW frequency.

    5.2.Tropospheric circulation response to SSWs

    Fig.9. Composite geopotential height anomalies (shading; units: gpm) at 500 hPa in (a) the JRA-55 reanalysis (1979–2014) and(b–d) three CMIP6 models (BCC-CSM2-MR, FGAOLS-f3-L, and FGOALS-g3) during day ?25 to ?15 (first column), day ?15 to?5 (second column), day ?5 to 5 (middle column), day 5 to 15 (fourth column), and day 15 to 25 (fifth column) relative to the onset date of vortex displacement SSWs. Black contours indicate that the composite height anomalies are significant at the 95%confidence level according to the t-test. (e) Composite wind anomalies at 500 hPa for NESM3 due to the unavailability of heights.The latitude range is 20°–90°N. Note that daily heights are unavailable for NESM3 and streamlines are exclusively shown for this model.

    On the one hand, tropospheric wave activities can force the development of SSWs; and on the other hand, SSWs can impact the troposphere after onset by inducing a NAM-like response from the stratosphere to the troposphere (i.e., easterly anomalies in the circumpolar region, Fig.5; or positive polar cap height anomalies; not shown). The two-way coupling is dominated by different processes before and after onset of SSWs. To test the variation of the tropospheric circulation, Fig.9 shows the evolution of height anomalies at 500 hPa in JRA-55 and models during displacement SSWs.As seen from JRA-55, a significant low anomaly center develops over the North Pacific from day ?25 to day ?5, and a significant high anomaly center forms in the subtropical central Pacific in observations (Fig.9a). Such a height anomaly distribution resembles a positive phase of the Pacific–North American (PNA) pattern, with another positive lobe in North Canada and another negative lobe in the eastern U.S.,respectively. The negative height anomalies over the North Pacific extend westward to the coastal region of China,which can be projected onto the negative phase of the western Pacific (WP) pattern (Wallace and Gutzler, 1981). The low anomaly center over the North Pacific has been reported as a precursor of a (following) weakening stratospheric polar vortex on both subseasonal and interannual time scales (Garfinkel et al., 2010; Dai and Tan, 2016; Hu et al.,2017; Rao et al., 2019). Around and after onset of SSWs,the tropospheric circulation anomalies decay and no significant circulation systems are observed (Fig.9a).

    Fig.10. Composite geopotential height anomalies (shading; units: gpm) at 500 hPa in (a) the JRA-55 reanalysis (1979–2014) and(b–d) three CMIP6 models (BCC-CSM2-MR, FGAOLS-f3-L, and FGOALS-g3) during day ?25 to ?15 (first column), day ?15 to?5 (second column), day ?5 to 5 (middle column), day 5 to 15 (fourth column), and day 15 to 25 (fifth column) relative to the onset date of vortex split SSWs. Black contours indicate that the composite height anomalies are significant at the 95% confidence level according to the t-test. (e) Composite wind anomalies at 500 hPa for NESM3 due to the unavailability of heights. The latitude range is 20°–90°N. Note that daily heights are unavailable for NESM3 and streamlines are exclusively shown for this model.

    The models simulate the tropospheric circulation anomaly pattern with different degrees of fidelity. The negative height anomaly center over the North Pacific in observations is biased to the Arctic in BCC-CSM2-MR from day?25 to day ?5, and the high anomaly center over the subtropical central Pacific is located further poleward (Fig.9b). In contrast, the negative height anomaly (i.e., the cyclone anomaly) center over the North Pacific before SSWs is well simulated in the other models (Figs. 9c–e). As the tropospheric wave perturbation propagates upward (denoted by the large eddy heat flux pulse in Fig.6a), the precursor decays and the low center moves westward to North Asia during day ?5 to day 5 in JRA-55 and the models (middle column in Fig.9).Although the height dipole (a high center over Iceland and a low center in the subtropics) in the Atlantic sector is not clearly observed in JRA-55, it is consistently simulated in the four CMIP6 models from day 15 to day 25, resembling the negative phase of the North Atlantic Oscillation (NAO).

    Similarly, the evolution of height anomalies at 500 hPa in JRA-55 and the models during vortex split SSWs is shown in Fig.10. Different from the observed tropospheric height pattern associated with displacement SSWs, the negative height anomaly center is situated further west from the Aleutian Islands, more resembling a WP pattern most of the time from day ?25 to day 5 in JRA-55 (Fig.10a). The circulation anomalies in the Pacific sector gradually decay after onset of SSWs, and a strong NAO is observed after onset of SSWs in observations.

    Models have different skills in simulating tropospheric circulation evolutions during SSWs. The low center develops only during day ?25 to day ?15 in BCC-CSM2-MR(Fig.10b), but decays faster than in JRA-55 and other models. A negative NAO is not clearly simulated after onset of SSWs for this model. In contrast, the WP-like pattern during day ?25 to day 5 is well captured by FGAOLS-f3-L,FGOALS-g3, and NESM3 (Figs. 10c–e). Models tend to underestimate the post-SSW tropospheric response amplitude in the Atlantic sector.

    6.Summary and discussion

    SSW is one of the most radical phenomena in the stratosphere–troposphere coupling system, and its successful simulation is a basic requisite for models emphasizing the role of the stratosphere. Based on a strict WMO SSW definition and a threshold-classification criterion, this paper studies the different statistical characteristics and impacts of vortex displacement and split SSWs. Several aspects are assessed for four Chinese CMIP6 models with daily AMIP outputs available. The main findings are as follows.

    In observations (represented by JRA-55), six or seven SSWs happen in 10 years (0.64 events per year). BCCCSM2-MR, FGOALS-f3-L, and FGOALS-g3 underestimate the SSW frequency nearly by half, while NESM3 doubles the SSW frequency. In addition, most SSWs appear in midwinter (January and February), but one-month climate drift for SSWs is consistently simulated in models(SSWs mainly appear in February and March).

    As shown in observations, displacement SSWs are nearly evenly distributed in December–March, while split SSWs primarily happen in January and February, determining the seasonal distribution of total SSWs. Both displacement and split SSWs tend to appear in late winter in models, so models cannot capture the seasonality of split (and therefore total) SSWs.

    On average, the composite strength of split SSWs is larger than that of displacement SSWs in both the reanalysis and most models, although it does not mean that a split SSW is necessarily stronger than a displacement SSW. The longer pulse of strong positive eddy heat flux before onset of split SSWs explains the stronger stratospheric signals than displacement SSWs for reanalysis and most models.An exception is that the composite of displacement SSWs is stronger than the composite of split SSWs in FGAOLS-g3,although the cause for such a bias is still unknown.

    Displacement SSWs and split SSWs have different impacts on lower troposphere (or the near surface), and the air temperature pattern is also different before and after the onset of SSWs. Before onset of displacement SSWs, cold anomalies cover North Eurasia and warm anomalies cover Arctic Canada in observations. After that, cold anomalies in North Eurasia gradually decay and the cold anomalies are confined to North America. Such a wavenumber-1-like temperature anomaly pattern (cold Eurasia and warm North America) in the Northern Hemisphere before onset of displacement SSWs are well forecasted by FGOALS-f3-L,FGOALS-g3, and NESM3 (the skill in BCC-CSM2-MR is relatively low). The cold anomalies in North America after onset of displacement SSWs are captured by all models.

    Similarly, warm anomalies also form in Arctic Canada before onset of split SSWs, but most parts of the two continents after onset of split SSWs are covered by cold anomalies in midlatitudes as seen from JRA-55. The models have different degrees of fidelity for the temperature anomaly pattern before onset of split SSWs, but the two cold continents are captured by all the models. The two cold continents after onset of split SSWs generally display a wavenumber-2-like temperature anomaly pattern in both the reanalysis and models.

    The tropospheric circulation precursors are different for both types of SSWs in observations. Before displacement SSWs, the North Pacific height anomalies are situated near the Aleutian Islands and the PNA pattern develops to its positive phase in JRA-55. However, before split SSWs, the negative height anomalies in the Pacific sector move westward to the coast of East Asia, and the WP pattern develops to its negative phase. The center of the negative height anomalies in the Pacific sector before onset of SSWs is also sensitive to the SSW type in the models. In observations, a negative NAO is observed only after split SSWs, but it is simulated after both types of SSWs in the models, with the response amplitude underestimated.

    To summarize, many aspects of SSWs can be well simulated by Chinese CMIP6 models, although some biases also exist. Rao et al. (2015) used a loose SSW definition and still found that far fewer SSWs were simulated in CMIP5 models. In addition, it is also noted that the CMIP6 models (i.e.,BCC-CSM2-MR, FGOALS-f3-L, and FGOALS-g3) simulate more SSWs than their CMIP5 versions (i.e., BCCCSM1-1, FGOALS-s2, and FGOALS-g2) (2–3 versus 1–2 SSWs per decade). A comparison of SSW simulation between CMIP5 and CMIP6 models is left for a future study. It is also noted that the significance levels in the models are relatively lower than in the reanalysis, although their composite patterns are very similar. This relatively low significance level is mainly attributed to a small sample size of SSWs in the models due to their underestimation of the SSW frequency and the deficiency of the model’s configuration and coarse (vertical) resolution in capturing the characteristics of the SSW (Charlton-Perez et al., 2013; Osprey et al.,2013; Cai et al., 2017; Rao et al., 2020b, c). Comparing with some previous studies (Cao et al., 2019; Liu et al.,2019), high-top models with a finer vertical resolution in the stratosphere (e.g., CESM1-WACCM, CESM2-WACCM) tend to simulate more SSWs. In our four Chinese CMIP6 models, the atmospheric top is around 1 or 2 hPa, so the middle and lower stratosphere is included in those coupled models. In contrast, it is also found that NESM3 has the highest top of the four models, which stops around 1 hPa. Actually, NESM3 also simulates the most SSWs out of the four models.

    Although the SSW frequency is underestimated in three models, the tropospheric precursors for SSWs are realistically simulated by most models, reflecting a better capture of bottom-up effect in models. The bias in the top-down effect of the stratospheric variability related to the SSW in low-top models might explain the deficiency in the simulation of the near-surface Arctic Oscillation pattern and East Asian winter climate (e.g., Wei et al., 2018; Gong et al., 2019),although the four models have a good simulation in some aspects of the tropospheric variability, including monsoon systems, Madden–Julian Oscillation, rainfall and typhoons (He et al., 2019; Li et al., 2020). However, SSWs are forced by different sources: some are caused by lower-tropospheric perturbation, while more are related to lower-stratospheric wave activities (de la Cámara et al., 2019; White et al.,2019). More efforts are still required for the modeling developers to improve several aspects of those models, including a higher model top to incorporate the whole stratosphere and even the mesosphere, a complete gravity wave parameterization from different sources, a chemical feedback to the stratosphere by adding a chemistry module, and improvements in other model components (e.g., ocean,land, and ice) to produce a better stratospheric response to boundary variations.

    Acknowledgements.

    The authors thank the BCC modeling group, two LASG modeling groups, and the NUIST modeling group for uploading their AMIP experiments to CMIP6. All the CMIP6 AMIP data are collected by ESGF (https://esgf-node.llnl.gov/projects/esgf-llnl/). All CMIP6 data used in this study are publicly available. The JRA-55 reanalysis is provided by the JMA and can be downloaded using FTP after a free registration(https://jra.kishou.go.jp/JRA-55/index_en.html#download). This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDA17010105) and the National Key R&D Program of China (Grant No.2016YFA0602104).

    国产aⅴ精品一区二区三区波| 欧美人与性动交α欧美软件| 日本一区二区免费在线视频| 男女之事视频高清在线观看| 国产91精品成人一区二区三区| 999久久久国产精品视频| 精品久久久久久久久久免费视频 | 久久精品成人免费网站| 欧美成人性av电影在线观看| 少妇裸体淫交视频免费看高清 | 91老司机精品| 夜夜夜夜夜久久久久| 欧美精品一区二区免费开放| 日韩视频一区二区在线观看| 日韩成人在线观看一区二区三区| 大型黄色视频在线免费观看| 欧美成人午夜精品| 亚洲欧美一区二区三区久久| 精品无人区乱码1区二区| 亚洲人成77777在线视频| 在线观看免费视频网站a站| 日韩精品青青久久久久久| 精品国产超薄肉色丝袜足j| 久久香蕉精品热| 亚洲一区二区三区不卡视频| 91老司机精品| 午夜久久久在线观看| 精品国产一区二区久久| 久久久久久久久久久久大奶| av在线天堂中文字幕 | 亚洲精品美女久久av网站| 精品国产一区二区三区四区第35| 亚洲专区中文字幕在线| 超碰97精品在线观看| 免费在线观看日本一区| 色尼玛亚洲综合影院| 成人国语在线视频| 亚洲人成电影观看| 视频在线观看一区二区三区| 日日干狠狠操夜夜爽| av在线天堂中文字幕 | 青草久久国产| 久久99一区二区三区| 精品久久久久久久毛片微露脸| 欧美日本中文国产一区发布| 欧美激情极品国产一区二区三区| 亚洲精品一二三| 国产精品香港三级国产av潘金莲| 国产精品影院久久| 午夜免费成人在线视频| 精品一区二区三区视频在线观看免费 | 一级片免费观看大全| 欧美激情高清一区二区三区| 久久久国产精品麻豆| 黄色视频不卡| 亚洲精品美女久久av网站| 88av欧美| 久久天堂一区二区三区四区| 久久精品aⅴ一区二区三区四区| 亚洲国产精品999在线| 亚洲成人免费av在线播放| 18禁观看日本| 久久欧美精品欧美久久欧美| 老汉色∧v一级毛片| 夜夜看夜夜爽夜夜摸 | 精品久久久久久,| 日本五十路高清| 大型黄色视频在线免费观看| 99久久久亚洲精品蜜臀av| 丝袜在线中文字幕| 一级a爱视频在线免费观看| 欧美精品啪啪一区二区三区| 咕卡用的链子| 91九色精品人成在线观看| 国产熟女xx| 久久精品亚洲精品国产色婷小说| 高清欧美精品videossex| 精品久久久精品久久久| 久久久久久久久中文| 日日夜夜操网爽| 搡老乐熟女国产| 黄色毛片三级朝国网站| 日日夜夜操网爽| 久久 成人 亚洲| 夜夜躁狠狠躁天天躁| 精品国产乱子伦一区二区三区| 看免费av毛片| 窝窝影院91人妻| 丁香六月欧美| 免费一级毛片在线播放高清视频 | 亚洲免费av在线视频| 国产色视频综合| 一个人观看的视频www高清免费观看 | 亚洲一区二区三区色噜噜 | 好男人电影高清在线观看| 激情在线观看视频在线高清| 日韩精品中文字幕看吧| 亚洲自偷自拍图片 自拍| 午夜91福利影院| 午夜日韩欧美国产| 亚洲午夜理论影院| 精品国产国语对白av| 欧美 亚洲 国产 日韩一| 一本大道久久a久久精品| 亚洲av成人一区二区三| 国产xxxxx性猛交| 村上凉子中文字幕在线| 天堂动漫精品| 18禁美女被吸乳视频| 两性夫妻黄色片| 久久久国产成人精品二区 | 一a级毛片在线观看| 亚洲熟妇熟女久久| а√天堂www在线а√下载| 午夜免费激情av| 精品一区二区三区四区五区乱码| 男女高潮啪啪啪动态图| 看片在线看免费视频| 亚洲专区国产一区二区| 好男人电影高清在线观看| 婷婷六月久久综合丁香| 久久久久久大精品| 国产成人影院久久av| 可以免费在线观看a视频的电影网站| 窝窝影院91人妻| 黄网站色视频无遮挡免费观看| 一区福利在线观看| 大陆偷拍与自拍| 在线观看免费午夜福利视频| 最新美女视频免费是黄的| 亚洲国产精品合色在线| 宅男免费午夜| 免费在线观看黄色视频的| 欧美日韩国产mv在线观看视频| 啦啦啦免费观看视频1| www.www免费av| 窝窝影院91人妻| 久久久国产欧美日韩av| 日韩高清综合在线| 看黄色毛片网站| 久久国产乱子伦精品免费另类| 日韩av在线大香蕉| 免费观看精品视频网站| 久久国产乱子伦精品免费另类| 久久亚洲真实| 久久青草综合色| 日本五十路高清| 亚洲aⅴ乱码一区二区在线播放 | videosex国产| 亚洲国产中文字幕在线视频| 欧美日韩亚洲高清精品| 中文字幕另类日韩欧美亚洲嫩草| 国产片内射在线| 国产高清视频在线播放一区| 日日干狠狠操夜夜爽| 国产精品一区二区在线不卡| 日韩欧美国产一区二区入口| 国产成人精品无人区| 久久国产精品人妻蜜桃| av中文乱码字幕在线| 交换朋友夫妻互换小说| 中文字幕av电影在线播放| 亚洲自拍偷在线| 中国美女看黄片| 丁香六月欧美| 亚洲精品久久成人aⅴ小说| 中文欧美无线码| 看黄色毛片网站| av在线播放免费不卡| 久久九九热精品免费| 精品高清国产在线一区| 日韩 欧美 亚洲 中文字幕| 国产精品久久视频播放| 男男h啪啪无遮挡| 男人舔女人下体高潮全视频| 别揉我奶头~嗯~啊~动态视频| 一本大道久久a久久精品| 1024香蕉在线观看| 国产97色在线日韩免费| 久久久久久免费高清国产稀缺| videosex国产| 久久久久精品国产欧美久久久| 国产欧美日韩一区二区三区在线| 夜夜爽天天搞| 成人精品一区二区免费| 80岁老熟妇乱子伦牲交| 亚洲人成伊人成综合网2020| av免费在线观看网站| 美女高潮到喷水免费观看| 大码成人一级视频| 99久久久亚洲精品蜜臀av| 一夜夜www| 国产日韩一区二区三区精品不卡| 久久久国产欧美日韩av| 真人一进一出gif抽搐免费| 欧美乱码精品一区二区三区| 国产主播在线观看一区二区| 精品无人区乱码1区二区| 欧美精品亚洲一区二区| 午夜福利在线免费观看网站| 欧美大码av| 人人妻,人人澡人人爽秒播| 日韩人妻精品一区2区三区| 免费在线观看日本一区| 美女午夜性视频免费| 成年人黄色毛片网站| 亚洲一区二区三区不卡视频| 一二三四社区在线视频社区8| 精品久久蜜臀av无| 欧美日韩亚洲高清精品| 久久久久久久精品吃奶| 中文字幕高清在线视频| 免费av毛片视频| 三级毛片av免费| 水蜜桃什么品种好| 另类亚洲欧美激情| 电影成人av| 欧美日韩亚洲国产一区二区在线观看| 欧美一区二区精品小视频在线| 亚洲一码二码三码区别大吗| 一区二区三区国产精品乱码| 99热只有精品国产| 日韩一卡2卡3卡4卡2021年| 在线观看免费视频日本深夜| 极品教师在线免费播放| 亚洲成a人片在线一区二区| 欧美精品啪啪一区二区三区| 亚洲视频免费观看视频| 欧美乱妇无乱码| 日日摸夜夜添夜夜添小说| 国产男靠女视频免费网站| 淫秽高清视频在线观看| 久久久精品国产亚洲av高清涩受| 水蜜桃什么品种好| 久久人妻福利社区极品人妻图片| 精品国产乱子伦一区二区三区| 啦啦啦免费观看视频1| av福利片在线| 亚洲欧美日韩另类电影网站| 叶爱在线成人免费视频播放| 亚洲 欧美一区二区三区| 日韩精品免费视频一区二区三区| 日韩一卡2卡3卡4卡2021年| 黄网站色视频无遮挡免费观看| 中文字幕人妻丝袜制服| 日日爽夜夜爽网站| 日韩欧美三级三区| 日韩欧美在线二视频| 精品人妻1区二区| 青草久久国产| 国产精品98久久久久久宅男小说| 亚洲伊人色综图| 这个男人来自地球电影免费观看| 国产精品久久久久成人av| 欧美黑人欧美精品刺激| 久久久国产一区二区| 欧美+亚洲+日韩+国产| 国产精品久久视频播放| 欧美日韩瑟瑟在线播放| 美女高潮到喷水免费观看| 男女高潮啪啪啪动态图| 美女午夜性视频免费| 三级毛片av免费| 亚洲 欧美 日韩 在线 免费| 人成视频在线观看免费观看| 欧美乱码精品一区二区三区| 成人免费观看视频高清| 亚洲精品久久午夜乱码| 亚洲,欧美精品.| 91麻豆av在线| 91成人精品电影| 久久午夜亚洲精品久久| 侵犯人妻中文字幕一二三四区| 91av网站免费观看| 国产1区2区3区精品| a级片在线免费高清观看视频| 中出人妻视频一区二区| 新久久久久国产一级毛片| tocl精华| 国产视频一区二区在线看| 日日干狠狠操夜夜爽| 看黄色毛片网站| 国产精品一区二区三区四区久久 | 久久久久久久精品吃奶| 国产免费男女视频| 黑丝袜美女国产一区| 免费观看精品视频网站| 国产成年人精品一区二区 | 19禁男女啪啪无遮挡网站| 黄色丝袜av网址大全| 亚洲国产精品999在线| 黄色片一级片一级黄色片| 夜夜看夜夜爽夜夜摸 | 香蕉久久夜色| 成人免费观看视频高清| 国产精品电影一区二区三区| 国产精品一区二区三区四区久久 | 又大又爽又粗| 一二三四在线观看免费中文在| 国产亚洲精品久久久久5区| 窝窝影院91人妻| 欧美激情 高清一区二区三区| 熟女少妇亚洲综合色aaa.| 亚洲少妇的诱惑av| 亚洲欧美激情综合另类| 日本免费一区二区三区高清不卡 | 亚洲欧美日韩高清在线视频| 99热只有精品国产| 精品第一国产精品| 两个人免费观看高清视频| 久久久久国产一级毛片高清牌| 99国产极品粉嫩在线观看| 别揉我奶头~嗯~啊~动态视频| 亚洲成人免费电影在线观看| 国产精品久久视频播放| 一级a爱视频在线免费观看| 久久精品国产亚洲av高清一级| 精品国产国语对白av| 日韩三级视频一区二区三区| 我的亚洲天堂| 亚洲av美国av| 亚洲在线自拍视频| 中文字幕人妻丝袜一区二区| 夜夜躁狠狠躁天天躁| 1024香蕉在线观看| 国产精品一区二区在线不卡| 亚洲精品一二三| 色尼玛亚洲综合影院| av在线播放免费不卡| www国产在线视频色| 男女午夜视频在线观看| 欧美黄色淫秽网站| 美女大奶头视频| 日本撒尿小便嘘嘘汇集6| 久久久久久人人人人人| 69精品国产乱码久久久| 欧美一区二区精品小视频在线| 制服诱惑二区| 亚洲一区二区三区色噜噜 | 欧美不卡视频在线免费观看 | 亚洲成人免费电影在线观看| 在线观看免费日韩欧美大片| 在线观看www视频免费| 丝袜人妻中文字幕| 亚洲 欧美 日韩 在线 免费| 制服诱惑二区| 一进一出抽搐gif免费好疼 | 国产一区在线观看成人免费| 成熟少妇高潮喷水视频| av欧美777| 国产成人欧美在线观看| 在线观看一区二区三区激情| av网站在线播放免费| 亚洲第一av免费看| 91麻豆av在线| 日本 av在线| 精品久久久久久成人av| 久久亚洲精品不卡| 伦理电影免费视频| 精品乱码久久久久久99久播| 国产在线精品亚洲第一网站| 人人妻人人爽人人添夜夜欢视频| 怎么达到女性高潮| 亚洲中文av在线| 久久 成人 亚洲| 亚洲精品粉嫩美女一区| 国产1区2区3区精品| 亚洲中文av在线| 国产1区2区3区精品| 亚洲中文av在线| 真人做人爱边吃奶动态| 露出奶头的视频| 在线国产一区二区在线| 国产精品自产拍在线观看55亚洲| 女人爽到高潮嗷嗷叫在线视频| 久久久精品欧美日韩精品| 欧美日韩国产mv在线观看视频| 大陆偷拍与自拍| 精品久久久精品久久久| 久久精品国产亚洲av高清一级| 老司机亚洲免费影院| 欧美成人性av电影在线观看| 国产精品98久久久久久宅男小说| 超色免费av| 日韩中文字幕欧美一区二区| 国产熟女午夜一区二区三区| 999久久久国产精品视频| 国产国语露脸激情在线看| 亚洲全国av大片| 在线观看免费午夜福利视频| 午夜免费成人在线视频| 精品久久蜜臀av无| 精品人妻在线不人妻| 国产高清videossex| 亚洲成人国产一区在线观看| 久久青草综合色| 美女高潮到喷水免费观看| 一区二区三区国产精品乱码| 久久人妻福利社区极品人妻图片| 日韩精品中文字幕看吧| 久久久久国内视频| 亚洲精品在线美女| 国产极品粉嫩免费观看在线| 国产男靠女视频免费网站| 在线av久久热| 久久狼人影院| 最近最新中文字幕大全电影3 | 悠悠久久av| 欧美日韩福利视频一区二区| 日本五十路高清| 性欧美人与动物交配| 国产xxxxx性猛交| 欧美中文综合在线视频| 久久精品国产清高在天天线| 夜夜躁狠狠躁天天躁| 久久香蕉国产精品| 欧美日韩国产mv在线观看视频| 日韩欧美免费精品| 不卡av一区二区三区| 一个人观看的视频www高清免费观看 | 91精品国产国语对白视频| 午夜福利在线免费观看网站| 欧美国产精品va在线观看不卡| 夜夜爽天天搞| 国产无遮挡羞羞视频在线观看| 69精品国产乱码久久久| 最好的美女福利视频网| 亚洲国产中文字幕在线视频| 88av欧美| 色综合婷婷激情| 欧美另类亚洲清纯唯美| 国产高清激情床上av| 天堂影院成人在线观看| 亚洲精品国产精品久久久不卡| 亚洲午夜理论影院| 男女下面插进去视频免费观看| 久热爱精品视频在线9| 一区二区三区精品91| 一进一出抽搐gif免费好疼 | 国产一卡二卡三卡精品| xxxhd国产人妻xxx| 99国产综合亚洲精品| 不卡av一区二区三区| 亚洲国产欧美一区二区综合| 麻豆久久精品国产亚洲av | 久热这里只有精品99| 国产精品一区二区免费欧美| 90打野战视频偷拍视频| 又黄又爽又免费观看的视频| 在线观看免费日韩欧美大片| 波多野结衣高清无吗| 十分钟在线观看高清视频www| 日韩精品免费视频一区二区三区| 人人澡人人妻人| 啦啦啦 在线观看视频| 国产黄a三级三级三级人| 黄色毛片三级朝国网站| 亚洲九九香蕉| 在线观看66精品国产| 欧美日韩视频精品一区| 免费看a级黄色片| 一级片'在线观看视频| a级片在线免费高清观看视频| 国产激情欧美一区二区| 国产亚洲精品久久久久5区| 国产精品98久久久久久宅男小说| 这个男人来自地球电影免费观看| 午夜激情av网站| 国产成人精品无人区| 国产精品久久久久成人av| 免费观看人在逋| av免费在线观看网站| 嫁个100分男人电影在线观看| 老司机靠b影院| 亚洲成a人片在线一区二区| 国产精品香港三级国产av潘金莲| 国产成人av激情在线播放| 免费高清在线观看日韩| 九色亚洲精品在线播放| 欧美日韩一级在线毛片| 久久国产精品男人的天堂亚洲| 久久午夜综合久久蜜桃| 久久久久久大精品| 欧美乱码精品一区二区三区| 麻豆av在线久日| 免费在线观看影片大全网站| 操美女的视频在线观看| av视频免费观看在线观看| 日韩大尺度精品在线看网址 | 色婷婷av一区二区三区视频| 男女做爰动态图高潮gif福利片 | 在线十欧美十亚洲十日本专区| 又大又爽又粗| 神马国产精品三级电影在线观看 | 天堂动漫精品| 午夜福利影视在线免费观看| 99riav亚洲国产免费| 91字幕亚洲| 一级片'在线观看视频| 99精品欧美一区二区三区四区| 久久精品国产综合久久久| 一级片免费观看大全| 女人被躁到高潮嗷嗷叫费观| 欧美成人免费av一区二区三区| 岛国视频午夜一区免费看| 国产亚洲精品综合一区在线观看 | 免费在线观看影片大全网站| 久久草成人影院| 欧美黄色片欧美黄色片| 女人爽到高潮嗷嗷叫在线视频| 亚洲五月天丁香| 亚洲成人久久性| 久久草成人影院| 欧美黄色片欧美黄色片| 成人永久免费在线观看视频| 亚洲五月婷婷丁香| 午夜影院日韩av| 久久午夜综合久久蜜桃| 99国产综合亚洲精品| 黑人巨大精品欧美一区二区mp4| 黄色片一级片一级黄色片| av欧美777| 日日爽夜夜爽网站| 国产精品久久久久久人妻精品电影| 免费在线观看影片大全网站| 欧美一级毛片孕妇| 91精品三级在线观看| 欧美 亚洲 国产 日韩一| 中文欧美无线码| 波多野结衣av一区二区av| 最新美女视频免费是黄的| 黄色视频,在线免费观看| 国产亚洲av高清不卡| 老司机午夜十八禁免费视频| 久久精品国产清高在天天线| 91九色精品人成在线观看| 18禁裸乳无遮挡免费网站照片 | 搡老熟女国产l中国老女人| 国产黄色免费在线视频| 黑人巨大精品欧美一区二区蜜桃| 国产av一区在线观看免费| 淫秽高清视频在线观看| 精品少妇一区二区三区视频日本电影| 日本vs欧美在线观看视频| 日韩欧美一区二区三区在线观看| netflix在线观看网站| 视频在线观看一区二区三区| 亚洲精品一区av在线观看| 搡老乐熟女国产| 俄罗斯特黄特色一大片| 一区二区三区激情视频| 久久精品国产99精品国产亚洲性色 | 在线视频色国产色| 精品一区二区三区av网在线观看| 少妇裸体淫交视频免费看高清 | 久久天躁狠狠躁夜夜2o2o| 97碰自拍视频| 中文字幕最新亚洲高清| 久久人妻福利社区极品人妻图片| 国产无遮挡羞羞视频在线观看| 老司机午夜福利在线观看视频| 国产欧美日韩精品亚洲av| 女性被躁到高潮视频| 免费在线观看黄色视频的| 国产又色又爽无遮挡免费看| 真人做人爱边吃奶动态| 中国美女看黄片| 91成人精品电影| 高潮久久久久久久久久久不卡| 高清黄色对白视频在线免费看| 日韩免费av在线播放| 美女扒开内裤让男人捅视频| 淫秽高清视频在线观看| 亚洲三区欧美一区| 成人国语在线视频| 99在线人妻在线中文字幕| 97人妻天天添夜夜摸| 国产精品秋霞免费鲁丝片| 欧美老熟妇乱子伦牲交| 国产精品乱码一区二三区的特点 | 国产1区2区3区精品| 1024视频免费在线观看| 亚洲自拍偷在线| 免费在线观看日本一区| 日韩人妻精品一区2区三区| 午夜免费成人在线视频| 91老司机精品| 免费搜索国产男女视频| 乱人伦中国视频| 国产在线精品亚洲第一网站| 男女下面插进去视频免费观看| 自线自在国产av| 精品第一国产精品| 亚洲av熟女| 国产深夜福利视频在线观看| 黑人欧美特级aaaaaa片| cao死你这个sao货| 人人妻,人人澡人人爽秒播| 精品久久久久久久久久免费视频 | 搡老熟女国产l中国老女人| 如日韩欧美国产精品一区二区三区| 男人舔女人下体高潮全视频| 亚洲国产欧美网| 香蕉国产在线看| av欧美777| 国产熟女xx| 一级作爱视频免费观看| 他把我摸到了高潮在线观看| 国产免费男女视频| 一进一出抽搐动态| 日韩欧美三级三区| 老熟妇乱子伦视频在线观看| 麻豆国产av国片精品| 国产深夜福利视频在线观看| 亚洲精品国产一区二区精华液|