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

    Effects of Drag Coefficients on Surface Heat Flux during Typhoon Kalmaegi (2014)

    2022-08-13 04:47:00LeiLIUGuihuaWANGZeZHANGandHuizanWANG
    Advances in Atmospheric Sciences 2022年9期

    Lei LIU, Guihua WANG, Ze ZHANG, and Huizan WANG

    1Department of Atmosphere and Ocean Sciences, Institute of Atmospheric Sciences,Fudan University, Shanghai 200433, China

    2College of Meteorology and Oceanology, National University of Defense Technology, Changsha 410073, China

    ABSTRACT The lack of in situ observations and the uncertainties of the drag coefficient at high wind speeds result in limited understanding of heat flux through the air-sea interface and thus inaccurate estimation of typhoon intensity in numerical models. In this study, buoy observations and numerical simulations from an air-sea coupled model are used to assess the surface heat flux changes and impacts of the drag coefficient parameterization schemes on its simulations during the passage of Typhoon Kalmaegi (2014). Three drag coefficient schemes, which make the drag coefficient increase, level off,and decrease, respectively, are considered. The air-sea coupled model captured both trajectory and intensity changes better than the atmosphere-only model, though with relatively weaker sea surface cooling (SSC) compared to that captured by buoy observations, which led to relatively higher heat flux and thus a stronger typhoon. Different from previous studies, for a moderate typhoon, the coupled simulation with the increasing drag coefficient scheme outputted an intensity most consistent with the observation because of the strongest SSC, reasonable ratio of latent and sensible heat exchange coefficients, and an obvious reduction in the overestimated surface heat flux among all experiments. Results from sensitivity experiments showed that surface heat flux was significantly determined by the drag coefficient-induced SSC rather than the resulting wind speed changes. Only when SSC differs indistinctively (<0.4°C) between the coupled simulations, heat flux showed a weak positive correlation with the drag coefficient-impacted 10-m wind speed. The drag coefficient also played an important role in decreasing heat flux even a long time after the passage of Kalmaegi because of the continuous upwelling from deeper ocean layers driven by the impacted momentum flux through the air-sea interface.

    Key words:buoy observations,surface heat flux,drag coefficient,numerical simulation,typhoon

    1.Introduction

    Considering the increasing damages made by the landfalls of tropical cyclones (TCs) and the large uncertainties of TC intensity in the warming ocean, it is important to accurately forecast the TC track and intensity. However, compared with the shrinking errors in TC track forecasts, predicted TC intensity lacks striking improvements over recent decades, which could be partly attributed to the limited understanding of energy exchange in the form of surface heat flux(i.e., sensible and latent heat fluxes)(Emanuel, 1986; Jaimes et al., 2015) between TCs and the ocean. Given that the surface heat flux is the primary energy source for TC intensification (Barnes and Dolling, 2013), it is closely correlated with the surface drag coefficient (Cd), which determines not only the momentum transport but also thermal and moisture fluxes through the exchange coefficients (Chand Cq) over the air–sea interface. Unfortunately, lacking knowledge of the processes between TCs and the ocean poses an obstacle to fully understanding the impact of Cdon surface heat flux.Within this context, how exactly Cdbehaves under strong cyclonic winds and what the effects are on surface heat flux and thus the development of TCs remain open questions.

    Due to the destructive winds and dense clouds of TCs,effective observations of complex air–sea interactions in the boundary layer are scarce. Therefore, Cdin high wind speed scenarios (>22 m s?1) is usually measured by extrapolation from low or medium wind speeds (Black et al., 2007). With GPS dropwindsonde data, Powell et al. (2003) revealed that Cdbecame uniform and then slightly decreased at wind speeds greater than 33 m s?1. A similar conclusion was drawn by Donelan et al. (2004) with laboratory extreme wind experiments. Shay and Jacob (2006) further confirmed a “saturation” Cdat 30 m s?1by equating internal wave–ocean fluxes to surface winds at the 10-m level. In addition, Vickery et al. (2009) utilized GPS dropwindsonde data and found that Cdappeared to level off and then decrease at wind speeds greater than 60 m s?1. Most recently, Zou et al. (2018) used two subsurface buoys during Typhoon Megi (2010) to estimate Cdat high wind speeds through two bottom-up methods based on the turbulence closure scheme and bulk model. They found that Cdincreased to a maximum at a critical wind speed of approximately 30 m s?1and then became uniform with wind speed. Although there are no common views of the threshold value at which Cdbegins to decrease, it is widely accepted that Cdcould not increase monotonically with increasing wind speed,which has been parameterized in the numerical models (Bell et al., 2012).

    Presently, the uncertainties of Cdat high wind speeds support the use of three different possible parameterization schemes in numerical models. The first and most widely used Cdtrend is monotonically increasing with wind speed based on the Charnock relationship (Charnock, 1955),though some studies have demonstrated that this trend is not suitable in strong TC scenarios (Green and Zhang, 2013;Zou et al., 2018). The second scheme, which is more consistent with observations (Powell et al., 2003; Donelan et al.,2004), maintains a constant sea surface roughness under high winds up to 40 m s?1. Such a scheme is currently predominant in TC simulations, where wind speeds reach up to 70 m s?1(Drennan et al., 2007; Jeong et al., 2012; Richter and Stern, 2014). The third scheme produces a reduced Cdat high wind speeds. Considering the effects of the wavy ocean surface, Liu et al. (2011) implemented the sea-sprayaffected aerodynamic roughness (z0, which is equivalent to Cd) in a coupled atmosphere–wave–ocean model and further confirmed that TCs weakened with the increasing Cdat lowto-medium wind speeds but intensified in high-wind conditions when Cddecreased. Taking the feedbacks of sea surface temperature (SST) and near-surface wind into consideration,whether this scheme adequately represents Cdin TC conditions remains to be investigated.

    Limited by insufficient in situ observations and the uncertainties in Cdunder high-wind conditions, investigations on the sensitivity of surface heat flux or TC intensity to Cdare still in their nascent stages (Montgomery et al., 2010; Green and Zhang, 2013; Smith et al., 2014). Using an axisymmetric numerical model with balanced boundary-layer dynamics,Emanuel (1986) presented a highly influential perspective wherein the TC intensity is proportional to the square root of the ratio of surface exchange coefficients of enthalpy, Ck,and momentum, Cd. With the first direct measurements of enthalpy flux in a hurricane boundary layer, Zhang et al.(2008) found that there was no statistically significant dependence of these bulk exchange coefficients on wind speed.They concluded that the enthalpy flux into the hurricane boundary layer, which is required to initiate and sustain hurricane development, may have to come from sources other than air–sea turbulent fluxes. Based on their conclusion, the Emanuel model assumptions should be revisited. Smith and Montgomery (2010) suggested that the increase of Cdaround an inertial regime leads to enhanced radial inflow and, thus, significantly intensified surface heat transfer.Zhang et al. (2015) also suggested that enhancing surface friction would enhance boundary layer inflow, which in turn strengthens the convergence of angular momentum and intensifies a TC. However, this could not lead to a positive relationship between Cdand the intensification rate of a TC since surface friction also enhances the momentum and heat dissipation to boundary layer winds, thus dynamically causing a negative effect on TC intensification (Stern et al., 2015; Heng and Wang, 2016; Wang and Heng, 2016). Montgomery et al.(2010) suggested a Cdof 2.0 × 10?3to be the threshold value determining the TC intensification rate and mature intensity in a cloud-resolved numerical model. With a Cdof less than 2.0 ×10?3, both intensification rate and mature intensity of the simulated TC increased slightly with increasing Cd, but the mature intensity decreased with either larger or zero Cd. However, some previous studies argued that both simulated TC intensification rate and mature intensity were insensitive to randomly perturbed Cdwith changes of over 60% (Thomsen et al., 2014; Peng et al., 2018). Chen et al.(2018) suggested the effect of Cdon TC evolution was negligible because Cdinduced moderate sea-surface cooling(SSC) despite the dominant role of Chand Cqin surface heat fluxes. Therefore, the effects of the uncertainties in Cdon both surface heat fluxes and TC intensity need to be better understood.

    Observations from five cross-shaped arrays of buoys and moorings in the northern South China Sea (SCS, Fig. 1)precisely captured the intensity change concerning the ocean responses and associated air–sea interactions of typhoon Kalmaegi (2014) (Zhang et al., 2018). Kalmaegi crossed the northeastern end of Luzon and entered the SCS at 1800 UTC 14 September with a speed of 30 km h?1, heading towards Hainan Island. It passed right over the array of buoys and moorings from 0000 UTC 15 September to 0000 UTC 16 September and thus provided us the opportunity to study the distribution and variation of heat flux under typhoon conditions. In tandem with the high-resolution air–sea coupled model, a series of sensitivity experiments are conducted to examine the effects of Cdon surface heat flux distribution and variation and the feedbacks on typhoon evolution. The rest of the paper is organized as follows: The data, detailed model configurations, and experimental designs used in this study are introduced in section 2. The model output is verified with both conventional and buoy observations in section 3. The sensitivity of the surface heat flux to both SSC and Cdis analyzed in section 4. Finally, conclusions and a brief discussion are given in section 5.

    Fig. 1. Comparisons between the best track data and simulated results of Kalmaegi (2014). The positions of the typhoon at 0000 UTC on each date are marked by circles and squares in (a), while five buoys/moorings are numbered and marked by triangles. (b) shows the time series of the maximum 10-m wind speed of Kalmaegi (2014) from 0600 UTC 12 to 0000 UTC 17 September.

    2.Data and experimental design

    2.1.Data

    Two sets of 6-hourly best-track data for Kalmaegi are obtained from the Joint Typhoon Warning Center (JTWC)and Shanghai Typhoon Institute, China Meteorological

    Administration (STI/CMA), including the location (in longitude and latitude) of the TC center and the maximum sustained 10-m wind speed. The daily SST data at 9-km resolution are obtained from the optimum interpolated SST(OI_SST) of Remote Sensing Systems.

    The in situ observations are from an array containing five cross-shaped buoys and four subsurface moorings(Fig. 1). Stations 1 and 4 are located to the right of Typhoon Kalmaegi (2014), while stations 2 and 5 are to the left. A detailed description of the observational design can be found in Zhang et al. (2016). Since the temperature data at station 1 was lost and the buoy at station 3 was adrift during the passage of Kalmaegi around 0000 UTC 15 September,analyses are mainly based on the data of stations 2, 4, and 5.The in situ observations include: (1) 10-m wind speed,2-m air temperature, and SST to verify the coupled model outputs; (2) subsurface temperature down to 40 m to estimate the ocean heat content and the thermal structure in the upper layer of the ocean; and (3) near-surface humidity to calculate the latent heat flux based on three Cdschemes.

    2.2.Model configuration

    The Weather Research and Forecasting (version 3.7.1 of WRF_ARW, hereinafter referred to as WRF) model (Skamarock et al., 2008) and the Stony Brook parallel ocean model (sbPOM) (Jordi and Wang, 2012), representing the atmospheric and oceanic components, respectively, were coupled to each other through the Model Coupling Toolkit(MCT) (Jacob et al., 2005; Larson et al., 2005). As a fully compressible, non-hydrostatic primitive equation atmospheric model in terrain-following vertical coordinates,WRF is widely used in regional mesoscale atmospheric modeling and operational numerical weather forecasts. The physical parameterization schemes in WRF are the same as those in (Zhang et al., 2019). WRF has 51 vertical half-sigma levels to the height of 50 hPa. There are 570×652 grid points in the x and y directions, respectively, with grid spacing of 9 km, and the outermost domain covered the northwest Pacific (1°?49°N, 100°–149°E, in Fig. A1.). The atmospheric initial and lateral boundary fields adopt the 6-hourly National Centers for Environmental Prediction (NCEP)Final (FNL) operational global analysis data on 1-degree by 1-degree grids.

    In the coupled model, ocean feedbacks are from the Stony Brook Parallel Ocean Model (sbPOM), a paralleled primitive equations ocean model with free surface and sigma coordinate. Vertically, 40 full-sigma levels are configured with a maximum depth of 5000 m. Therein, 16 levels are set in the upper 200 m to better resolve the ocean mixed layer. To reduce the computational cost and avoid unnecessary interpolation error, sbPOM is set with the same domain and horizontal resolution as the WRF model. The initial and lateral boundary conditions of temperature, salinity, currents, and sea level are acquired from the Hybrid Coordinate Ocean Model (HYCOM) + The Navy Coupled Ocean Data Assimilation (NCODA) system Global 1o/12 analysis,which assimilates the satellite altimeter observations and in situ SST, as well as the in situ vertical temperature and salinity profiles from the Expendable bathythermographs (XBTs),ARGO floats, and moored buoys.

    During the coupled simulations, the time steps of atmospheric and oceanic components are set to 30 s and 180 s,respectively. The WRF and sbPOM models interact every 1800 s, which means that the WRF model is integrated over every 60 steps to make it coincide with the sbPOM model after every 10 integration steps. The atmosphere component transfers atmospheric forces (including longwave and shortwave radiation flux, sensible heat flux, latent heat flux, rainfall, and wind stress) to the ocean. Simultaneously, SST and current feedbacks are transferred to the atmospheric component to adjust the bottom conditions and the wind-current shear (Liu et al., 2013). The specific model settings are listed in Table 2.

    2.3.Experimental design

    Theoretically, Cdis a function of wind speed, sea surface state, and atmospheric stability, defined as:

    where k is the von Kármán constant and zrefis the reference height (often 10 m). It not only impacts the amount of momentum flux transferred into water columns but also determines the heat-flux exchange coefficients in air–sea interactions(Donelan et al., 2004):

    where zhand zqare the thermal roughness length and the moisture roughness length, respectively.

    To involve all the Cdbehaviors mentioned above, sensitivity experiments are designed using three main schemes(Fig. 2). The first two Cdparameterization schemes are chosen to represent the increasing and leveling off Cdwith TC wind speed, which are listed in the WRF name-list file as the parameter isftcflx= 0 and 1, respectively. Option 0 is inherited from the Mesoscale Model 5 (MM5) model and yields a monotonic increase of Cdwith wind speed(Charnock, 1955). It is expressed as a function of z0and has been widely used in earlier simulations of TCs with the MM5 model (hereinafter referred to as Charnock55):

    where u?is frictional velocity and g is the acceleration due to gravity. This scheme is called S0, indicating a rougher sea surface with increasing wind speed.

    When the parameter isftcflx is set as 1 (the green line in Fig. 2), Cdincreases with increasing wind speed but saturates at the speed of 33 m s?1, which is more consistent with observations (Powell et al., 2003; Donelan et al., 2004). The momentum roughness length z0is defined as:

    where

    This scheme is called S1, and all other settings are the same as S0.

    Considering the effect of sea spray in high-wind conditions, the third scheme used in this study was proposed by Liu and Guan (2007) (hereinafter referred to as Liu07 and shown by the blue line in Fig. 2). In this scheme, when 10-m wind speed is below 25 m s–1,z0andCdincrease steadily with increasing wind speed. Subsequently,z0gradually reaches its peak value at a wind speed of around 33 m s–1.From thereon,z0(same tendency asCd) becomes uniform or decreases with increasing wind speed. This scheme is given by:

    where α?=cp/u?, ω =min(1,αcr/κu?), and αcr=0.64ms?1.α?denotes the wave age, andcpdenotes the heat capacity at constant pressure for dry air. This scheme is called S2.

    Denoted by a black dashed line in Fig. 2a, both thermal and moisture roughness length (zhandzq, respectively) are set to 10–4to guarantee the independent effects ofCdon changes of both sensible heatChand latent heatCqfluxes[Eqs. (2) and (3)]; i. e.,Chequates toCqin the sameCdscheme. To compare the relative effects ofCdand SSC on surface heat fluxes, six experiments using threeCdschemes respectively are designed and classified into two groups(Table 1). One group, using the atmosphere-only model, is identified with the letter U. Another group is identified with the letter C, denoting the air–sea coupled experiments. All experiments share the same domains and the same atmospheric conditions, except for theCdscheme and ocean responses. Experiments US0 and CS0 adopt the S0Cdscheme, while US1 and CS1 are conducted with the S1Cdscheme, correspondingly. These two groups are designed in order to compare their results with those of past studies(Chen et al., 2018). To analyze the effects of the sprayaffected sea surface on the heat fluxes (Fig. 2), experiments using spray-involving scheme S2 are conducted with both WRF-only and coupled models and are named US2 and CS2, respectively. All the simulations are integrated for 120 hours from 0000 UTC 12 September to 0000 UTC 17 September, covering the period during which Typhoon Kalmaegi(2014) quickly developed and passed through five crossshaped buoy and mooring arrays (Fig. 1a).

    Table 1. List of Cd schemes and experimental designs.

    Table 2. Summary of the model settings.

    3.Typhoon Kalmaegi (2014) and control simulations

    3.1.Typhoon Kalmaegi (2014)

    Kalmaegi (2014) formed in favorable environmental conditions over the Northwest Pacific Ocean and developedinto a tropical storm on 12 September. It then took a northwestward trajectory to Luzon Island and became a Category 1 typhoon on 14 September, peaking at an intensity of 36 m s?1(Fig. 1a). Though Kalmaegi (2014) weaken0.8ed slightly during its first landfall over the Philippines, it quickly regained an intensity of 38 m s?1after entering the SCS on 15 September and began heading toward the southeast coast of China. As a fast-moving typhoon (defined as having a translation speed of > 5.6 m s?1at a latitude between 15°N and 20°N ) (Holland, 1993; Lin et al., 2009), Kalmaegi(2014) passed the mooring arrays with a translation speed of above 8.3 m s?1, with its maximum wind velocity increasing gradually to around 42 m s?1in the early morning of 16 September. Regarding the maximum 10-m wind speed(VMAX, Fig. 1b), some discrepancies exist during this period between the two best track data sets. The VMAX values from CMA maintained 42 m s?1after a quick intensification and were generally stronger than the VMAX values from JTWC. After making landfall in Guangxi Province,Kalmaegi (2014) decayed quickly (not shown).

    3.2.Verification of the simulations

    The recorded tracks of Typhoon Kalmaegi (2014) from JTWC and CMA are largely overlapping, especially when Kalmaegi (2014) enters the SCS (Fig. 1a). Like the best track data, all the simulated tracks pass the array of mooring buoys with a relatively slower translation speed than that seen in the observations, but this does not affect our findings. The overlapping tracks also suggest that the track of Kalmaegi (2014) is mainly determined by the steering flow and environmental systems and is not sensitive to SSC andCd(Green and Zhang, 2013; Chen et al., 2018).

    Fig. 2. Plots of (a) aerodynamic roughness z0, (b) exchange coefficient of drag Cd ,sensible heat Ch and water vapor Cq,and (c) exchange coefficient ratio Ch/C d associated with 10-m wind speed for S0, S1 and S2.

    Similar to what is seen in the best track data, the VMAX results of the six experiments reveal two different stages during the lifetime of Kalmaegi (2014) (Fig. 1b).Before 1800 UTC 14 September, all simulated TCs intensify steadily until Typhoon Kalmaegi (2014) makes its first landfall. After entering the SCS (around 0000 UTC 15 September), Kalmaegi (2014) resumes intensifying in all simulations, reaching its maximum intensity at 0000 UTC 16 September. In general, the coupled model reproduces weaker maximum intensities than the WRF-only simulations, except for those using the S1 scheme; CS1 and US1 simulate similar wind speeds (shown as the green lines in Fig.1b). Compared with other experiments, the coupled simulation with the Charnock55 scheme (CS0) produces the proximal intensity to the best tracks.

    To further validate the coupled model and interpret the effect of Cdon ocean response, SSTs from the coupled simulations are compared with the WindSat SSTs (Wentz et al.,2013) shown in Fig. 3. Evident cooling appears to the right of the track, with a minimum value of 26°C in the WindSat data, which is 1°C–2°C colder than the surrounding waters.Though all coupled simulations capture SSC to the right of the track, both magnitude and coverage of the SSC are weaker than the observation. This could also be a reason for a simulated TC that is stronger than the best track data suggest. Discrepancies also exist between the three coupled experiments, in which SSC in CS0 is more evident and consistent with the WindSat data, indicating stronger ocean responses with increasing Cd.

    Fig. 3. SST (shaded) at 0000 UTC 17 September and typhoon track (black solid lines) from (a) CS0, (b) CS1, and (c)CS2, and (d) SST (shaded) at 0000 UTC 17 September from Windsat and typhoon tracks from CMA (thick solid line)and JTWC (thin solid line), respectively. Black triangles denote the positions of the five buoys/moorings.

    3.3.Comparisons with the buoy observations

    The simulated 10-m wind speed and differences between SST and 2-m air temperature are compared with the buoy observations in Figs. 4a1–a3 and Figs. 4b1–b3,respectively. To reduce the instantaneous error in the comparisons, mean values from the nine model grid cells nearest to the buoy location are adopted. Since the sizes of the simulated TCs (represented by the radius of maximum 10-m wind speed) are similar (not shown), their effects are not considered here. For 15 September, the coupled model reproduces a tendency of the 10-m wind similar to that of the buoy observations, though in relatively higher values. To the right of the track, the buoy at station 4 captured a VMAX of about 22 m s–1at 1200 UTC 15 September, which is much stronger than that at stations 2 and 5. The peak value at station 2 occurred at 0000 UTC 16 September because of its further west location and the intensifying typhoon. Comparing all simulations, the 10-m wind speed of CS0 is the weakest due to the increasing Cd.

    For the temperature differences (? T) between 2-m height and the sea surface, relatively large discrepancies exist among model results and buoy observations. At stations 2 and 5, the modeled ?Ts are similar and steady over time with a minimum value of –1°C, whereas the observed ?T decreased rapidly to –2.5°C on 15 September and varied significantly afterwards. This suggests that the near-surface air cooling, mostly induced by heat dissipations, evaporation,and heavy rainfall under the eyewall, could possibly be underestimated in the numerical model, especially when SSC is indistinctive under a fast-moving typhoon. However, ? T tendencies at station 4 were different from those at stations 2 and 5, which were to the left of the track. Positive ? T values of about 1.5°C were observed throughout this stage. The observed values are about 1.0°C higher than the results of experiments CS0 and CS2 and are even higher than the results of experiment CS1. This could be attributed to the underestimation of the SSC relative to the observations when the wind speed is weaker than 30 m s?1(Fig. 3), especially when using the leveling off Cdscheme. In other words, SSC plays the key role in determining ?T under typhoon conditions on the right side of the TC but is most likely underestimated by the air–sea coupled model, even with the increasing Cdscheme.

    Fig. 4. Time series of 10-m wind speed (a1–a3), temperature differences between 2 m and sea surface (b1–b3), and latent heat flux (c1–c3) from the coupled simulation cases (CS0–2) and observational data (Obs) at stations 2, 4, and 5 from 0000 UTC 15 September to 0000 UTC 17 September. For (c1–c3), bold lines indicate coupled-model results (CS0–2), and buoy observations calculated with the three Cd schemes (S0–2, named Obs_S0, Obs_S1, and Obs_S2, respectively) are shown in light-colored solid lines.

    Latent heat flux (LHF), known as the primary energy source for TCs, is also investigated (Figs. 4c1–c3). LHF is written as:

    where ρ is the air density and Lvdenotes the specific heat capacity of air at constant pressure. ?q=(qs?qa), where qais the specific humidity at a reference height of 2 m above the sea surface and qsis the saturation humidity at the sea surface. Uais the wind velocity at the same level of qa. Since the buoys could not measure LHF directly, the changes of which was calculated with 10-m wind speed and δq from the buoy observations according to Eq. (10) while Cqestimated from the parameterization schemes. As shown in Figs. 4c1–c3, the calculated LHFs with the three Cqschemes are almost overlapping, indicating the insensitivity of LHF to Cq. The modeled LHFs exhibit synchronous changes but are higher than the observations, especially at stations 2 and 5, which is possibly due to the relatively stronger 10-m wind speeds (Figs. 4a1–a3) and higher qaover the warmer sea surface. The modeled LHF at station 4 reaches a maximum value that is 100 W m?1higher and about 9 h earlier than the buoy observations and then quickly decreases after 1800 UTC 15 September. The LHF of CS0 shows more consistencies with the observations due to the stronger SSC compared to the other experiments. The above analyses also suggest that the Cd-induced ocean responses associated with LHF could be different on both sides of the typhoon track. Supported by earlier work, “right bias” is a typical phenomenon for fast-moving TCs in which wind vectors turn clockwise and strengthen the wind velocity to the right of the track. Therefore, the ocean responses could be enlarged due to the uncertainty of Cd, resulting in evident discrepancies in LHF.

    4.Sensitivity of surface heat flux to the drag coefficient

    Determining both momentum (e.g., wind stress) and enthalpy flux through the air–sea interface, Cdnot only directly changes near-surface wind speed, but also drives the mixing and upwelling processes in the ocean upper layer and thus the SSC under TCs. Both SSC and wind speed are known as the main reasons for changing surface heat flux.Therefore, the relative effects of the Cd-induced SSC and wind speed on the surface heat flux should be examined.

    From the very beginning, the girl s family objected strongly to her dating this guy, saying that it had got to do with family background and that the girl would have to suffer for the rest of her life if she were to be with him.

    4.1.Effects of Cd-induced SSC and wind speed on heat flux

    Since CS1 and US1 simulate similar wind speeds(Fig. 1b) for 15 September, the daily-averaged sensible and latent heat fluxes from these two experiments are firstly compared with the NCEP data to verify the effects of SSC(Fig. 5). Although the heat fluxes from the NCEP data are characterized by homogeneous distributions due to the relatively low resolution (1° × 1°), both magnitude and pattern of the data are more consistent with the results from CS1.US1, with fixed SST, overestimates the heat fluxes, especially around the five mooring buoys, further emphasizing the importance of ocean feedbacks in TC simulation (Elsberry et al., 1976; Price, 1981; Shay et al., 1992).

    Fig. 5. Daily averaged sensible (upper row) and latent heat fluxes (bottom row) from (a, d) CS1, (b, e) US1, and (c, f) NCEP(shaded, units: W m?2) on 15 September. Black triangles denote the positions of the five buoys/moorings.

    Fig. 6. Daily averaged (a) sensible heat flux differences (shaded, units: W m?2) with SST differences (blue lines,units: °C) and (b) latent heat flux differences (shaded, units: W m?2) with 10-m wind speed differences (black lines,units: m s?1) between US1 and CS1 on 15 September. Black triangles are the positions of the five buoys/moorings.Red squares are the positions of typhoon centers in CS1.

    Figure 6 shows the differences in sensible and latent heat fluxes (S/LHD), SST, and 10-m wind speed between US1 and CS1. The S/LHD are characterized by negative values, especially to the right of the track, which is consistent with the remarkable SST cooling of –0.6°C to –1.0°C and slightly weakened wind speed (< –1 m s?1) in CS1. Positive S/LHD are located where SSC is smaller than –0.4°C, showing an indistinctive relationship with the 10-m wind speed(about 0.4 m s?1, not shown). Note that the variation of LHF is three times larger than that of SHF, although their patterns are similar. LHF is changed significantly by SSC (Wada et al., 2018), which determines the intensification of Typhoon Kalmaegi (2014) in the following stages.

    Figure 7 compares the different effects from Cd-induced SSC and wind speed on the LHD at stations 2, 4,and 5. Shown by the red lines, the differences between the SST-fixed experiments US0 and US1, US1 and US2, respectively, denote the effect of Cdon LHD through wind changes. Similarly, comparisons of C0 minus C1 and C1 minus C2, respectively, reveal the effect of different Cdschemes on LHD. Three green lines show the difference of LHD between the coupled simulations and the uncoupled ones, including C0 (CS0 minus US0), C1 (CS1 minus US1),and C2 (CS2 minus US2), respectively, indicating the impacts of Cd-induced SSC on LHD. At station 4, significant LHD (shown by green lines in Fig. 7) appears at around 1200 UTC 15 September with a magnitude much larger than what could be the effects from Cd-induced wind speed and different Cdschemes, suggesting that the SSC, relating to the ocean status, is the predominant factor determining the air–sea surface heat fluxes. The changes of LHF through Cd-impacted wind speed is almost equivalent to those induced by different Cdschemes, indicating that the near-surface wind speed only has a significant effect when SSC is inconspicuous. The results of Jaimes et al. (2015) also support the finding that the enthalpy fluxes are closely related to the upper oceanic thermal structure but independent of wind speed, which could be altered by the pressure field when SSC occurs. Compared with station 4, station 2 was located in the front left of the typhoon. Except for the strengthening SSC under the typhoon, cold water to the right rear of Kalmaegi (2014) was transported cyclonically to station 2,resulting in obvious LHD 12 hours later than at station 4.Located to the rear left of the typhoon, SSC at station 5 was not significant, and the LHD induced by SSC and wind speed was therefore of equal magnitude.

    Fig. 7. Time series of differences in latent heat fluxes (LHD)from Cd-induced wind speed (US0–US1; US1–US2), SSC (C0=CS0–US0; C1=CS1–US1; C2=CS2–US2), and SSC change(C0–C1; C1–C2) on stations 2, 4, and 5 from 0000 UTC 12 September to 0000 UTC 17 September.

    4.2.Effects of different C d schemes on heat flux

    As shown by the blue lines in Fig. 7, the LHD induced by the three Cdparameterization schemes in the coupled model are of similar magnitude as the wind-impacted LHD(shown by red lines), suggesting there is less of an effect by different Cdschemes on enthalpy flux than by wind-driven evident SSC. To further understand the sensitivity of enthalpy flux under the storm to the three Cdschemes, differences in daily surface heat flux between CS experiments are compared in Fig. 8.

    During the passage of typhoon Kalmaegi on 15 September, SST differs slightly between the coupled simulations(less than –0.2°C) to the right of the track, whereas the difference of 10-m wind speed maintains –1 m s?1to –3 m s?1(Figs. 8c and d), at a similar magnitude as that in Fig. 6b. Different from the pattern seen in Fig. 6, the weakened 10-m wind overlaps negative S/LHD zones randomly on both sides of the track rather than to the right of it, supporting the findings of the wind effects on LHF with weak SSC suggested in section 4.1.

    Compared with Figs. 8b and d, the negative S/LHD zones of CS0 and CS1 (Figs. 8a and c) are more notable,which could be explained by the differences of Cdvalues and their impacts during this period. Before 1200 UTC 15 September, the VMAX of CS0 (Fig. 1b) accelerates to approximately 33 m s?1, with Cdbeing around 3 × 10?3(Fig. 2).Simultaneously, Cd(or Ch/Cd) in CS1 gradually rises but is limited to 2.5 × 10?3(or Ch/Cdat 0.7). The difference of Cd(or Ch/Cd) between CS0 and CS1 thus enlarges with the increasing 10-m wind speed (Figs. 2b and c) to about 0.9 ×10?3at 0000 UTC on 16 September. In contrast, the difference of Cd(or Ch/Cd) between CS1 and CS2 decreases first with increasing wind speed and then increases moderately to 0.6 × 10?3, being 0.3 × 10?3smaller than the difference of Cdbetween CS0 and CS1. Since the increased surface friction could decelerate the near-surface wind speed immediately and thus inhibit the enthalpy flux, which is much faster than the ocean cooling, smaller S/LHD show up, suggesting high sensitivity of surface heat fluxes to Cdschemes.

    During the early period when the VMAX is still weaker than 30 m s?1, positive LHD between CS0 and CS1 show up in front of the TC (Fig. 9a). As Kalmaegi (2014) intensifies,Cdincreases quickly with increasing wind speed, especially for CS0. LHD between CS0 and CS1 exhibit negative values in most areas east of the track, supporting the above findings that LHF relates closely to both wind speed and Cdwhen SST changes moderately. In contrast, evident SSC dominates LHD between CS1 and US1 (Fig. 9b), which are negative in all directions, though both TCs are of similar intensity, especially to the northeast of the center. Moreover, LHD induced by Cdparameterization scheme differences (Fig. 9a)are one-third or less than LHD caused by SSC (Fig. 9b), suggesting the wind-induced SSC is the more dominant factor for LHF compared to wind speed.

    4.3.Effects of Cd on the thermal structure of the ocean mixed layer

    The downward momentum transport through an air–sea interface, determined by Cdand wind speed, usually induces cooling of the ocean mixed layer (OML) when the SST becomes 0.5°C colder than the water below (e.g.,Black et al. 2007) and vertical mixing and upwelling of deeper water occurs. Following Zhang et al. (2016), sea temperature of the upper 40 m is extracted from all coupled simulations and compared with the buoy observation at station 4.Figures 10a?d present the temperature differences between run time and 0000 UTC 15 September. Time series of temperature in the ocean upper 160 m are preseted in Fig. A2 .Before 0600 UTC 15 September, the simulated sea temperature below 20 m increases by about 0.5°C, while the surface layer cools slightly, indicating weak vertical mixing in the shallow OML. Later on, evident cooling occurs in the OML as the typhoon intensifies. As shown by Fig. 2b, although the wind speed of CS0 is weakest at station 4,Cdis largest there throughout 15 September. Contrarily, Cdin CS1 is the smallest at station 4, with the wind speed being weaker than 33 m s?1, within the same period. The cooling OML(Figs. 10a–c) is proportionate to Cdon 15 September for both CS0 and CS1, and CS0 produces the most significant cooling with the longest duration, while CS1 exhibits the weakest effects. Compared with the OML pattern at station 4, the magnitude of the cooling in CS0 is most consistent with the buoy observation. Note that the observed maximum cooling occurred at the depth of 20 m, different from the sea surface in a coupled model. This could be attributed to the incorrectly resolved solar heating in the model or horizontal diffusion caused by mesoscale vortexes and highlights the urgent need in developing an advanced air–sea coupled model. The variation of ocean heat capacity (OHC) with time could further reflect the effect of Cdon the thermal status of the OML (Fig. 10e) (Jaimes et al., 2015). The model outputs a similar tendency with, but relatively higher than, the observed ?OHC, which could be caused by the existing cold eddies observed by the altimeter at station 4 (not shown). In CS0, OHC decreases gradually by 100 KJ cm?2after 0600 UTC 15 September, while the OHC of CS1 simultaneously decreases only by about 50 KJ cm?2.

    Fig. 9. Azimuth–time Hovm?ller diagrams of the normalized differences of latent heat flux between (a) CS0 and CS1 and (b) CS1 and US1 (shaded, units: W m?2). Normalization is done by the azimuthally averaged latent heat fluxes in CS1 from Rmax to 4Rmax over the period 0000 UTC 15 September to 0000 UTC 17 September. Note the different label bars of the two snapshots.

    To further elucidate the effects of Cdon the slow recovery ocean responses and the surface heat fluxes evolution, the same analyses shown in Fig. 8 are made for 16 September and shown in Fig. 11. Since the maximum intensity of Kalmaegi (2014) on 0000 UTC 16 September varies greatly between the three experiments, the differences of SST, wind field, and enthalpy flux between the CSs are thus significant near 114°E, 19°N. The notable positive differences of SST,combined with the heat flux differences between CS1 and CS2 to the right of the center (Fig. 11b) are related to the quickly decreasing Cdwith peak wind speed of higher than 40 m s?1. Across the mooring array, however, the daily averaged 10-m wind speeds differ little, while the SST differences among the three experiments remain larger than 0.3°C (Figs.11a and b) and are well correlated with the variations of heat fluxes. Therefore, the distribution of surface heat flux could be greatly affected by choice of Cdparameterization scheme through the longer lasting SSC, even after the passage of a TC.

    Fig. 10. Time series of temperature differences in the ocean upper 40 m from CS0 (a), CS1 (b), CS2 (c), and buoy observations (Obs) at station 4 (d) and ocean heat capacity (OHC) differences (e) between 0000 UTC 17 September and 0000 UTC 15 September.

    Fig. 11. Same as Fig. 8, but for 16 September.

    5.Conclusion and discussion

    Choosing a suitable Cdparameterization scheme is of great importance for numerical models, as Cdnot only controls the surface wind stress directly but also affects the thermal status of the upper ocean layers. With observed Cdand Ckprofiles being incorporated in numerical models, better TC intensity and structure could be simulated, especially over the ocean (Ming and Zhang, 2016). However, the uncertainties of Cdin high-wind scenarios make the estimation of surface heat fluxes in numerical simulations challenging, limiting improvement of typhoon intensity forecasting. Benefiting from buoy observations (Zhang et al., 2016) and the fully-coupled WRF-sbPOM model, a series of Cd-sensitivity numerical simulations were conducted in this study to understand the impacts of Cdon surface heat fluxes during Typhoon Kalmaegi (2014).

    The air–sea coupled model was first confirmed to perform better in simulating the intensity of Kalmaegi with more reasonable SSC than the SST-fixed WRF model by comparing output with two best track datasets and the satelliteobserved SST. Sensitivity experiments with three Cdparameterization schemes produced evident differences in near-surface winds and air–sea heat flux, suggesting the Cdscheme to be an influential factor impacting numerical simulations.Due to the significant SSC associated with an evident reduction in heat flux, the coupled model with the increasing Cdparameterization scheme reproduced the weakest TC, which was most consistent with the best track data. The reason is twofold: 1) as the coupled model tended to overestimate the intensity of a Category 1 typhoon, the increasing Cdscheme could effectively enhance the surface friction and SSC, therefore reducing the intensity; 2) the ratio of Ch/Cdin the CS0 experiment was between 0.5 and 0.65, which had been calculated with the first direct measurements of enthalpy flux in the hurricane boundary layer (Zhang et al., 2008) and proven to be suitable in estimating the air–sea fluxes relative to the OHC changes (Jaimes et al., 2015). Nevertheless, the Ch/Cdratios in CS1 and CS2 were greater than 0.7 as wind accelerated, especially between 15 September and 16 September, when the Ch/Cdratio exceeded 0.8 in CS2. It was significantly greater than the steady value of 0.65 in Zhang et al.(2008). With Cdleveling off or decreasing in the high-wind scenario, SSC was inhibited. OML cooling and heat loss through the air–sea interface was therefore weakened, resulting in a stronger TC.

    Compared with the buoy observations in the SCS, the simulated 10-m wind speeds exhibited a similar tendency but at a relatively higher magnitude than the buoy observations. The ?T between 2-m air temperature and SST was also compared between observations and model outputs. Discrepancies existed on different sides of the TC track. To the left of Kalmaegi (2014), observed ?T declined more obviously than the model results, suggesting that the near-surface air cooling could be likely underestimated in the numerical models. Rather, the significantly higher ?T observed at station 4 indicated the underestimation of SSC to the right of TC track by the coupled model. Due to the warmer initial SST in the WRF-sbPOM model, the simulated LHF was generally higher than buoy observations before 1200 UTC 15 September, especially at stations 2 and 5, but it then decreased quickly to the observed magnitude at station 4 after the appearance of evident SSC.

    As two important factors determining heat flux, SST and surface winds are significantly affected by the surface status denoted by Cdand Ck(or Ch/Cd; Ming and Zhang,2016). But the sensitivity of surface heat flux to these two parameters is different. Compared with the near-surface wind fields, the heat exchange was found to be more sensitive to Cd-induced SSC. With intense SSC, LHF was significantly inhibited around the core region of the typhoon, contributing to a weaker TC. However, when SST differed inconspicuously among the coupled simulations, the heat fluxes tended to be positively associated with the change of 10-m wind speeds. Therefore, different Cd(or Ch/Cd) schemes won’t significantly affect the enthalpy fluxes through wind-induced processes, unless SSC is indistinctive. Wadler et al. (2021)and Zhang et al. (2015) also emphasized the importance of measuring atmospheric and oceanographic parameters in understanding TC intensity and structure changes. Moreover, the Cdschemes examined here also affected the distributions of heat flux through continuous cooling of the deeper ocean, even long after the passage of TC. A larger Cdpromoted more momentum to be transported into the water columns, resulting in the continuous upwelling of the deeper colder water to the upper layer and the cooling of the OML (Fig. A3). This process was more evident to the right of the typhoon, consistent with the typically stronger wind field to the right of TCs. Therefore, a suitable Cdscheme is essential for accurate numerical simulation of TCs.

    The above results also indicate the complexity of choosing Cdschemes in numerical simulations, especially for a moderate typhoon, as the threshold value of Cdwhen becoming uniform or decreasing is uncertain. Powell et al. (2003)detected 33 m s?1as the threshold value. Vickery et al.(2009), however, suggested that the threshold value is up to 60 m s?1. Since the intensity of Kalmaegi (2014) was weaker than 60 m s?1, observations under stronger TC conditions were not obtained. Therefore, more numerical simulations and observations on the threshold values of Cdshould be conducted to draw a concrete conclusion on the speed of wind corresponding to Cdleveling off or decreasing.

    Another interesting phenomenon is the different intensification rates of the typhoon between the experiments with(the CSs cases) and without (the US cases) ocean feedbacks from 0000 UTC to 1800 UTC 15 September. Similar to the results found in earlier studies using an atmosphere-only model (Wang and Xu, 2010; Zeng et al., 2010), different Cds in the US experiments rarely changed the intensification rates of the TC (dashed lines in Fig. 1b) but did lead to different maximum intensities of the mature TC. However, in the air–sea coupled model, the intensification rates of Kalmaegi(2014) were different using the three Cdparameterization schemes when the wind speed became 35 m s?1or stronger.Both VMAX and the intensification rate in CS2 are undoubtedly the largest because of the decreased Cdin high-wind conditions, while VMAX and the intensification rate in CS0 were the smallest due to the increasing sea surface roughness throughout the simulation. The results found here are different from earlier findings when using single atmospheric models(Wang and Xu, 2010; Zeng et al., 2010). This could be attributed to the enhanced enthalpy fluxes with complex air–sea conditions in the boundary layer. Suggested in previous studies (Zhang et al., 2017; Wadler et al., 2018; Zhang and Rogers, 2019), the surface heat flux-induced boundary layer recovery, which has been proven to regulate the thermodynamic structure of the boundary layer, plays an important role in TC development, including enhancing both the inflow and convergence and resisting the weak and dry downdraft near the core region, which in turn leads to stronger and more symmetric deep convection in TCs. Therefore, the role of SSC-impacted heat flux should be considered in future numerical simulations.

    Fig. A1. Configuration of the coupled model domain. The color shades denote land and ocean topography in meters.

    Fig. A2. Time series of temperature in the ocean upper 160 m.

    Fig. A3. The schematic diagram of the processes between the wind-induced SSC and its feedbacks throughout the sea surface with different Cd schemes.

    As is suggested in earlier studies, wave breaking and sea spray could significantly change the wave-current fields and enthalpy fluxes at high wind speeds, which are closely related to Cd(Andreas and Emanuel, 2001; Moon et al.,2004; Bao et al., 2011). Therefore, a better understanding of the oceanic effects on Cdduring the early stage of a typhoon is needed to develop a suitable Cdparameterization scheme in the coupled model. Besides, the discrepancies between in situ observations and model output also highlight the urgent need for advanced air–sea coupled models and more observation facilities in the future.

    Acknowledgements.The authors are grateful to two anonymous reviewers for their constructive comments on this work and recommending several beneficial references. The authors thank Prof. Dake CHEN for providing the in situ mooring data, and Prof.Weimin ZHANG and Jianfang FEI for their help in improving the manuscript. This study has been supported by the National Natural Science Foundation of China under Grant Nos. 41775053,41976003, and 42192552, and the National Key Research and Development Program of China under Grant Nos. 2019YFC1510001 and 2019YFC1510102. Additional support has been provided by the National Program on Global Change and Air-Sea Interaction(GASI-IPOVAI-04). We acknowledge model support from http://www2.mmm.ucar.edu/wrf/users/ and http://imedea.uib-csic.es/users/toni/sbpom/, provided by NCAR’s Antoni JORDI and Dong-Ping WANG, respectively. The FNL data for the WRF model was downloaded from https://rda.ucar.edu/datasets/ds083.2/. The ocean topographic data, ETOPO1, was acquired from the one-arc-minute global relief model of Earth’s surface (https://maps.ngdc.noaa.gov/viewers/wcs-client/). The ocean field for sbPOM was from https://www.HYCOM.org/data/glbu0pt08/expt-91pt1.

    APPENDIX

    嫩草影院精品99| 国精品久久久久久国模美| 国产高清不卡午夜福利| 国产乱人视频| 亚洲熟女精品中文字幕| 少妇的逼好多水| 成人午夜高清在线视频| 免费黄网站久久成人精品| 日韩欧美国产在线观看| 国产精品99久久久久久久久| 国产爱豆传媒在线观看| 人人妻人人澡人人爽人人夜夜 | 啦啦啦韩国在线观看视频| kizo精华| 国产一级毛片七仙女欲春2| 国内揄拍国产精品人妻在线| 舔av片在线| 18禁在线无遮挡免费观看视频| 最近2019中文字幕mv第一页| 国产精品一区www在线观看| 国产精品伦人一区二区| a级一级毛片免费在线观看| 在线a可以看的网站| 久久久久久九九精品二区国产| 国产一区二区在线观看日韩| 国产一区亚洲一区在线观看| 国产老妇伦熟女老妇高清| 春色校园在线视频观看| 日韩 亚洲 欧美在线| 亚洲国产精品成人久久小说| 久久精品国产亚洲网站| 久久久久久久久中文| 国产精品一区二区在线观看99 | 午夜久久久久精精品| 久久这里有精品视频免费| 成人毛片60女人毛片免费| 亚洲成人中文字幕在线播放| 精品一区在线观看国产| 亚洲国产最新在线播放| 欧美性猛交╳xxx乱大交人| 亚洲欧美精品自产自拍| 人妻制服诱惑在线中文字幕| videos熟女内射| 国产亚洲av嫩草精品影院| 国产女主播在线喷水免费视频网站 | 亚洲av电影在线观看一区二区三区 | 乱人视频在线观看| 亚洲成人精品中文字幕电影| 三级经典国产精品| 狂野欧美白嫩少妇大欣赏| 久久97久久精品| 超碰av人人做人人爽久久| 国产男人的电影天堂91| 久久久久精品久久久久真实原创| 国产精品99久久久久久久久| 国产亚洲av片在线观看秒播厂 | 久久久久免费精品人妻一区二区| 国内揄拍国产精品人妻在线| 国产精品爽爽va在线观看网站| 亚洲欧美成人综合另类久久久| 亚洲成人av在线免费| 成年女人在线观看亚洲视频 | 日产精品乱码卡一卡2卡三| 色综合亚洲欧美另类图片| 亚洲国产欧美在线一区| 黄片无遮挡物在线观看| 久久久久久久久久成人| 精品一区二区三卡| 色视频www国产| 午夜福利视频1000在线观看| 国产av不卡久久| 国产精品1区2区在线观看.| 丝袜喷水一区| 亚洲成人一二三区av| 国产亚洲精品久久久com| av在线播放精品| 国内精品一区二区在线观看| 国产日韩欧美在线精品| 搞女人的毛片| 亚洲精品,欧美精品| 国产精品一区二区三区四区久久| 亚洲欧美日韩卡通动漫| 欧美丝袜亚洲另类| 91狼人影院| 一个人免费在线观看电影| 国产 一区精品| 亚洲伊人久久精品综合| 男人和女人高潮做爰伦理| 免费少妇av软件| 日韩国内少妇激情av| 久久国产乱子免费精品| 天美传媒精品一区二区| 在线免费观看不下载黄p国产| 国产成人精品一,二区| 岛国毛片在线播放| 欧美日韩精品成人综合77777| 精品一区在线观看国产| 久久久久性生活片| 99久久精品一区二区三区| av.在线天堂| av福利片在线观看| 欧美日韩精品成人综合77777| 色综合站精品国产| 五月天丁香电影| 熟妇人妻不卡中文字幕| 三级经典国产精品| 精品久久久久久久久亚洲| 丰满乱子伦码专区| 99热这里只有是精品在线观看| av专区在线播放| 国产有黄有色有爽视频| 国产不卡一卡二| 中国国产av一级| 丰满少妇做爰视频| 中文字幕制服av| 国产高清三级在线| 国产精品一区www在线观看| 国产探花在线观看一区二区| 午夜激情福利司机影院| 国产精品一二三区在线看| 在线免费十八禁| 能在线免费观看的黄片| 亚洲最大成人手机在线| 久久久色成人| 蜜桃亚洲精品一区二区三区| 久久精品国产亚洲网站| 亚洲aⅴ乱码一区二区在线播放| 国产精品人妻久久久影院| 久久这里只有精品中国| 久久人人爽人人片av| 久久午夜福利片| 免费看美女性在线毛片视频| 一个人观看的视频www高清免费观看| 日本猛色少妇xxxxx猛交久久| 校园人妻丝袜中文字幕| 国产成人福利小说| 九色成人免费人妻av| 亚洲乱码一区二区免费版| 日韩精品有码人妻一区| 亚洲精品成人久久久久久| 亚洲在久久综合| 国产精品久久久久久久电影| 国产精品一区二区性色av| 两个人的视频大全免费| 少妇人妻一区二区三区视频| 1000部很黄的大片| 美女脱内裤让男人舔精品视频| 欧美日本视频| 亚洲不卡免费看| 亚洲国产成人一精品久久久| 18禁动态无遮挡网站| 老女人水多毛片| 99久国产av精品国产电影| 久久久久精品性色| 久久午夜福利片| 国产人妻一区二区三区在| 九九爱精品视频在线观看| 九草在线视频观看| 热99在线观看视频| 卡戴珊不雅视频在线播放| 99热这里只有精品一区| 精品国产露脸久久av麻豆 | 日本色播在线视频| 精品国产三级普通话版| 国产免费视频播放在线视频 | 欧美三级亚洲精品| 精品久久久久久电影网| 777米奇影视久久| 日日撸夜夜添| 日韩大片免费观看网站| a级毛片免费高清观看在线播放| 日本免费a在线| 国产伦理片在线播放av一区| 日本黄色片子视频| 国内精品宾馆在线| 成人亚洲精品av一区二区| 高清毛片免费看| 女人被狂操c到高潮| 国产视频内射| 女人被狂操c到高潮| 亚洲成人久久爱视频| 免费黄色在线免费观看| 免费看日本二区| 直男gayav资源| 久久久久久久久久久丰满| 国产在线男女| 国产伦在线观看视频一区| 黄色配什么色好看| 国产精品av视频在线免费观看| 啦啦啦中文免费视频观看日本| 亚洲精品成人av观看孕妇| 国产伦理片在线播放av一区| 精品一区二区免费观看| 春色校园在线视频观看| 久久这里有精品视频免费| 午夜福利在线在线| 性插视频无遮挡在线免费观看| 免费黄网站久久成人精品| 在线a可以看的网站| 精品午夜福利在线看| 国产精品一区www在线观看| 成人午夜精彩视频在线观看| 国产 一区精品| 一二三四中文在线观看免费高清| 国产高清三级在线| 国产亚洲5aaaaa淫片| 乱人视频在线观看| 蜜桃亚洲精品一区二区三区| 国产熟女欧美一区二区| 久久久成人免费电影| 国产av国产精品国产| 色尼玛亚洲综合影院| 成人亚洲欧美一区二区av| 97超视频在线观看视频| 最近2019中文字幕mv第一页| 亚洲,欧美,日韩| 男人舔女人下体高潮全视频| 欧美成人一区二区免费高清观看| 中文在线观看免费www的网站| 特级一级黄色大片| 黄片无遮挡物在线观看| 啦啦啦中文免费视频观看日本| 国产老妇伦熟女老妇高清| 精品国内亚洲2022精品成人| 国产亚洲5aaaaa淫片| 一夜夜www| 午夜精品国产一区二区电影 | 综合色av麻豆| 99热这里只有是精品在线观看| 搞女人的毛片| av又黄又爽大尺度在线免费看| 国产老妇伦熟女老妇高清| 免费看不卡的av| 国产精品福利在线免费观看| 亚洲欧美日韩卡通动漫| 内地一区二区视频在线| 国产精品一区二区三区四区免费观看| 亚洲一级一片aⅴ在线观看| 性色avwww在线观看| 国产免费又黄又爽又色| 中文字幕av成人在线电影| 超碰av人人做人人爽久久| 亚洲人成网站在线播| 大陆偷拍与自拍| 亚洲av免费在线观看| 在线免费十八禁| 日日摸夜夜添夜夜添av毛片| 91av网一区二区| 亚洲精品色激情综合| 亚洲国产色片| 中文资源天堂在线| 国产精品av视频在线免费观看| 一区二区三区乱码不卡18| 最近最新中文字幕大全电影3| 一区二区三区高清视频在线| 免费看美女性在线毛片视频| 九九爱精品视频在线观看| 欧美xxⅹ黑人| 嫩草影院新地址| 国产在视频线在精品| 大话2 男鬼变身卡| 深爱激情五月婷婷| 国产在视频线在精品| 国内少妇人妻偷人精品xxx网站| 欧美性感艳星| 女的被弄到高潮叫床怎么办| 亚洲国产高清在线一区二区三| 精品久久国产蜜桃| 成人性生交大片免费视频hd| 97超碰精品成人国产| 日韩中字成人| 亚洲国产精品成人久久小说| 少妇高潮的动态图| 女人久久www免费人成看片| 免费av毛片视频| 国产人妻一区二区三区在| 国产探花极品一区二区| 国产精品国产三级国产av玫瑰| 欧美xxxx性猛交bbbb| 国产探花极品一区二区| 色吧在线观看| 亚洲精品第二区| 日韩欧美一区视频在线观看 | 欧美性感艳星| 亚洲不卡免费看| 亚洲电影在线观看av| 久久久精品94久久精品| 久久草成人影院| 日日摸夜夜添夜夜添av毛片| 七月丁香在线播放| 少妇人妻精品综合一区二区| 天天一区二区日本电影三级| 欧美精品国产亚洲| 婷婷色av中文字幕| 亚洲欧美精品自产自拍| 亚洲av电影在线观看一区二区三区 | 精华霜和精华液先用哪个| 禁无遮挡网站| 熟妇人妻久久中文字幕3abv| 最近最新中文字幕免费大全7| 国产人妻一区二区三区在| 99re6热这里在线精品视频| 色吧在线观看| 麻豆av噜噜一区二区三区| 99视频精品全部免费 在线| 国产成人精品婷婷| 免费黄网站久久成人精品| 国产 一区精品| 国产乱来视频区| 久久久亚洲精品成人影院| 国产精品嫩草影院av在线观看| 亚洲经典国产精华液单| 午夜福利视频精品| 亚洲经典国产精华液单| 日韩欧美三级三区| 免费看av在线观看网站| 日韩在线高清观看一区二区三区| 精品一区二区三卡| av专区在线播放| 国产探花在线观看一区二区| 免费人成在线观看视频色| 高清视频免费观看一区二区 | 亚洲国产欧美在线一区| av在线播放精品| 国产单亲对白刺激| 免费人成在线观看视频色| 亚洲av二区三区四区| 少妇高潮的动态图| av在线观看视频网站免费| 亚洲av福利一区| 天天躁日日操中文字幕| 国产美女午夜福利| 直男gayav资源| 日日摸夜夜添夜夜添av毛片| 丰满人妻一区二区三区视频av| 99热全是精品| 寂寞人妻少妇视频99o| 最近视频中文字幕2019在线8| 午夜久久久久精精品| 国产黄色免费在线视频| 99久久精品国产国产毛片| 高清视频免费观看一区二区 | 精品亚洲乱码少妇综合久久| 黄片无遮挡物在线观看| 中文乱码字字幕精品一区二区三区 | 免费黄频网站在线观看国产| 日韩欧美精品v在线| 91在线精品国自产拍蜜月| 国产色爽女视频免费观看| 久久国产乱子免费精品| 国产精品久久久久久精品电影小说 | 亚州av有码| 麻豆精品久久久久久蜜桃| 美女脱内裤让男人舔精品视频| 一级二级三级毛片免费看| 日本三级黄在线观看| 亚洲激情五月婷婷啪啪| 禁无遮挡网站| 成人亚洲欧美一区二区av| 听说在线观看完整版免费高清| 观看免费一级毛片| 永久网站在线| 亚洲欧美日韩卡通动漫| 日日摸夜夜添夜夜添av毛片| 九九久久精品国产亚洲av麻豆| 校园人妻丝袜中文字幕| 欧美97在线视频| eeuss影院久久| 国产精品美女特级片免费视频播放器| h日本视频在线播放| 精品一区二区免费观看| 亚洲精品中文字幕在线视频 | 天天躁夜夜躁狠狠久久av| 欧美zozozo另类| 小蜜桃在线观看免费完整版高清| av又黄又爽大尺度在线免费看| 成人亚洲精品av一区二区| 男人舔奶头视频| 高清欧美精品videossex| 亚洲性久久影院| 国产一区亚洲一区在线观看| 久久精品国产鲁丝片午夜精品| 国产午夜精品久久久久久一区二区三区| 在现免费观看毛片| 国产黄片美女视频| 在线观看美女被高潮喷水网站| 精品久久久噜噜| 国产黄频视频在线观看| 国产成人aa在线观看| 精品国产露脸久久av麻豆 | 一级黄片播放器| 一个人看的www免费观看视频| 一级a做视频免费观看| 啦啦啦啦在线视频资源| 欧美激情国产日韩精品一区| ponron亚洲| 免费电影在线观看免费观看| 性插视频无遮挡在线免费观看| 国产色爽女视频免费观看| 女人久久www免费人成看片| 最新中文字幕久久久久| 亚洲精品久久午夜乱码| 女的被弄到高潮叫床怎么办| 久久精品国产亚洲av涩爱| 亚洲精品日韩av片在线观看| or卡值多少钱| 三级国产精品片| 午夜视频国产福利| 日韩一本色道免费dvd| 亚洲精品乱久久久久久| 欧美一区二区亚洲| 人妻一区二区av| 视频中文字幕在线观看| 白带黄色成豆腐渣| 亚洲精品国产av成人精品| 麻豆久久精品国产亚洲av| av女优亚洲男人天堂| 男女啪啪激烈高潮av片| 少妇人妻精品综合一区二区| 91av网一区二区| 九九爱精品视频在线观看| 国产不卡一卡二| 尤物成人国产欧美一区二区三区| 亚洲自偷自拍三级| 免费看不卡的av| 日本av手机在线免费观看| 免费大片黄手机在线观看| 亚洲婷婷狠狠爱综合网| 永久网站在线| 麻豆乱淫一区二区| 最新中文字幕久久久久| 99热这里只有是精品50| 丰满乱子伦码专区| 人妻制服诱惑在线中文字幕| 一本久久精品| 亚洲激情五月婷婷啪啪| 如何舔出高潮| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 老师上课跳d突然被开到最大视频| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 中文资源天堂在线| 少妇猛男粗大的猛烈进出视频 | 最近最新中文字幕免费大全7| 免费观看精品视频网站| 久久综合国产亚洲精品| 免费av不卡在线播放| 在现免费观看毛片| 人人妻人人澡欧美一区二区| 国产免费视频播放在线视频 | 春色校园在线视频观看| 男女啪啪激烈高潮av片| 久热久热在线精品观看| 色综合亚洲欧美另类图片| 夜夜爽夜夜爽视频| or卡值多少钱| 黄色一级大片看看| 久久国内精品自在自线图片| 日产精品乱码卡一卡2卡三| 18禁动态无遮挡网站| 亚洲av.av天堂| 欧美97在线视频| 青春草亚洲视频在线观看| 高清视频免费观看一区二区 | 1000部很黄的大片| 干丝袜人妻中文字幕| 亚洲美女视频黄频| 成人美女网站在线观看视频| 精品一区在线观看国产| 伦理电影大哥的女人| 一本久久精品| 久久久午夜欧美精品| 国产91av在线免费观看| 国产欧美另类精品又又久久亚洲欧美| 99久久精品一区二区三区| 国产黄色小视频在线观看| 少妇熟女aⅴ在线视频| 亚洲va在线va天堂va国产| 国产乱人视频| 国产v大片淫在线免费观看| 国产成人免费观看mmmm| 久久久精品免费免费高清| 久久精品久久精品一区二区三区| 免费观看在线日韩| 色视频www国产| 菩萨蛮人人尽说江南好唐韦庄| 天堂网av新在线| 中国国产av一级| 国产精品1区2区在线观看.| 男人爽女人下面视频在线观看| 免费高清在线观看视频在线观看| 亚洲久久久久久中文字幕| 亚洲av福利一区| 国产伦一二天堂av在线观看| 精品久久久精品久久久| 午夜福利在线在线| 久久久久久久久久久丰满| 有码 亚洲区| 少妇熟女aⅴ在线视频| 久久国内精品自在自线图片| 久久99蜜桃精品久久| 精品久久久久久久久久久久久| 夫妻午夜视频| 欧美激情久久久久久爽电影| 亚洲精品日本国产第一区| 国产精品无大码| 乱码一卡2卡4卡精品| or卡值多少钱| 最近中文字幕2019免费版| 午夜爱爱视频在线播放| 在线观看一区二区三区| 成人欧美大片| 丝袜美腿在线中文| 国产免费一级a男人的天堂| 麻豆av噜噜一区二区三区| 七月丁香在线播放| 欧美高清成人免费视频www| 色播亚洲综合网| av在线蜜桃| 国产v大片淫在线免费观看| 国产一区亚洲一区在线观看| 99re6热这里在线精品视频| 亚洲欧美清纯卡通| 黄片无遮挡物在线观看| 日本一本二区三区精品| 一个人免费在线观看电影| 日本-黄色视频高清免费观看| 美女高潮的动态| 美女黄网站色视频| 亚洲国产日韩欧美精品在线观看| 亚洲国产最新在线播放| 久久99精品国语久久久| av线在线观看网站| 淫秽高清视频在线观看| 精品不卡国产一区二区三区| 国产精品一二三区在线看| 99久久中文字幕三级久久日本| 国产伦精品一区二区三区四那| 免费观看性生交大片5| 欧美xxⅹ黑人| 欧美日韩亚洲高清精品| 有码 亚洲区| 国产白丝娇喘喷水9色精品| 亚洲欧美一区二区三区国产| 免费播放大片免费观看视频在线观看| 秋霞伦理黄片| 国产精品久久久久久av不卡| 久久热精品热| 亚洲国产精品国产精品| 日韩强制内射视频| 高清午夜精品一区二区三区| 男人狂女人下面高潮的视频| 久久久久久久久久黄片| 欧美人与善性xxx| 国产精品嫩草影院av在线观看| 亚洲精品成人久久久久久| 啦啦啦韩国在线观看视频| 少妇猛男粗大的猛烈进出视频 | av一本久久久久| 欧美最新免费一区二区三区| 秋霞在线观看毛片| 在线免费观看的www视频| 国产探花极品一区二区| 成人av在线播放网站| 亚洲美女视频黄频| 亚洲国产日韩欧美精品在线观看| 高清欧美精品videossex| 免费观看a级毛片全部| 夜夜爽夜夜爽视频| 国产av国产精品国产| 两个人的视频大全免费| 久久久久久久大尺度免费视频| 久久久久久国产a免费观看| 国产男女超爽视频在线观看| 51国产日韩欧美| av国产免费在线观看| 国产激情偷乱视频一区二区| 美女内射精品一级片tv| 2022亚洲国产成人精品| 又黄又爽又刺激的免费视频.| 最近视频中文字幕2019在线8| 亚洲成人av在线免费| 免费看av在线观看网站| 日本猛色少妇xxxxx猛交久久| 黄色日韩在线| 天堂俺去俺来也www色官网 | 国产在视频线在精品| a级毛片免费高清观看在线播放| 国产av码专区亚洲av| 麻豆精品久久久久久蜜桃| 日韩精品有码人妻一区| 国产精品麻豆人妻色哟哟久久 | 美女主播在线视频| 中文字幕亚洲精品专区| 深夜a级毛片| 亚洲精品一区蜜桃| 国产精品国产三级国产av玫瑰| 成人一区二区视频在线观看| 国产男人的电影天堂91| 日日摸夜夜添夜夜添av毛片| 91久久精品国产一区二区成人| 一级毛片黄色毛片免费观看视频| 97精品久久久久久久久久精品| 国产精品一区二区性色av| 99热这里只有是精品在线观看| 热99在线观看视频| 国产探花在线观看一区二区| 国国产精品蜜臀av免费| 97热精品久久久久久| 熟妇人妻不卡中文字幕| 七月丁香在线播放| 国产又色又爽无遮挡免| 女人十人毛片免费观看3o分钟| 在线a可以看的网站| 一区二区三区乱码不卡18| ponron亚洲|