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

    High-resolution Simulation of an Extreme Heavy Rainfall Event in Shanghai Using the Weather Research and Forecasting Model: Sensitivity to Planetary Boundary Layer Parameterization

    2021-01-05 06:47:06RuiWANGYitingZHUFengxueQIAOXinZhongLIANGHanZHANGandYangDING
    Advances in Atmospheric Sciences 2021年1期

    Rui WANG, Yiting ZHU, Fengxue QIAO*,3, Xin-Zhong LIANG, Han ZHANG, and Yang DING

    1Key Laboratory of Geographic Information Science, Ministry of Education,East China Normal University, Shanghai 200241, China

    2School of Geographic Sciences, East China Normal University, Shanghai 200241, China

    3Institute of Eco-Chongming, Shanghai 200062, China

    4College of Atmospheric Science, Nanjing University of Information Science and Technology, Nanjing 210044, China

    5Department of Atmospheric & Oceanic Science, University of Maryland, College Park, Maryland 20742, USA

    6Earth System Science Interdisciplinary Center, University of Maryland, College Park, Maryland 20742, USA

    7Shanghai Central Meteorological Observatory, Shanghai 200030, China

    ABSTRACT In this study, an extreme rainfall event that occurred on 25 May 2018 over Shanghai and its nearby area was simulated using the Weather Research and Forecasting model, with a focus on the effects of planetary boundary layer (PBL) physics using double nesting with large grid ratios (15:1 and 9:1). The sensitivity of the precipitation forecast was examined through three PBL schemes: the Yonsei University Scheme, the Mellor-Yamada-Nakanishi Niino Level 2.5 (MYNN)scheme, and the Mellor-Yamada-Janjic scheme. The PBL effects on boundary layer structures, convective thermodynamic and large-scale forcings were investigated to explain the model differences in extreme rainfall distributions and hourly variations. The results indicated that in single coarser grids (15 km and 9 km), the extreme rainfall amount was largely underestimated with all three PBL schemes. In the inner 1-km grid, the underestimated intensity was improved; however,using the MYNN scheme for the 1-km grid domain with explicitly resolved convection and nested within the 9-km grid using the Kain-Fritsch cumulus scheme, significant advantages over the other PBL schemes are revealed in predicting the extreme rainfall distribution and the time of primary peak rainfall. MYNN, with the weakest vertical mixing, produced the shallowest and most humid inversion layer with the lowest lifting condensation level, but stronger wind fields and upward motions from the top of the boundary layer to upper levels. These factors all facilitate the development of deep convection and moisture transport for intense precipitation, and result in its most realistic prediction of the primary rainfall peak.

    Key words: PBL parameterization, extreme rainfall, high resolution

    1. Introduction

    Accurate and quantitative forecasting of heavy rainfall events in the warm season over East China is a longstanding challenge for regional climate and/or weather prediction models. Difficulty arises not only from the complicated multi-scale interactions between convective activities and different weather systems (Sun et al., 2010; Fu et al.,2017), but also from the limitations of the model itself, such as the initial state, model resolutions and large uncertainties in representing the related precipitation physical processes(Qiao and Liang, 2017; Srinivas et al., 2018). Among these processes, the planetary boundary layer (PBL) is directly affected by the underlying surface and responds to its changes through vertical mixing and turbulent eddies that transport heat, momentum and moisture to the atmosphere above. The vertical eddy transport in the PBL determines the profiles of low-level moisture, temperature and winds(Beljaars and Viterbo, 1998), which have an impact on convection and precipitation development by surface flux exchanges, thermodynamic instability and low-level forcings (Roebber et al., 2004; Hu et al., 2010; Coniglio et al.,2013; Holtslag et al., 2013; Clark et al., 2015; Cohen et al.,2015). However, large uncertainties exist in representing the PBL processes in current mesoscale models, because different assumptions and various kinds of tuning parameters are required in PBL parameterization schemes to describe the vertical subgrid-scale (SGS) fluxes due to eddy transport (Wyngaard, 2004). To date, different kinds of PBL schemes have been proposed, but no consensus exists yet with respect to their applications in model predictions.

    Previous studies have focused on the development of PBL schemes and their impacts on the simulations of boundary layer structures under various atmospheric conditions(Shin and Hong, 2011; Efstathiou et al., 2013; Dong et al.,2019). Generally, there are two categories for the PBL schemes in the WRF model based on their turbulence closure assumptions. The first category is the nonlocal closure,such as the Yonsei University (YSU) (Hong et al., 2006)and the Asymmetric Convection Model 2 (ACM2) (Pleim,2007, b) schemes. These conventional nonlocal closure PBL schemes determine the nonlocal SGS transport through a mass-flux term (Pleim, 2007a, b) or a gradient-adjustment gamma term (Hong et al., 2006) and adopt multiple vertical levels to simulate the effect of larger eddies in the convective boundary layers (Stensrud, 2007). They are typically characterized by overly vigorous vertical mixing and tend to simulate deep, dry and warm boundary layers (Coniglio et al.,2013; Burlingame et al., 2017). The second is local closure,such as the 1.5-order Bougeault-Lacarrere (Bougeault and Lacarrere, 1989), Mellor-Yamada-Janic (MYJ) (Janji?,1994) and Mellor-Yamada Nakanishi Niino (MYNN)Level 2.5 and 3.0 (Mellor and Yamada, 1982; Nakanishi and Niino, 2006, 2009; Olson et al., 2019) schemes. These are generally turbulent kinetic energy (TKE) closure schemes and adopt only adjacent vertical levels to estimate the turbulent fluxes through a full range of atmospheric turbulent regimes. These local closure schemes, such as MYNN and MYJ, more likely simulate shallow, cool and moist boundary layers due to relatively weak vertical mixing (Burlingame et al., 2017). However, controversy still exists regarding their application under convectively unstable conditions.For instance, Shin and Hong (2011) concluded that nonlocal closure PBL schemes such as YSU and ACM2 are favorable in unstable conditions, while Burlingame et al. (2017)suggested that local closure PBL schemes such as MYJ and the Quasi-normal Scale Elimination (QNSE) scheme (Sukoriansky et al., 2005) show advantages over nonlocal closure schemes in detecting convection initiation.

    Many studies have likewise demonstrated the model sensitivity to the choice of PBL scheme in the simulation of heavy precipitation related to tropical cyclones (Dong et al.,2019) and severe midlatitude convections (Shin and Hong,2011; Efstathiou et al., 2013). For instance, Wang et al.(2013) reported that the YSU scheme predicted the intensity, track and associated precipitation more realistically for a weak typhoon [Muifa (2011)] than the MRF (Medium-Range Forecast) scheme. However, Liu et al. (2017) argued that the MYJ scheme generated a better simulation of the rapid intensification process of Hurricane Katrina (2005) offshore, before its landfall, than the YSU scheme. Dong et al.(2019) showed that the YSU and MYNN2.5 schemes outperformed the MYJ and QNSE schemes in reproducing the precipitation after the landfall of Typhoon Fitow (2013). For severe convection and heavy precipitation in the midlatitudes, Efstathiou et al. (2013) suggested that the YSU scheme produces enhanced vertical mixing and moisture transport, and more realistically simulates the strong convection and intense precipitation over the Chalkidiki peninsula in northern Greece. However, Srinivas et al. (2018) examined an extreme heavy rainfall event in December over India and found that, in the 1-km cloud-resolving grid, the MYNN2.5 scheme showed advantages over the YSU and MYJ schemes in simulating the strong convection and associated heavy rainfall. Therefore, large uncertainty still exists with respect to the effects of PBL schemes on heavy rainfall predictions. In particular, as the model resolution is increased to explicitly represent cumulus convection, the effects of PBL schemes become as important as those of the microphysics parameterizations (Braun and Tao, 2000; Li and Pu, 2008).Therefore, this study attempts to explore the sensitivity of PBL schemes in high-resolution WRF simulations based on an extreme rainfall event over Shanghai and its nearby area on 25 May 2018. This event involved the interaction of multi-scale systems from the synoptic circulation, mesoscale low-level jet (LLJ) and low-level shear line to local convective activities, which provides a great test bed for examining the boundary layer processes and their interaction with the large-scale environment and mesoscale convections.

    The objectives of this study are twofold: (1) to investigate the different sensitivities of PBL schemes in extreme rainfall simulations over the mesoscale (15-, 9-km) and high-resolution (1-km) grids with different means of treating convection; and (2) to identify the model differences and explore the underlying mechanisms that affect the interaction of boundary layer processes with the large-scale forcings and convection development.

    The rest of the paper is organized as follows: Section 2 briefly describes the extreme rainfall event in Shanghai. Section 3 describes the design of the WRF model experiments and their physics settings, as well as the three PBL parameterization schemes (YSU, MYJ, MYNN2.5) used in this study.Section 4 compares the spatiotemporal characteristics of heavy rainfall simulations using the different PBLs over varying grid spacings. Section 5 focuses on exploring the physical mechanisms for model differences using the three PBL schemes in the inner 1-km grid simulations. Section 6 summarizes our conclusions.

    2. Observed synoptic conditions and mesoscale features

    Figure 1 presents the double nesting grids for the WRF experiments and the geographic distribution of daily rainfall observations from the Automated Weather Stations(AWS) network on 25 May 2018. There are 5207 sites distributed in the inner domain, D2. The daily rainfall was accumulated from 0000 LST to 2300 LST (LST denotes Local Standard Time, which is 8 h earlier than UTC). Of the 5207 sites,379 stations registered a daily rainfall accumulation exceeding 100 mm. Among them, daily rainfall for two stations[Chongming (173.8 mm) and Baoshan (132.6 mm) in Shanghai] even broke the 50-year record for May, as shown in Fig. 1c. The rainband covered by these 379 AWS sites is the focus of our study and referred to below as the core rainband, which can roughly be divided into two parts (east and west) by the longitudinal line of 120°E (outlined by the red boxes in Fig. 1b).

    To determine the dividing line, we fully examined the observed hourly variations of extreme rainfall over the stations within the core rainband. As shown in Fig. 2, the coastal stations east of 120°E basically had three rainfall peaks, whereas the stations in the western inland area had two peaks without the early-morning secondary one. The core rainband had three rainfall peaks, with the primary peak at 1200 LST and two secondary peaks at 0700 LST and 1800 LST, respectively. The primary and late-afternoon peaks were contributed by both the east and west parts, while the early-morning peak was primarily attributable to the east part. This implies that the extreme rainfall in these two subregions was caused by different precipitation systems or mechanisms, which will be further analyzed below.

    Fig. 1. (a) WRF experiment domain configurations. (b) Geographic distribution of the observed 24-h accumulated rainfall amounts on 25 May 2018 from AWSs. (c) Location of Shanghai.

    Fig. 2. Hourly variation of AWS-observed precipitation on 25 May averaged over the core rainband (AVE 379) and averaged over the two typical regions in Fig. 1b: East region (East 227)and West region (West 152) (units: mm h-1).

    Figure 3 presents the evolution of three-hourly AWS accumulated rainfall and wind fields, moisture convergence,as well as geopotential heights at 850 hPa from the National Centers for Environmental Prediction (NCEP) Final (FNL)Operational Global Analysis data at 1° grids at 0800 LST,1400 LST and 2000 LST 25 May, and 0200 LST 26 May.Hereafter, the LLJ is defined by wind speeds greater than 12 m sat 850 hPa. At 0800 LST, the southwesterly LLJ was strong and transported a great amount of water vapor to the Anhui-Shanghai region, producing intense moisture convergence and a heavy rainfall center with the presence of a low-level warm shear line along the east coast region. A warm shear line is a discontinuous line between southerly and easterly winds with the cyclonic shear in the lower troposphere (Yan and Yao, 2019). This explains the early-morning rainfall peak in the east rainband. Thereafter, the LLJ gradually weakened due to its diurnal variation (Bonner and Paegle, 1970), but the moisture convergence and uplifting motion associated with the low-level warm shear lines were still conducive to the development of mesoscale convection in the east coastal region. At the same time, a shortwave trough at 500 hPa was propagating eastward (figure not shown), and the low-pressure vortex at 850 hPa ahead of the trough began to affect the west rainband. The strong convergence at the back of the low vortex intensified the uplifting motion and mesoscale convection, resulting in the primary rainfall peak at 1200 LST for the west rainband. At 1400 LST, the LLJ and associated low-level moisture convergence decreased, and then, up until 2000 LST, the nocturnal LLJ intensity began to increase but its northern terminus moved south, and the moisture convergence and rainfall at the back of low vortex were intensified. Therefore,the diurnal variation of the LLJ and associated moisture advection contributed more to the secondary peaks in the early morning (0700 LST) and late afternoon (1800 LST),and the strong uplifting motion caused by the low-level warm shear line (east rainband) and the convergence at the back of the low 850-hPa vortex (west rainband) made significant contributions to the development and maintenance of mesoscale convection and heavy precipitation at 1200 LST.

    Fig. 3. Geopotential height (black solid contours; units: gpm), wind field (black vectors), LLJ(wind speed ≥12 m s-1; thick blue vectors) at 850 hPa, and moisture convergence [red solid contours; units: 10-7g (cm2 hPa s)-1] at 850 hPa from the FNL data and three-hourly AWS accumulated rainfall (colored dots; units: mm), at (a) 0800 LST, (b) 1400 LST and (c) 2000 LST 25 May 2018, and (d) 0200 LST 26 May 2018.

    3. Model configuration and experimental designs

    The model used in this study is WRF, version 3.9.1.1(Skamarock et al., 2008). All domains used one-way nesting without feedback and the same vertical discretization of 51 levels, including 11 layers within the boundary layer(below 1500 m). The model was driven by the NCEP FNL Operational Global Analysis Data (1° × 1° grid) as initial and boundary conditions for the outer domains. The 30-arcsecond USGS terrain data and 21-category, 2-m resolution MODIS land use/cover data were prescribed for static surface conditions. The model physics included the Kain-Fritsch (KF) cumulus scheme for the outer grids, the Morrison microphysics scheme (Morrison et al., 2009), the Noah Land Surface Model (Tewari et at., 2004), and the RRTMG longwave and shortwave radiation scheme (Iacono et al., 2008).

    Liang et al. (2019) put forward a pragmatic and effective approach in precipitation forecasting to avoid the challenge in representing convection across the gray zone. The convective gray zone refers to the model grid spacing around 1-10 km, where parameterized and resolved convective clouds could exist simultaneously. In the gray zone,many widely used cumulus parameterization assumptions become invalid and it is difficult to realistically represent the convection (Yano et al., 2010). Liang et al. (2019) found that double-nesting simulations using the WRF model with a large grid ratio (15:1 or 9:1) outperformed the traditional triple nesting with a middle (3- or 5-km) grid in forecasting the extreme rainfall of Jiangsu province. In particular, the outer grid using the KF scheme (Kain, 2004) to parameterize cumuli at 15 km with explicitly resolving convection at 1 km produced the best forecast of hourly rainfall variation.Therefore, we followed the same double-nesting experiments using the KF cumulus scheme in the outer coarse(15- and 9-km) grids and explicitly resolving convection at the inner 1-km grid. Our focus, however, is to examine the PBL effects on the prediction of extreme rainfall.

    Three PBL schemes were used and the related diffusion equations and eddy diffusivity coefficient computations are listed in Table 1. These schemes are all based on the K-gradient transport theory, which determines the turbulent flux by multiplying the eddy moment or heat diffusivity coefficient with the vertical gradient of grid-mean variables. Their main differences are the closure assumptions used to define the eddy diffusivity coefficient and the nonlocal effect for the energy exchange between model layers.The YSU scheme is a first-order closure scheme that determines the K-profile in the mixed layer by introducing the PBL height [Eq. (3)] but depends on the mixing length and Richardson number above the entrainment zone [Eqs. (1)and (2)]. A counter gradient transport term [Eq. (7)] related to surface buoyancy flux is included for the nonlocal effect,especially for unstable boundary layers (Hong et al., 2006).As higher-order schemes, both MYNN (which hereafter refers to MYNN 2.5 in WRFv3.9.1.1) and MYJ (Janji?,1994) are based on TKE closure to parameterize the eddy diffusivity but without considering the nonlocal effect in the unstable layer. They differ in the diagnostic equations for dimensionless stability function (S) and turbulent length scale (l). The MYJ (Janji?, 1994) scheme represents a non-singular implementation of the Mellor-Yamada level-2.5 turbulence closure model by adding limitation to the turbulent length scale. The MYNN scheme was considered to better represent the mixing in the convective boundary layer than MYJ because it considers the buoyancy effect on the stability function and turbulent length scale (Srinivas et al.,2018). In WRFv3.9.1.1, the MYNN scheme includes several options to improve the coupling of the PBL scheme with radiation (icloud_bl=1) and microphysics (bl_mynn_cloudmix=1), and two options (bl_mynn_edmf=1, bl_mynn_mixlength=2) to use the cloud-specific and scale-aware mixing length following Ito et al. (2015). As suggested by the WRF physics documentation (Wang and Bruyere, 2017),the YSU scheme adopts the MM5 (Jiménez et al., 2012) surface layer scheme, and the MYJ and MYNN2.5 schemes both use the Monin-Obukhov (Monin and Obukhov, 1954)surface layer scheme.

    Table 2 summarizes the configurations for two groups of double-nesting simulations. The mesoscale grids (15 and 9 km) used KF-parameterized convection, while the 1-km grid used fully explicit convection (EC). All experiments were initialized at 0800 LST 24 May 2018, and integrated for 48 h up to 0800 LST 26 May 2018. For convenience, D1 and D2 are denoted as the outer and inner domains, respectively. Nesting grid configurations are denoted directly by the grid spacing in km in sequential order. For example, 15-1 km denotes a double nesting configuration between D1 at 15 km and D2 at 1 km. The model results were bilinearly interpolated to the AWS stations to facilitate the comparison.

    The evaluation methods include pattern correlations(COR; Barnston, 1992), root-mean-square errors (RMSE;Barnston, 1992), threat score (TS; Wilks, 2011), and bias score (BS; Wilks, 2011) for different precipitation intensity thresholds (i.e., 0.1, 10, 25, 50, 100 mm, representing light,moderate, large, heavy, and extreme rain respectively). Larger COR and smaller RMSE indicate a better forecast of the spatial pattern of rainfall amount. The TS is also called the Critical Success Index, and a larger TS indicates higher predictive skill (perfect = 1) for the corresponding precipitation intensity. The bias score (Bias) denotes the frequency of rainfall forecasts compared with observations, and BS = 1 indicates an ideal prediction of the relative area size of corresponding precipitation intensity. The equations of COR,RMSE, TS, and BS are as follows:

    Table 1. Diffusion equations and vertical diffusivity coefficient computations for the three PBL parameterization schemes used in the model experiments.

    Table 2. Two groups of model experiments using different configurations of grid nesting and PBL parameterization schemes as well as convection treatments. 15-1 km and 9-1 km denote the experiments in 1-km grid nested with the 15 km- and 9 km-grid, respectively.

    4. Sensitivity of heavy rainfall simulations to PBL schemes

    This section examines the sensitivity of heavy rainfall simulations to three PBL schemes in the mesoscale grids(15 and 9 km) and the high-resolution grids (1 km). The model evaluation will be focused on the spatial distributionof 24-h accumulated precipitation and hourly variations of regional mean precipitation over the core rainband (24-h accumulated precipitation > 100 mm).

    Table 3. Four categories of the occurrence/non-occurrence for a rain event at each threshold. Four categories of hit (a), miss (b),false alarm (c) and correct non-rain forecast (d) are used to refer to the occurrence/non-occurrence of a rain event at each threshold.

    4.1. 24-h accumulated precipitation

    Fig. 4. Spatial distributions of 24-h accumulated rainfall on 25 May 2018 simulated by WRF in D1 (15 km, 9 km)and D2 1-km grids with three different PBL schemes (YSU, MYJ and MYNN) (units: mm).

    Figure 4 compares the 24-h accumulated precipitation from 0000 LST 25 May to 2300 LST 25 May using different grid spacing (15 km, 9 km, 1 km) and PBL schemes(YSU, MYJ, and MYNN). For the D1 15-km grid, all PBL schemes failed to capture the location of the observed rain band and the extreme rainfall amount was largely underestimated. As grid spacing was reduced to 9 km, these dry biases for all PBL schemes were reduced, but this still underestimated the west rainband and overestimated the coverage of the east rainband. The model differences between these three PBL schemes were not significant, implying that the extreme rainfall simulations at coarser grids are not sensitive to the choice of PBL scheme. In the inner D2 1-km grid simulations, the model biases in the outer coarser grids were greatly reduced and the simulation of rainfall structure was refined, but the differences in using the three PBL schemes became more significant. Specifically, using the YSU scheme in 1-km grids generally underestimated the east rainband nearby the Yangtze River Delta. The general pattern of the core rainband was better represented in the experiments using the MYJ and MYNN schemes, but the former produced a rainband that deviated too far north.

    Figure 5 further compares the TS and BS statistics for different daily rainfall intensity thresholds, as well as COR and RMSE for all 5207 stations in the inner 1-km grid from all the double-nesting simulations. In all the 1-km grid simulations, TS scores for light rain were generally higher than those for moderate-heavy rain. This is similar to the study of Liang et al. (2019), in which the WRF model using the KF cumulus parameterization scheme at 15- and 9-km grids produced higher TSs for light-moderate rain but lower scores for all other categories in a monthly prediction over Jiangsu Province. Thus, we speculate that this might be because the WRF model has limitations in representing physical (microphysics, PBL, etc.) and dynamic processes, leading to lower skill in predicting the convective cells than the general rainfall patterns. More importantly, the three PBL schemes showed notable differences in predicting the intensity of extreme rainfall and the size of its relative area. For instance, in the double nesting of 9-1 km simulations, the MYNN scheme was superior in predicting the intensity and relative area of extreme rainfall, with the highest TS and more realistic BS value; while the YSU scheme only performed better in predicting the intensity of light rain, and even performed worse in predicting the intensity and relative area of extreme rainfall, with the lowest TS and BS values. Using the MYJ scheme showed advantages in predicting the moderate rain, with the highest TS and a more realistic BS (close to 1), but had medium predictive skill for the intensity and area of extreme rainfall. Furthermore, double nesting of 9-1 km performed slightly better than 15-1 km in predicting the entire precipitation distribution, with systematically higher COR and lower RMSE.

    Therefore, among all the double-nesting simulations,using MYNN in the 1-km grid nested with the 9-km grid produced the best predictive skill for the entire precipitation distribution, intensity, and size of the core rainband. This can be demonstrated by the largest COR for all the stations, as well as TS and BS (0.68, 0.34, 0.93), for the extreme rainfall threshold compared to YSU (0.63. 0.12, 0.40) and MYJ(0.76, 0.24, 0.60).

    4.2. Hourly rainfall variations

    Figure 6 compares the hourly variations of observed and simulated rainfall averages over the two subregions(Fig. 1b) and the entire core rainband in the inner 1-km grid as well as the outer D1 coarser grids (15 and 9 km). In observations, the rainfall amount was gradually increasing from 0000 LST and reached an early-morning peak at around 0700 LST. The primary rainfall peak occurred at around 1200 LST and then weakened up until 1600 LST, after which the rainfall increased again followed by a late-afternoon secondary peak at 1800 LST due to the intensification of the nocturnal LLJ and the mesoscale convective systems(Fig. 2). In the outer 15-km or 9-km grid (Figs. 6a and b),all three simulations systematically underestimated the primary and secondary rainfall peaks. The YSU scheme showed slight advantages over the MYNN and MYJ schemes in capturing the peak time as well as the intensity,while the MYJ scheme had the worst performance by producing a 1-h earlier peak with the lowest intensity.

    Fig. 5. The (a) TS, (b) BS and (c) COR/RMSE for all 5207 stations with parameterized PBLs (YSU, MYJ and MYNN) for the inner D2 grids.

    Fig. 6. Hourly variations of observed (OBS) and simulated rainfall averages over the (a-d) core rainband and (e-h)two subregions [(e, f) East; (g, h) West] on 25 May 2018 by the WRF model using three different PBL schemes(units: mm h-1).

    In the inner 1-km grid nested with the outer 15-km or 9-km grid (Figs. 6c and d), consistent advantages existed for the MYNN scheme over the other schemes in capturing the primary rainfall peak and intensity. The YSU scheme was capable of capturing the primary rainfall peak time, but with lower intensity than that of the MYNN scheme. However,using the MYJ scheme would still produce a 1-h earlier rainfall peak with the lowest intensity in the double nesting of 15-1 km, which was improved in the double nesting of 9-1 km.

    By analyzing the hourly rainfall simulations over two subregions, we found that using the MYNN scheme in the double nesting of 9-1 km also reproduced the primary rainfall peaks well over both subregions (Figs. 6f and h).However, using the MYNN scheme in double nesting of 15-1 km would overestimate the peak intensity over the east rainband (Fig. 6e), also producing an earlier and weaker peak over the west rainband (Fig. 6g) compared to the observations. Furthermore, the largely underestimated rainfall peaks simulated by the MYJ scheme for both 15-1 km and 9-1 km were all mainly due to the dry biases over the west rainband. Therefore, we will focus on the double nesting of 9-1 km using the MYNN scheme, and explore the possible causes for its superiority in simulating the primary rainfall peak intensity. The analyses mainly include the PBL effects on the boundary layer structure, convective thermodynamic,and large-scale forcing variations.

    5. Understanding the PBL sensitivities

    5.1. Simulations of boundary layer structure

    Fig. 7. Time-height cross sections of the vertical eddy transport of (a-c) equivalent potential temperature (unit:K s-1) and (d-f) water vapor [units: kg (cm hPa s)-1] averaged over the core rainband on 25 May in the inner 1-km grids (9-1) simulations using YSU, MYJ and MYNN.

    Figure 7 compares the time-height sections of vertical eddy transport of equivalent potential temperature and water vapor mixing ratio. They are averaged over the core rainband from the double nesting of 9-1 km using the three PBL schemes. The simulated vertical transport of heat and water vapor in the boundary layer all gradually became stronger from the early morning to noon, corresponding to the time when the primary rainfall peak occurred. The major difference existed with the vertical transport and mixing process at the top of the boundary layer between the nonlocal(YSU) and local (MYNN & MYJ) PBL schemes. The YSU scheme includes the nonlocal effect with the counter gradient transport and produced stronger vertical mixing below 2000 m than that of the local MYNN scheme before 1200 LST. The MYJ scheme showed intermediate vertical mixing between the YSU and MYNN schemes, as supported by Srinivas et al. (2018). This would cause differences in the vertical distribution of heat and water vapor from the boundary layer to the lower troposphere, and thus affect the stability of the atmosphere as well as the development of convection and precipitation, which will be discussed as follows. It is noteworthy that there was high vertical equivalent potential temperature and moisture transport at 2000 LST in the experiment with the MYNN scheme, which may explain the high precipitation at 2100 LST shown in Fig. 6.

    Figure 8 compares the hourly variation of vertical distributions of equivalent potential temperature, water vapor mixing ratio, horizontal wind speed, and vertical velocity averaged over the core rainband from the surface to the lower troposphere. The vertical profiles of these corresponding variables at 1200 LST in Fig. 9 are also combined to discuss the different effects of the PBL schemes. For the thermal structures, the three PBL schemes all produced an inversion layer near the surface and the depth was decreasing from the early morning till noon, but the local MYNN scheme produced the shallowest inversion layer at 1200 LST around 1000 m compared to the nonlocal YSU scheme (1200 m) or the local MYJ scheme (1300 m) (Fig. 9a). Meanwhile, they all produced a warm layer in the middle-upper boundary layer, resulting in an unstable layer above. For the water vapor distribution, the three PBL schemes showed relatively little difference, but the water vapor was still more evenly distributed within the inversion layer (< 1000 m) during morning till noon in the nonlocal YSU scheme. The local MYJ and MYNN schemes both showed a more concentrated water vapor center in the bottom layer, but the MYJ scheme produced slightly higher water vapor mixing ratios from the top of the boundary layer to the middle troposphere at 1200 LST (Fig. 9b). For the wind structure, the local MYNN scheme produced stronger horizontal wind speed and upward motion than the nonlocal YSU scheme from the lower troposphere to the upper levels. The wind fields simulated by the MYJ scheme were strong at the lower troposphere; however, they rapidly reduced from the lower troposphere to upper levels (Figs. 9c and d).However, the nonlocal YSU scheme produced the lowest wind speed from the low-level troposphere (2500 m;13.60 m s) to upper levels (6000 m; 17.42 m s) owing to its strongest vertical mixing at the top of the unstable boundary layer (Figs. 7a and d).

    Fig. 8. Time-height sections of the (a-c) equivalent potential temperature (units: K), (d-f) water vapor mixing ratio (units:g kg-1), (g-i) horizontal wind speed (unit: m s-1), and (j-l) vertical velocity (units: m s-1) averaged over the core rainband 1-km grids (9-1) using YSU, MYJ and MYNN.

    Therefore, the local MYNN scheme produced the shallowest and most humid inversion layer in the bottom layer,the least warm middle boundary layer, but stronger horizontal wind and upward motion from the top of the boundary layer to the upper levels. YSU produced a deeper inversion layer than MYNN in the bottom layer, a relatively warm and humid middle boundary layer, and the weakest wind fields due to the strongest vertical mixing at the top of the boundary layer. The local MYJ scheme produced an inversion layer with a depth between the other two schemes, but the warmest middle boundary layer, and the wind fields of MYJ at the boundary layer top were strong but quickly weakened above 5000 m. The MYJ scheme had slightly higher wind speeds in the boundary layer compared to the others, but caution should be taken here considering the uncertainty of wind observations and computational errors. In the following we focus on the effects of these different boundary layer structures on the development of convection.

    5.2. Simulations of convective thermodynamic variations

    Fig. 9. Profiles of (a) equivalent potential temperature (θe, units: K), (b) water vapor mixing ratio (QV, units:g kg-1), (c) horizontal wind speed (wspd, units: m s-1), and (d) vertical velocity (W, units: m s-1) from 100 to 6000 m over the core rainband at 1200 LST in the 1-km (9-1) grids using YSU, MYJ and MYNN.

    Figure 10 compares the hourly variations of simulated convective available potential energy (CAPE) and convective inhibition energy (CIN) averaged over the core rainband from the double nesting of 9-1 km simulations with the three PBL schemes. CAPE is defined as the accumulated buoyancy energy from the level of free convection to the equilibrium level, which represents the net kinetic energy that a rising air parcel can gain from the environment and hence reflects the potential strength of convective systems and associated precipitation. CIN is defined as the accumulated negative buoyant energy from the air parcel’s starting position to the level of free convection, which represents the inhibition energy that the rising air parcel has to overcome to become free convection. Under stormy conditions,the PBL height could be defined by some alternative methods such as using a cloud base or lifting condensation level(LCL) (Wisse and de Arellano, 2004; Stull, 2011), and a higher LCL usually denotes a stronger vertical mixing.Table 4 also lists the simulated values of CAPE, CIN and LCL, as well as the horizontal wind speed and upward motions at the top of the boundary layer (around 1500 m),the low-level troposphere (around 2500 m), and the middle troposphere (around 6000 m) at 1200 LST for a quantitative comparison.

    The nonlocal YSU scheme produced the lowest CAPE peak at around 0800 LST, which was two hours earlier than the primary rainfall peak. Although the YSU scheme produced comparable CIN with the MYNN scheme, it systematically had the highest LCL compared to the other two schemes. This implies that, compared to the local MYNN scheme, the air parcel in the YSU scheme needs to be uplifted to a higher altitude to achieve condensation; nevertheless, the unstable energy in the upper boundary layer is insufficient to sustain the development of convection.

    Fig. 10. Hourly variations of simulated (a) CAPE (units: J kg-1), (b) CIN(units: J kg-1) and (c) LCL (units: m) averaged over the core rainband at 1-km(9-1) grids on 25 May using three PBL schemes.

    Table 4. Simulated values of CAPE, CIN and LCL, as well as the horizontal wind speed, vertical velocity and water vapor mixing ratio,at the top of boundary layer (around 1500 m), the low-level troposphere (around 2500 m), and the middle troposphere (around 6000 m) at 1200 LST.

    However, the local MYJ scheme systematically produced the largest CAPE and CIN during the morning through the afternoon, indicating that the convection was strongly suppressed and the CAPE was slowly released. The LCL in the MYJ scheme was also higher than in the MYNN scheme. The larger CIN, the slower release rate of CAPE,and the higher LCL with the MYJ scheme all determine that it is more difficult for the convection to initiate and develop compared to using the MYNN scheme. The atmospheric instability, low-level moisture convergence and vertical motion are the prerequisites for the development and maintenance of deep convection and mesoscale convective systems that often lead to heavy rainfall events (Srinivas et al.,2018). Therefore, in the following, we discuss the PBL effects on the large-scale forcings, including the low-level moisture supply and the upward motions.

    5.3. Simulations of large-scale forcing variations

    Fig. 11. Time-height sections of (a-c) horizontal water vapor mixing ratio flux [units: kg (cm hPa s)-1] and (d-f)vertical water vapor mixing ratio flux [units: kg (cm hPa s)-1] and vertical velocity (black solid contours; units: m s-1)averaged over the core rainband on 25 May in the inner 1-km (9-1) grids using YSU, MYJ and MYNN.

    Figure 11 compares the time-height sections of horizontal and vertical fluxes of water vapor mixing ratio and vertical velocity averaged over the core rainband from the surface to the troposphere in the double nesting of 9-1 km simulations with the three PBL schemes. The water vapor transport simulated by these three PBL schemes mainly differed in its contributions to the primary rainfall peak at noon. The local MYNN scheme systematically produced the strongest horizontal water vapor transport, as well as upward motion from the top of the boundary layer to the middle troposphere during 1000-1200 LST. This can also be identified from Table 4, in which the MYNN scheme shows slightly stronger wind speeds and vertical velocity at the boundary layer top and the middle troposphere at 1200 LST compared to YSU. Although the MYJ scheme produced higher wind speeds and upward motion at 1200 LST at the top of the boundary layer, they weakened rapidly and produced less moisture transport to the upper levels. Therefore, the MYJ scheme also predicted a lower intensity of the primary rainfall peak compared to the MYNN scheme.

    Figure 12 presents longitude-height sections of moisture convergence and zonal-vertical wind averaged over 31°-33°N in the inner 1-km grid from the double nesting of 9-1 km simulations at 1200 LST and 1400 LST using the three PBL schemes. Here, we chose to demonstrate the moisture convergence at 1400 LST because of the availability of ECMWF reanalysis data. In the ECMWF data, at 1400 LST, there was strong low-level moisture convergence over the region 120°-122°E associated with strong upward motion from the lower to upper troposphere, corresponding to the east rainband. MYNN produced stronger upward motion than the other PBL schemes, which affected the amount of vertical moisture transport and produced the strongest low-level moisture convergence over this region.Weak low-level moisture convergence and upward motion in the lower troposphere were also shown in the ECMWF data over the region 117°-118°E, corresponding to the relatively weak rainfall in the west rainband. However, using the MYNN scheme could only partly capture the upward motion and moisture convergence over this region, and with underestimated intensity; and yet, the other two schemes could barely capture the upward motion, and the MYJ scheme even produced weak moisture divergence in the lower troposphere. This also explains why the experiments with the YSU and MYJ schemes largely underestimated the extreme rainfall over the west rainband. Besides, at 1200 LST, MYNN also produced stronger low-level moisture convergence and upward motion in the lower troposphere over the regions 120°-122°E and 117°-118°E.

    Fig. 12. Longitude-height sections of moisture convergence [shaded; units: 10-7 g (cm2 hPa s)-1] and U-W wind (black vectors; units: U × 100 W m s-1) averaged over 31°-33°N on 25 May in the inner 1-km nested with 9-km grid using three PBL schemes at 1200 LST (a-c) and 1400 LST (d-f).

    Therefore, the nonlocal YSU scheme, with the strongest vertical mixing, produced the deepest inversion layer in the bottom layer with the highest LCL, the lowest CAPE, and the weakest wind fields, as well as low-level moisture transport. The MYJ scheme, with intermediate vertical mixing, produced larger CIN, slower release of CAPE, and higher LCL compared to the MYNN scheme, which suppressed the development of convection. The wind fields in the MYJ scheme at the low-level troposphere were the strongest but quickly weakened upwards, resulting in less moisture transport and convergence in the lower troposphere, especially over the west rainband. However, the MYNN scheme, with the weakest vertical mixing, produced the shallowest and most humid inversion layer in the bottom layer with the lowest LCL, but stronger wind fields and upward motions from the boundary layer top to upper levels. These all facilitated the development of deep convection and moisture transport for intense precipitation.

    6. Conclusions

    High-resolution predictions of an extreme heavy rainfall event were carried out in this study using the WRF model based on double nesting with two different grid ratios(15:1 and 9:1). The sensitivity of the precipitation forecast was examined using three PBL parameterization schemes,including the nonlocal first-order YSU scheme and the local higher-order MYJ and MYNN schemes. The extreme precipitation event occurred in Shanghai and its nearby area on 25 May 2018 and it involved the interaction of multi-scale systems from the synoptic circulation, mesoscale LLJ, and low-level shear line to local convective activities. Various statistical measures were adopted to quantitatively evaluate the PBL effects on daily rainfall distributions and the hourly variations of extreme rainfall. The impacts of the PBL schemes on the simulations of the boundary layer structure, convective thermodynamic, and large-scale forcing variations were further analyzed. The main conclusions are as follows:

    First, in the outer mesoscale grids (15 km and 9 km),all experiments failed to capture the location of the observed rainband and the extreme rainfall amount was largely underestimated. The model differences between the three PBL schemes were not significant. In the inner 1-km grid simulations, the model biases in the outer coarser grids were greatly reduced and rainfall-structure simulation was refined. The east rainband was generally underestimated with the YSU scheme, while using both the MYJ and MYNN scheme could capture the general pattern of the core rainband, but the former produced a rainband that deviated further north. Among all the double nesting simulations,using MYNN in the 1-km grid with explicitly resolving convection nested with the 9-km grid with the KF cumulus scheme, produced the best predictive skill for the entire precipitation distribution and the intensity, as well as the size of the core rainband.

    Second, in the outer mesoscale grids (15 km and 9 km),all three simulations systematically underestimated the primary and secondary rainfall peaks. In the inner 1-km grid nested with the outer 15-km or 9-km grid, the MYNN scheme showed consistent advantages over the other two PBL schemes in capturing the primary rainfall peak and intensity. Using the MYNN scheme in the double nesting of 15-1 km, however, the peak intensity over the east rain band was overestimated, and an earlier as well as weaker peak over the west rainband was produced. Nevertheless, the MYNN scheme in the double nesting of 9-1 km showed superiority over the other 1-km simulations in predicting the primary rainfall peaks over both the west and east core rainbands.

    Third, the three PBL schemes differ in their vertical mixing processes, leading to different vertical profiles of temperature, water vapor and wind fields from the boundary layer to the lower troposphere, which affect the atmospheric stability and moisture transport for the development of convection and precipitation. Among the PBL schemes, for this extreme rainfall event, the local MYNN scheme showed the weakest vertical mixing and produced the shallowest as well as the most humid inversion layer in the bottom layer and the lowest LCL, which are conducive to the development of convection. At the same time, the MYNN scheme produced stronger horizontal winds and upward motion from the top of the boundary layer to the upper levels than the nonlocal YSU scheme, which facilitates the moisture transport and convergence in the lower troposphere for intense precipitation. The combined effects led to the strongest and most realistic rainfall peak in the MYNN scheme. However, the nonlocal YSU scheme, with the strongest vertical mixing, produced a deep inversion layer,the highest LCL and the lowest CAPE, and weaker wind fields, which tend to inhibit the development of convection and moisture transport for heavy rainfall. This explains the largely underestimated primary rainfall peak in the YSU scheme. On the other hand, the MYJ scheme, with intermediate vertical mixing, produced larger CIN and a higher LCL with a slower release of CAPE compared to the MYNN scheme, which also suppressed the development of convection. Although the horizontal winds of the MYJ scheme were comparably strong with the MYNN scheme in the upper boundary layer, it weakened rapidly in the lower troposphere, resulting in less moisture transport than the MYNN scheme. This especially led to the underestimated rainfall amount over the west rainband.

    Although our study was based on one single-day case over Shanghai and its adjacent region, the model performances of these three PBL schemes in this case were highly similar to those of previous high-resolution model studies,such as Srinivas et al. (2018), but over different regions. On the other hand, the PBL sensitivities with the varying model resolution and convection treatment also provide a better understanding of the interaction of the boundary layer with the large-scale environment and the organization of mesoscale convective systems. In future work we will conduct high-resolution simulations over a longer period and under different climate regimes to verify the typical characteristics of different PBL schemes. Moreover, the degradation of 1-km simulations downscaled from the coarser grids (15 km and 9 km) in predicting the weak secondary rainfall peak requires further study through improving the cumulus or microphysics representations.

    This research was supported by the National Natural Science Foundation of China (Grant No.41730646), National Natural Science Foundation for Young Scientists of China (Grant No. 41605079), and the National Key R&D Program of China (Grant No. 2018YFC1507702). The FNL Operational Global Analysis data at 1° grids provided by NCEP and used in this study can be found at https://rda.ucar.edu/datasets/ds083.2/index.html. The ERA-Interim Global Atmospheric Reanalysis Data (0.125° × 0.125° grids) provided by the ECMWF and used to analyze the evolution of horizontal wind speed can be found at https://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/. The model simulations were conducted at the ECNU Multifunctional Platform for Innovation 001 facilities, and data were deposited in the ECNU public data server (IP 49.52.29.112). The views expressed are those of the authors and do not necessarily reflect those of the sponsoring agencies.

    老司机福利观看| 水蜜桃什么品种好| 日本欧美视频一区| 精品国产乱码久久久久久男人| 十八禁人妻一区二区| 免费久久久久久久精品成人欧美视频| 老熟妇乱子伦视频在线观看| avwww免费| 欧美日韩乱码在线| 老司机午夜福利在线观看视频| 日日夜夜操网爽| 999久久久精品免费观看国产| 精品国产一区二区三区四区第35| 啪啪无遮挡十八禁网站| 首页视频小说图片口味搜索| 免费日韩欧美在线观看| 亚洲成人精品中文字幕电影 | 中文字幕最新亚洲高清| 亚洲精品一卡2卡三卡4卡5卡| 免费少妇av软件| 90打野战视频偷拍视频| 国产区一区二久久| 国产成+人综合+亚洲专区| 乱人伦中国视频| 琪琪午夜伦伦电影理论片6080| 制服诱惑二区| 中文字幕高清在线视频| 国产男靠女视频免费网站| 91成年电影在线观看| 久99久视频精品免费| 91老司机精品| 精品电影一区二区在线| 琪琪午夜伦伦电影理论片6080| 99香蕉大伊视频| 黄色怎么调成土黄色| 国产亚洲精品久久久久久毛片| 国产高清激情床上av| 亚洲av熟女| 老司机午夜福利在线观看视频| 中文字幕精品免费在线观看视频| 日韩精品免费视频一区二区三区| 交换朋友夫妻互换小说| 欧美另类亚洲清纯唯美| 亚洲成人免费电影在线观看| 长腿黑丝高跟| 精品无人区乱码1区二区| 国产乱人伦免费视频| 日韩免费av在线播放| 在线观看一区二区三区激情| 日本 av在线| 国产一区二区激情短视频| 欧美人与性动交α欧美精品济南到| 国内久久婷婷六月综合欲色啪| 老熟妇乱子伦视频在线观看| 亚洲成人国产一区在线观看| 久久午夜亚洲精品久久| 国产1区2区3区精品| 女人被躁到高潮嗷嗷叫费观| www.999成人在线观看| 久久精品91无色码中文字幕| 多毛熟女@视频| 成年版毛片免费区| 国产精品免费一区二区三区在线| 又大又爽又粗| 日韩大尺度精品在线看网址 | 精品乱码久久久久久99久播| 欧美成狂野欧美在线观看| 国产aⅴ精品一区二区三区波| 又紧又爽又黄一区二区| 亚洲成av片中文字幕在线观看| 成人手机av| 91老司机精品| 色哟哟哟哟哟哟| 天堂√8在线中文| 99在线视频只有这里精品首页| 高潮久久久久久久久久久不卡| 欧美亚洲日本最大视频资源| 在线天堂中文资源库| 久9热在线精品视频| 淫妇啪啪啪对白视频| 亚洲一区二区三区色噜噜 | 伦理电影免费视频| 12—13女人毛片做爰片一| 深夜精品福利| 日本vs欧美在线观看视频| 韩国av一区二区三区四区| 自线自在国产av| 日本黄色日本黄色录像| 无人区码免费观看不卡| 不卡av一区二区三区| 国产又爽黄色视频| 一边摸一边抽搐一进一小说| 精品一区二区三区av网在线观看| 黄色 视频免费看| 亚洲精品国产精品久久久不卡| 视频区图区小说| 搡老熟女国产l中国老女人| 欧美日本亚洲视频在线播放| 男女下面插进去视频免费观看| 国产激情久久老熟女| 侵犯人妻中文字幕一二三四区| 在线免费观看的www视频| 真人一进一出gif抽搐免费| 欧美最黄视频在线播放免费 | 国产av又大| 亚洲 欧美一区二区三区| 自线自在国产av| 亚洲九九香蕉| 男女高潮啪啪啪动态图| 一进一出抽搐gif免费好疼 | 亚洲欧美激情在线| 国产精品日韩av在线免费观看 | 两性午夜刺激爽爽歪歪视频在线观看 | 亚洲精品粉嫩美女一区| 大码成人一级视频| 久久久久久免费高清国产稀缺| 国产成年人精品一区二区 | a级毛片黄视频| 村上凉子中文字幕在线| 精品乱码久久久久久99久播| а√天堂www在线а√下载| 欧美另类亚洲清纯唯美| 12—13女人毛片做爰片一| 精品久久久久久久毛片微露脸| 久久久国产成人精品二区 | 亚洲五月婷婷丁香| 人妻久久中文字幕网| 波多野结衣一区麻豆| 在线看a的网站| 中文字幕人妻熟女乱码| 高清毛片免费观看视频网站 | 两人在一起打扑克的视频| 色尼玛亚洲综合影院| 国产成+人综合+亚洲专区| ponron亚洲| 亚洲第一青青草原| 老鸭窝网址在线观看| 大码成人一级视频| 99久久99久久久精品蜜桃| 国产高清视频在线播放一区| 久久精品国产清高在天天线| 在线看a的网站| 国产区一区二久久| 日日摸夜夜添夜夜添小说| 激情视频va一区二区三区| 国产午夜精品久久久久久| 国产成人一区二区三区免费视频网站| 欧美日韩精品网址| 成人亚洲精品一区在线观看| 伦理电影免费视频| 亚洲伊人色综图| 天堂中文最新版在线下载| www.熟女人妻精品国产| 97碰自拍视频| 男女下面插进去视频免费观看| 久久国产精品人妻蜜桃| 亚洲 国产 在线| 狠狠狠狠99中文字幕| 国产高清国产精品国产三级| 久久久久久人人人人人| 亚洲aⅴ乱码一区二区在线播放 | av网站免费在线观看视频| 久久人妻av系列| 日韩欧美国产一区二区入口| 啪啪无遮挡十八禁网站| 国产成人免费无遮挡视频| 国产在线观看jvid| www.www免费av| 麻豆一二三区av精品| 99re在线观看精品视频| 成人黄色视频免费在线看| 国产1区2区3区精品| 精品国产国语对白av| 中文字幕精品免费在线观看视频| 99久久综合精品五月天人人| 日日爽夜夜爽网站| 动漫黄色视频在线观看| 精品一品国产午夜福利视频| 久久人妻福利社区极品人妻图片| 后天国语完整版免费观看| 天堂影院成人在线观看| 久久久久亚洲av毛片大全| 免费人成视频x8x8入口观看| 窝窝影院91人妻| 亚洲第一av免费看| 午夜免费激情av| 久久久久久久久中文| avwww免费| 69精品国产乱码久久久| 欧美精品啪啪一区二区三区| 女生性感内裤真人,穿戴方法视频| 首页视频小说图片口味搜索| 午夜两性在线视频| 中文字幕人妻丝袜制服| 99精品在免费线老司机午夜| 亚洲片人在线观看| 美女 人体艺术 gogo| 欧美不卡视频在线免费观看 | 日本三级黄在线观看| 12—13女人毛片做爰片一| 天堂中文最新版在线下载| 人人澡人人妻人| 在线视频色国产色| 亚洲人成电影免费在线| 岛国在线观看网站| 日韩视频一区二区在线观看| 人妻丰满熟妇av一区二区三区| 一级a爱片免费观看的视频| 成人手机av| 老司机深夜福利视频在线观看| 人人妻人人爽人人添夜夜欢视频| 亚洲男人的天堂狠狠| 午夜视频精品福利| cao死你这个sao货| 啦啦啦在线免费观看视频4| 日韩精品免费视频一区二区三区| 国产精品秋霞免费鲁丝片| 97人妻天天添夜夜摸| 国产午夜精品久久久久久| 精品久久久久久成人av| 日本五十路高清| 一区二区三区精品91| 国产成+人综合+亚洲专区| 亚洲成av片中文字幕在线观看| 中出人妻视频一区二区| 日韩精品青青久久久久久| 久久精品91无色码中文字幕| 不卡av一区二区三区| 99在线人妻在线中文字幕| 精品日产1卡2卡| 18禁观看日本| 男人操女人黄网站| 色婷婷久久久亚洲欧美| 久久久精品欧美日韩精品| 久久人妻av系列| 国产在线观看jvid| 免费不卡黄色视频| 亚洲人成网站在线播放欧美日韩| 久久人妻福利社区极品人妻图片| 男男h啪啪无遮挡| 国产av一区二区精品久久| 精品国产一区二区久久| 黄网站色视频无遮挡免费观看| 黄频高清免费视频| 黑人操中国人逼视频| 日韩欧美三级三区| 高清在线国产一区| 啦啦啦免费观看视频1| 久99久视频精品免费| 变态另类成人亚洲欧美熟女 | 少妇的丰满在线观看| 黄色 视频免费看| 亚洲国产欧美一区二区综合| 欧美 亚洲 国产 日韩一| 91精品三级在线观看| 熟女少妇亚洲综合色aaa.| 999久久久国产精品视频| 久久精品成人免费网站| 亚洲av五月六月丁香网| а√天堂www在线а√下载| 真人做人爱边吃奶动态| 91精品三级在线观看| 法律面前人人平等表现在哪些方面| www.自偷自拍.com| 久久性视频一级片| 亚洲五月色婷婷综合| 一夜夜www| 成人精品一区二区免费| 天堂√8在线中文| 男女高潮啪啪啪动态图| 久久精品aⅴ一区二区三区四区| x7x7x7水蜜桃| 国产欧美日韩一区二区三区在线| 久久九九热精品免费| 亚洲色图av天堂| 99精国产麻豆久久婷婷| 午夜日韩欧美国产| 亚洲熟妇熟女久久| 我的亚洲天堂| 少妇被粗大的猛进出69影院| 高清av免费在线| 国产三级在线视频| 日本一区二区免费在线视频| 午夜视频精品福利| 久久久国产成人免费| 在线观看免费视频日本深夜| 国产成人精品无人区| 在线观看舔阴道视频| 啦啦啦 在线观看视频| 国产精品一区二区在线不卡| 国产真人三级小视频在线观看| 日韩一卡2卡3卡4卡2021年| 999久久久精品免费观看国产| 啪啪无遮挡十八禁网站| 激情在线观看视频在线高清| 日日干狠狠操夜夜爽| 国产亚洲精品综合一区在线观看 | 国产精品乱码一区二三区的特点 | 最新在线观看一区二区三区| 亚洲精品在线观看二区| 老司机午夜十八禁免费视频| 乱人伦中国视频| 久久香蕉激情| 变态另类成人亚洲欧美熟女 | 黑人猛操日本美女一级片| 欧美日本亚洲视频在线播放| 桃色一区二区三区在线观看| 午夜影院日韩av| 日本wwww免费看| 老汉色∧v一级毛片| 国产黄色免费在线视频| 国产亚洲av高清不卡| av网站免费在线观看视频| 免费在线观看视频国产中文字幕亚洲| 成在线人永久免费视频| 色综合婷婷激情| 亚洲精品美女久久久久99蜜臀| 国产高清国产精品国产三级| 久久中文看片网| 男女高潮啪啪啪动态图| 免费av毛片视频| av网站在线播放免费| 波多野结衣av一区二区av| 人人妻人人澡人人看| 99久久久亚洲精品蜜臀av| 黄片大片在线免费观看| 日本黄色视频三级网站网址| 91国产中文字幕| 国产精品免费视频内射| 久久亚洲真实| 水蜜桃什么品种好| 18美女黄网站色大片免费观看| 亚洲色图av天堂| av天堂在线播放| av福利片在线| 嫩草影视91久久| 99在线视频只有这里精品首页| 天堂俺去俺来也www色官网| 中国美女看黄片| 国产精品99久久99久久久不卡| 久久99一区二区三区| 成人永久免费在线观看视频| 天堂中文最新版在线下载| 久久精品亚洲精品国产色婷小说| 黄频高清免费视频| 久久久久国产精品人妻aⅴ院| 午夜亚洲福利在线播放| 757午夜福利合集在线观看| 性色av乱码一区二区三区2| 国产熟女xx| 精品久久久久久久毛片微露脸| 在线观看免费高清a一片| 亚洲精品美女久久av网站| 久久久久久久久免费视频了| 天堂√8在线中文| 亚洲欧美日韩高清在线视频| 久久久水蜜桃国产精品网| 免费观看精品视频网站| 琪琪午夜伦伦电影理论片6080| 国产精品一区二区在线不卡| 欧美午夜高清在线| 国产精品免费一区二区三区在线| 日本三级黄在线观看| 涩涩av久久男人的天堂| a在线观看视频网站| 热99国产精品久久久久久7| 国产真人三级小视频在线观看| 精品国产超薄肉色丝袜足j| 91成年电影在线观看| 伦理电影免费视频| 亚洲av成人不卡在线观看播放网| 大香蕉久久成人网| 亚洲中文日韩欧美视频| 欧美日韩av久久| 亚洲精品美女久久久久99蜜臀| 级片在线观看| 啦啦啦免费观看视频1| 老司机午夜十八禁免费视频| 亚洲精品成人av观看孕妇| 欧美老熟妇乱子伦牲交| 无人区码免费观看不卡| 欧美日韩瑟瑟在线播放| 亚洲精品粉嫩美女一区| 国产成人欧美在线观看| 色老头精品视频在线观看| 亚洲成人国产一区在线观看| 色尼玛亚洲综合影院| 日韩 欧美 亚洲 中文字幕| 这个男人来自地球电影免费观看| 国产成+人综合+亚洲专区| 人人妻人人爽人人添夜夜欢视频| 欧美另类亚洲清纯唯美| 婷婷精品国产亚洲av在线| 日韩有码中文字幕| 久热爱精品视频在线9| 99国产综合亚洲精品| 午夜免费鲁丝| 亚洲第一青青草原| 成人影院久久| 国产不卡一卡二| 操美女的视频在线观看| 一边摸一边做爽爽视频免费| 69av精品久久久久久| 精品久久久精品久久久| 亚洲av美国av| 日韩欧美在线二视频| 欧美一区二区精品小视频在线| av超薄肉色丝袜交足视频| 在线观看免费高清a一片| 午夜福利,免费看| 夫妻午夜视频| av天堂在线播放| 男人的好看免费观看在线视频 | 精品福利观看| 欧美激情高清一区二区三区| 欧美最黄视频在线播放免费 | 首页视频小说图片口味搜索| 国产av在哪里看| 欧美日韩一级在线毛片| 色哟哟哟哟哟哟| 亚洲aⅴ乱码一区二区在线播放 | 91成人精品电影| 日韩大码丰满熟妇| 亚洲中文日韩欧美视频| 1024香蕉在线观看| 中文字幕色久视频| 亚洲 国产 在线| 久久九九热精品免费| 777久久人妻少妇嫩草av网站| 久久中文字幕一级| 99精品在免费线老司机午夜| 另类亚洲欧美激情| 亚洲成国产人片在线观看| 999久久久国产精品视频| 超色免费av| 精品福利观看| 天堂中文最新版在线下载| 男男h啪啪无遮挡| 国产精品一区二区三区四区久久 | 琪琪午夜伦伦电影理论片6080| 欧美日韩国产mv在线观看视频| 51午夜福利影视在线观看| 欧美精品亚洲一区二区| 黄片大片在线免费观看| 侵犯人妻中文字幕一二三四区| 一区二区三区国产精品乱码| 亚洲国产精品999在线| 免费久久久久久久精品成人欧美视频| 亚洲第一青青草原| 一区二区三区国产精品乱码| 在线看a的网站| 一边摸一边做爽爽视频免费| 一进一出抽搐动态| 在线观看一区二区三区激情| 国产片内射在线| 日本黄色日本黄色录像| 亚洲一区二区三区不卡视频| 1024视频免费在线观看| 国产成人欧美| 色婷婷av一区二区三区视频| 一a级毛片在线观看| 制服人妻中文乱码| 亚洲精品国产精品久久久不卡| 国产亚洲精品一区二区www| 人成视频在线观看免费观看| 天堂俺去俺来也www色官网| 国产极品粉嫩免费观看在线| 在线视频色国产色| 国内久久婷婷六月综合欲色啪| 亚洲熟妇中文字幕五十中出 | 久久久久久久久久久久大奶| 精品国产超薄肉色丝袜足j| 无人区码免费观看不卡| 国产又爽黄色视频| 欧美大码av| 国产三级黄色录像| 欧美 亚洲 国产 日韩一| 视频在线观看一区二区三区| 亚洲精品成人av观看孕妇| 国产欧美日韩一区二区精品| 国产av在哪里看| 欧美中文综合在线视频| 国产精品一区二区三区四区久久 | 亚洲一区高清亚洲精品| 亚洲国产精品999在线| 级片在线观看| av有码第一页| 国产精品乱码一区二三区的特点 | 午夜免费鲁丝| 99国产综合亚洲精品| 一二三四在线观看免费中文在| 99香蕉大伊视频| 欧美av亚洲av综合av国产av| 中文字幕最新亚洲高清| 久久久久久久午夜电影 | 两性夫妻黄色片| 97碰自拍视频| 亚洲精品国产一区二区精华液| 亚洲久久久国产精品| 欧美日韩国产mv在线观看视频| 母亲3免费完整高清在线观看| 亚洲精品国产区一区二| 日本撒尿小便嘘嘘汇集6| 在线观看一区二区三区| 人人妻人人添人人爽欧美一区卜| 亚洲九九香蕉| 黑人巨大精品欧美一区二区mp4| 人人妻人人澡人人看| av在线播放免费不卡| 黑人欧美特级aaaaaa片| av在线天堂中文字幕 | 国产伦人伦偷精品视频| 亚洲三区欧美一区| 免费不卡黄色视频| 久久天堂一区二区三区四区| 国产成人精品久久二区二区免费| 免费搜索国产男女视频| 午夜视频精品福利| 精品日产1卡2卡| 久久精品影院6| 一本大道久久a久久精品| 桃红色精品国产亚洲av| 新久久久久国产一级毛片| 正在播放国产对白刺激| 成人亚洲精品一区在线观看| 午夜激情av网站| 日韩有码中文字幕| 在线观看一区二区三区| 可以在线观看毛片的网站| 国产成人精品久久二区二区91| 麻豆久久精品国产亚洲av | 美女高潮喷水抽搐中文字幕| 一本综合久久免费| 男人的好看免费观看在线视频 | 国产日韩一区二区三区精品不卡| 嫩草影院精品99| 一级黄色大片毛片| 精品无人区乱码1区二区| 少妇被粗大的猛进出69影院| 在线国产一区二区在线| 国产精品 国内视频| 成人国产一区最新在线观看| 日本三级黄在线观看| 90打野战视频偷拍视频| 又紧又爽又黄一区二区| 日韩一卡2卡3卡4卡2021年| 久久天堂一区二区三区四区| 亚洲欧美一区二区三区久久| 嫁个100分男人电影在线观看| 免费看a级黄色片| 久久热在线av| 国产伦人伦偷精品视频| 国产成人欧美| 熟女少妇亚洲综合色aaa.| 国产精品九九99| 国产精品98久久久久久宅男小说| 精品国产一区二区久久| 亚洲视频免费观看视频| netflix在线观看网站| 嫁个100分男人电影在线观看| 亚洲成人精品中文字幕电影 | 精品国产美女av久久久久小说| 99国产精品99久久久久| 免费看十八禁软件| 91成人精品电影| 午夜免费鲁丝| 日韩av在线大香蕉| 日韩精品青青久久久久久| 成人影院久久| 免费女性裸体啪啪无遮挡网站| 亚洲一区二区三区色噜噜 | 国产亚洲精品第一综合不卡| 久久人人爽av亚洲精品天堂| 久久人妻熟女aⅴ| 免费在线观看影片大全网站| 亚洲精品国产一区二区精华液| 久久精品人人爽人人爽视色| 一边摸一边做爽爽视频免费| 国产单亲对白刺激| 啦啦啦 在线观看视频| 国产精品1区2区在线观看.| 极品人妻少妇av视频| 岛国视频午夜一区免费看| 欧美精品啪啪一区二区三区| 日本一区二区免费在线视频| 午夜免费观看网址| 午夜福利影视在线免费观看| 国产欧美日韩一区二区精品| 一区二区日韩欧美中文字幕| 91在线观看av| 午夜精品在线福利| 一区福利在线观看| 欧美日韩国产mv在线观看视频| 制服诱惑二区| 一区二区日韩欧美中文字幕| 好看av亚洲va欧美ⅴa在| 成人手机av| 亚洲片人在线观看| 搡老熟女国产l中国老女人| 欧美色视频一区免费| 两个人免费观看高清视频| 午夜91福利影院| 99在线视频只有这里精品首页| 黄色女人牲交| 日韩国内少妇激情av| 中文字幕人妻熟女乱码| 免费av毛片视频| 国产精品乱码一区二三区的特点 | 女人爽到高潮嗷嗷叫在线视频| 波多野结衣一区麻豆| 国产精品久久久人人做人人爽| 国产精品久久久久成人av| 99久久国产精品久久久| 深夜精品福利| 亚洲欧美精品综合一区二区三区|