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

    Effectsof Surface Flux Parameterization on the Numerically Simulated Intensity and Structureof Typhoon Morakot(2009)

    2016-08-12 03:41:33JieMINGandJunZHANG
    Advances in Atmospheric Sciences 2016年1期

    JieMINGand Jun A.ZHANG

    1Key Laboratory ofMesoscale SevereWeather/MOE and SchoolofAtmospheric Sciences,Nanjing University,Nanjing 210023

    2Hurricane Research Division,Atlantic Oceanographic and Meteorological Laboratory, NationalOceanographic and Atmospheric Administration,Miami,FL,USA

    3Rosenstiel SchoolofMarine and Atmospheric Science,University ofMiami,Miami,FL,USA

    Effectsof Surface Flux Parameterization on the Numerically Simulated Intensity and Structureof Typhoon Morakot(2009)

    JieMING?1and Jun A.ZHANG2,3

    1Key Laboratory ofMesoscale SevereWeather/MOE and SchoolofAtmospheric Sciences,Nanjing University,Nanjing 210023

    2Hurricane Research Division,Atlantic Oceanographic and Meteorological Laboratory, NationalOceanographic and Atmospheric Administration,Miami,FL,USA

    3Rosenstiel SchoolofMarine and Atmospheric Science,University ofMiami,Miami,FL,USA

    The effects of surface flux parameterizationson tropical cyclone(TC)intensity and structure are investigated using the Advanced Research Weather Research and Forecasting(WRF-ARW)modeling system w ith high-resolution simulations of Typhoon Morakot(2009).Numerical experiments are designed to simulate Typhoon Morakot(2009)w ith different formulations of surface exchange coefficients for enthalpy(CK)and momentum(CD)transfers,including those from recent observationalstudiesbased on in situ aircraftdata collected in Atlantic hurricanes.The resultsshow that the simulated intensity and structure are sensitive to CKand CD,but the simulated track isnot.Consistentw ith previous studies,the simulated storm intensity is found to bemore sensitive to the ratio of CK/CDthan to CKor CDalone.The pressure–w ind relationship isalso found to be influenced by the exchange coefficients,consistentw ith recentnumericalstudies.This paperemphasizes the importanceof CDand CKon TC structuresimulations.The resultssuggest that CDand CKhavea large impacton surface w ind and flux distributions,boundary layer heights,thewarm core,and precipitation.Compared to available observations, theexperimentw ith observed CDand CKgenerally simulated better intensity and structure than theotherexperiments,especially over the ocean.The reasons for the structuraldifferencesamong the experimentsw ith different CDand CKsetups are discussed in the contextof TC dynam icsand thermodynam ics.

    Typhoon Morakot,surface flux parameterization,exchange coefficients,boundary layer

    1. Introduction

    Tropical cyclones(TCs)obtain heat and moisture from the ocean and transfermomentum back to the ocean at the air–sea interface.Thus,surface fluxesof sensible and latent heatplay a very important role in the developmentandmaintenance of TCs.Malkus and Riehl(1960)found that strong surface fluxes can increase the equivalent potential temperature in the TC eyewall,which is linked to the decrease in m inimum sea levelpressure at the TC center(i.e.,increase in TC intensity).Early theoretical studies found that the TC intensity is sensitive to the selection of exchange coefficients for air–sea momentum and heat transfers(Ooyama,1969; Rosenthal,1971;Emanuel,1986).Rotunno and Emanuel (1987)showed that the increase in the exchange coefficient for enthalpy transfer(CK)and the decrease in the drag coefficient(CD)caused the increase in TC intensity in idealized numericalmodels.Based on comparisons ofmodel predictionsw ith observations fora number of hurricanes,Emanuel (1995)concluded that the TC intensity ismost sensitive to the ratio of CKover CD,and this ratiomust lie in the range 0.75–1.5 to achieve consistentmodelsimulationsw ith observations.Full-physical numericalmodel simulations of TCs also suggest that the TC intensity is sensitive to CDand CK(e.g.,Braun and Tao,2000;Nolan etal.,2009).

    Early observationalstudies(e.g.,Smith,1980;Largeand Pond,1981;Geernaertetal.,1988)investigated CDin weak tomoderatew ind conditionsw ith 10m w ind speed(U10)<25 m s?1in theopen ocean.They found that CDincreasesalmost linearly w ith U10.The relationship between CDand U10has the form 1000CD=0.063U10+0.61,whereboth CDand U10have been adjusted to neutral stability(Sm ith,1980).However,results from the Coupled Boundary Layer and Air–Sea Transfer(CBLAST)experiment indicate that CDdoesnot increasew ith U10with no limitation(Black etal.,2007;French et al.,2007).Donelan et al.(2004)estimated CDin awave tank and found that CDreaches amaximum for U10at~33 m s?1and then levels off at higher w ind speeds.They attributed this leveling-off of CDto the increasing shelteringof w inds behind steeperwaves at higherw ind speeds.The above-mentioned studies support the result of Powell et al. (2003),who estimated CDusing the profi lemethod by fi tting hundreds of GPS dropsonde w ind profi les collected in hurricanes to a logarithmic relationshipw ith height.Ming et al.(2012)showed that CDestimated based on observational data in typhoonsalso behavessim ilarly as in the Powelletal. (2003)study.

    Due to instrumentation limitations associated w ith sea saltand/orseaspray,directmeasurementsofsurfaceenthalpy fluxeshavebeen very difficultover theocean.Theonly available in situ observations over the past30 years(DeCosmo et al.,1996;Drennan etal.,2007)show that the exchange coefficient for latentheat transfer(CE)doesnotdepend on the w ind speed.Sim ilarly,the bulk sensible heat exchange coefficient(CH),or Stanton number,has also been found to be independent of wind speed.Note that CKis generally assumed to be equal to CEand CH,which is supported by the CBLAST result(Zhang et al.,2008).Wave tank observations also support the result from the CBLAST experiment, that CKis nearly constantw ith U10up to~40m s?1(Haus etal.,2010).Belletal.(2012)also found thatvaluesof CKestimated based on budget analyses of dropsonde and radar data in theboundary layerof theeyewall regionwere close to those reported by Zhang etal.(2008).

    Most previous numerical studies on the impact of surface exchange coefficients on TC intensity and structure focused on Atlantic hurricanes or used idealized simulations (e.g.,Montgomery etal.,2010;Bryan,2012;Baoetal.,2012; Green and Zhang,2013).How sensitive the simulated intensity and structure inWestPacific typhoonsare to CDand CKis notwell documented.The objective of this study is to investigate theeffectsof CDand CKon the simulated intensity and structureof Typhoon Morakot(2009).Emphasisisplaced on simulations of structure,such as boundary layer heights,the warm core,and precipitation.Typhoon Morakot(2009)possessed unusual structures,including a large eye and strong outer circulation;thus,this storm was atypical compared to most hurricanes in the Atlantic basin.Thismotivates us to evaluate the choice of CDand CKin the numerical simulations to reproduce theunusualstructuresof TyphoonMorakot (2009).

    2. Typhoon Morakot(2009)

    Typhoon Morakot(2009)generated from a tropical depression over the Philippine Sea on 2 August2009.Itgradually intensified to a tropical storm and was assigned the name Morakot on 3 August.A fter that,itmoved westward towards Taiwan.Morakot(2009)intensified to typhoon intensity at1800 UTC 5 August,w ith the location of its center at(23.1°N and 129.3°E)and amaximum w ind speed of 33.4 m s?1.Morakot(2009)continued to move westward and made landfall in central Taiwan at 1200 UTC 7 August.One day later,Morakot(2009)turned tomove northwestward and went back over water into the Taiwan Strait. Subsequently,Morakot(2009)weakened to a severe tropical storm andmade itssecond landfallin Xiapu,Fujian Province, on the east coast of China at0820 UTC 9 August(Fig.1a). It gradually weakened as itmoved northward.Roughly 24 hours later,itweakened to adepression during the landfall in Zhejiang Province,China.

    Fig.1.The(a)tracksand time seriesofsimulated(b)maximum surfacew ind speed and(c)minimum central pressure(hPa)of Morakot(2009)from 0000UTC 5August to 0000UTC 10August 2009 from the three experiments and the corresponding best-track analysisby JMA.

    TyphoonMorakot(2009)broughtrecord-breaking torrential rainfallover Taiwan(Chien and Kuo,2011).This catastrophic storm sadly claimedmore than 600 livesand causedmore than 200 people to be classed asm issing,creating a total cost of damage estimated at US$3.3 billion(Zhang et al.,2010).It also affected more than 11 m illion people throughouteastern China and damaged thousandsof homes. Asa result,Typhoon Morakot(2009)has received attention in severalnumericalstudies.

    Schwartz et al.(2012)studied Typhoon Morakot(2009) using the Advanced Research Weather Research and Forecasting(WRF-ARW)modeling system w ith dataassimilation of them icrowave radiances using a cyclic,limited-area ensemble adjustment Kalman fi lter.They found that the track, intensity,and precipitation forecastswere improved after assim ilatingmicrowave radiances.Xie and Zhang(2012)performed ensemble simulations using WRF-ARW to investigate the dynam ics and predictability of the record-breaking rainfall and flooding event in Taiwan induced by Typhoon Morakot(2009).They found that the large typhoon circulation and the southwesterly monsoon flow transported abundantmoisture into southern Taiwan,which,along w ith the influence of the complex high terrain,produced the heavy rainfall.

    Wang etal.(2012)used theCloud-Resolving Storm Simulator(CReSS)to study the dynam ics related to themotion of Typhoon Morakot(2009)over Taiwan.Their simulations showed that the reduced moisture content induced the decrease in the rainfalland the increase in thestorm translation speed.They also pointed out that the asymmetric precipitation in Typhoon Morakot(2009)played an important role in itsvery slow motion upon leaving Taiwan,the lengthening of theheavy-rainfallperiod,and the increaseof the total rainfall amount.Furthermore,their results emphasized the potential contribution of asymmetric heating to the slowdown of the typhoon motion in the presence of complex terrain or in a monsoon environment.Wang etal.(2013)demonstrated that CReSS successfully simulated and reproduced both the distribution and tim ing of theheavy rainfall in Taiwanw ith high accuracywhen TyphoonMorakot(2009)passed by.The realtime forecast integrations of the CReSSmodelalso showed high-quality quantitative precipitation forecasts.

    Hall et al.(2013)used the Advanced Regional Prediction System to simulate Typhoon Morakot(2009)when it made landfallover Taiwan.They investigated themesoscale structure ofMorakot(2009)and emphasized the roleof deep convection on the rainfall simulation.They identified that relatively large-amplitude wave structures developed in the outereyewall,known asvortex Rossbywaves(VRWs).They found that the strong asymmetry of the convection was associated w ith wavenumber1(WN1)VRWs,while theWN2 and WN3 VRWs were associated w ith the development of the deep convective band in Morakot’s(2009)southwestern quadrant.Huang et al.(2014)used the WRFmodel to explicitly simulate Typhoon Morakot(2009)and found that the simulated rain rateand precipitation efficiency(PE)over the CentralMountain Range(CMR)were highly correlated. They found also that the PE and the processesof vapor condensation and raindrop evaporation were strongly influenced by orographic lift.They also found that the increase in PE over the CMR comparedw ith thatover theoceanwasdue to an increase in the ice-phase deposition ratiowhen the liquidphase condensation reduced as theairon the leesidesubsided andmoved downstream.They emphasized the effects of the terrain on the simulationsof precipitation.

    3. Experimentaldesign

    Asmentioned above,the objective of this paper is to investigate the sensitivity of the simulated track,intensity and structureof TyphoonMorakot(2009)to thesurfaceexchange coefficients.We used the ARW modeling system(version 3.2;Skamarock et al.,2008)to conduct numerical experiments.Three experiments were performed using the ARW modelw ith twomodel domains.Themodelhorizontalgrid resolution was 4.5 km for the parent domain(D1)and 1.5 km for the nested domain(D2).The two domains covered 600×400 and 421×421 grid points,respectively(Fig.1a). D2 was an automatic vortex-follow ingmoving nestand the center of the domain was always located at the center of the storm.In all theexperiments,35 vertical(σ)levelswereused from the surface to themodel top at50 hPa.

    D1 was initialized at0000 UTC 4 August2009 and was run for 6 days.D2 was started at 0000 UTC 5 August after a one-day spin-up.The outer domain was run in a fourdimensional data assim ilation(FDDA)mode to provide the bestpossible large-scaleconditions for the innerdomain.The innernested domainwas runwithout FDDA.The Japan Meteorological Agency(JMA)6-hourly gridded regionalanalysesat20×20 km horizontal resolution w ith 20 pressure levelswere used as the initial and boundary conditions for the ARW model.The JMA analyseswereproduced using amultivariate three-dimensionaloptimum interpolationmethod to combine the fi rst-guess fields from JMA’s regional spectral model(RSM)with observations from a variety of platforms (JMA,2007;Hosomi,2005).Note that the JMA analysesdid notcapture the realstructureand intensity ofMorakot(2009). Compared to the best track,the error of them inimum sealevelpressurewas13.7hPa,and thatof themaximum surface w ind speedwas10m s?1.The track of the storm in the JMA analyses also has a northwestward bias.Follow ing Ming et al.(2009),we appended a good vortex w ith the appropriate intensity(close to the best track)into the initial conditions. We ran D01 only from 0000UTC to 1600UTC 4 August,to make the storm intensity sim ilar to the best track.The threedimensional vortex was extracted and inserted back into the initial condition using the position of the best track,which corrected the biasof the storm centerat the initial time.Under the new initial condition,the error of them inimum sea level pressurewas 0.6 hPa and thatof themaximum surface w ind speedwas5.6m s?1.

    The physics schemes used in the numericalexperiments included the Purdue Lin m icrophysics scheme(Lin et al., 1983;Chen and Sun,2002),the Yonsei University boundary layer scheme(Noh etal.,2003;Hong et al.,2006),the Rapid Radiative TransferModel longwave radiation scheme (Mlawer et al.,1997),and the Dudhia shortwave radiationscheme(Dudhia,1989).Cumulus parameterization schemes werenotused forboth domains.

    Table 1.Thedifferentroughness lengths(z0)and thermal roughness lengths(z0q)in the three experiments.

    In this study,three experiments(see Table 1)were designed to exam ine the impactof CDand CKon the typhoon intensity and structuresimulations.Asmentioned earlier,surface fluxes were parameterized through CDand CK.In the surface-layer schemes of WRF-ARW,CDand CKwere parameterized using surfacemomentum(z0)and thermal(z0q) roughness lengths.The drag coefficient in the neutral conditionwas defined as:

    where k is the von K′arm′an constant.The enthalpy exchange coefficientwasdefined as:

    In the control experiment(referred to as CTL hereafter),z0depended on thew ind speed using the Charnock(1955)relationship,in the form of:

    where u?is the friction velocity,g is the gravitationalacceleration.This formula causes z0to increase continuously w ith the w ind speed,w ith no limit.For z0qthe Carlson–Boland scheme(Carlson and Boland,1978)wasused:

    where xkais a constant that equals 2.4×10?5.Since z0qvariesmore slow ly w ith the friction velocity than z0,CKincreasesmore slow lyw ith thew ind speed than CD(Dudhiaet al.,2008).

    In the second experiment(referred to as TC1 hereafter), an alternate CDformulation based on the high-wind laboratory experimentof Donelan etal.(2004)was adopted.The roughness lengthwas defined as:

    where 1.27×10?7z02.85×10?3(units:m).Values of CDbased on Eqs.(1)and(5)are smaller than those from the Charnok relationship.CDincreases almost linearly with the 10 m w ind speed,up to a maximum of 0.0024 at 35 m s?1(Davis et al.,2008).In addition,the thermal roughness lengthwascalculated using the ramped formula(Dudhia etal.,2008)as follows:

    Thissetup of z0and z0qcauses CKto increasealmost linearly w ith w ind speed,and CKbecomes larger than CDatw ind speedsof>50m s?1.Note that the ratio CK/CDexceeds1 as thew ind speed becomes larger than 50m s?1,although the condition for themaximum surfacew ind speed of>50m s?1wasnevermet in oursimulations(Fig.1b).

    In the third experiment(referred to as TC2 hereafter),the drag formulationwasstillbased on the resultofDonelan etal. (2004),but the thermal roughness length was adopted based on theGarratt formula(Garratt,1992,p.102)as follows:

    where T is air temperature at the surface,γis the kinematic viscosity of air(units:m2s?1).In this setup,CKis related to u?and T.CKdecreasesw ith increasingw ind speed in smooth flow,while it isa constant in rough flow orhigh-wind conditions(Garratt,1992,p.102).Note that CKusing the Garratt formula is closer to recent observations from CBLAST and thewave tank experimentmentioned in the introduction.

    All simulations in the three experimentswere run from 0000 UTC 5 August to 0000 UTC 10 August 2009,which covered almost the entire lifecycle of Typhoon Morakot (2009),from the intensifying stage to the landfalling stage. A ll experiments were conducted w ith the same initial vortex and boundary conditions,so the differences in the experiments were solely related to the different formulas of the exchange coefficientsused in theWRF-ARW model.

    The exchange coefficients and ratio of CK/CDas a function of the surfacewind speed are shown in Fig.2.In CTL, the Charnok relationship wasused such that CDincreased almost linearly w ith the 10 m w ind speed,because of the increasing surface roughness.CKalso increased slow ly w ith the w ind speed.In TC1 and TC2,when the Charnok relationship was changed to the new formulation of CDbased on Donelan etal.(2004),CDincreasedw ith thew ind speed up to~33m s?1,then levelled off.In CTL,CKwas based on the Carlson–Boland formula and was larger than that based on the ramped formula used in TC1 forw ind speeds of<30m s?1.On the other hand,the CKquantities produced by these two formulaswere close to one another for w ind speeds of>30m s?1.The Garratt formula for the thermal roughness lengthwasused in TC2,where the CKwas smaller than thatin CTL and TC1 forw ind speedsof>20m s?1.In TC2,the ratioof CK/CDwas larger than thatin CTL forw ind speedsof<20m s?1,andmarginally smaller forwind speeds of>20 m s?1.Note that the ratio of CK/CDin TC1 was larger than that in the other two experiments forallw ind speeds.

    Fig.2.The simulated(a)drag coefficient(CD),(b)enthalpy exchange coefficient(CK)and(c)exchange coefficient ratio (CK/CD)as a function of 10 m w ind speed in the whole of domain 2 from the three experiments at 1800 UTC 6 August 2009.

    4. Results

    4.1. Track and intensity

    Thesimulated track of Typhoon Morakot(2009)from the threeexperiments from 0000UTC 5August to 0000UTC 10 Augustare comparedw ith the JMA best track in Fig.1a.As mentioned earlier,TyphoonMorakot(2009)movedwestward in the fi rst two days,and turned northwestward aftermaking landfall in Taiwan.It then continued tomovenorthward,and made a second landfall in Fujian Province.The simulated storm tracks in the three experimentswere quite sim ilar,although they deviated from the best track.A ll three experiments reproduced the observed west-northwestward storm motion.Overall,the simulated tracksare notsensitive to the different formulasof exchange coefficients.

    The simulatedmaximum surfacew ind andm inimum sea levelpressure from the threeexperimentsare comparedw ith the best track from JMA in Figs.1b and 1c.Although all threeexperimentswere initializedw ith thesame initialcondition,the intensity forecastsweredifferent,especially after the maximum w ind speed reached 30m s?1at the forecast time of 0000 UTC 6 August.In the next two days,all three experiments simulated maximum surface w ind speed thatwas larger than observed.Themaximum surfacew ind speed simulated in CTL and TC2 were closer to the best track.The difference in the simulated m inimum sea level pressure,on the other hand,was clearer than that in themaximum surfacew ind speed among the threeexperiments.The simulated m inimum sea level pressure in TC2 was closest to the best track among the threeexperiments.

    Figure 3 plots the relationship between the maximum w ind speed and m inimum sea level pressure from the three experiments and the best track of JMA.It appears that the pressure–w ind relationship is sensitive to CDand CK.The result is consistentw ith a recentnumerical study carried out by Green and Zhang(2013),who found that the pressure–w ind relationship was sensitive to the selection of different surface layer schemes inWRF simulationsof Hurricane Katrina(2005).Our result is also consistent w ith Bao et al. (2012),who ran different surface-layer physics in idealized simulationsof theHurricaneWeatherand Research Forecasting(HWRF)modeland confi rmed that thepressure–w ind relationship was sensitive to CDand CK.

    The simulated storm intensity and pressure–w ind relationship in TC2 were closest to the best track among the threeexperiments,especially forw ind speedsof>30m s?1. This resultisencouraging because the CDand CKused in TC2 were closer to recentobservations than those in theother two experiments.Our resultgenerally indicates thatwhen CKis larger,the storm obtainsmoresensibleand latent fluxes from the underlying ocean,resulting in larger surfacew ind speed and lower central pressure.On the other hand,when CDis larger,surface friction is larger,which tends to reduce the surfacew ind speed.Itwas pointed outby(Montgomery et al.,2010)that larger surface friction can also lead to a larger gradientw ind imbalance in the boundary layer,which couldlower the central pressure.Thisexplainswhy the simulated m inimum sea level pressure in CTL was lower than that in TC2 but themaximum wind speed in CTL was smaller than that in TC2.

    Fig.3.Scatterplotofm inimum sea levelpressure vsmaximum surface w ind speed.The blue diamonds represent CTL,the red circles represent TC1,and the green crosses represent TC2. Theopen squares representthepressure–w ind relationship from JMA.

    4.2. Evolution ofprimaryand secondary circulationsand thewarm core

    Time–radius Hovm¨oller diagrams of azimuthally averaged tangentialw ind velocity(Vt)at the altitude of 2 km and radialwind velocity(Vr)at the altitude of 250m,from 0000 UTC 6August to 1200UTC 7August,aredisplayed in Figs. 4a–c and 4d–f,respectively.As can be seen,the evolution of theaxisymmetric Vtand Vrin the three experimentswassimilar.In the fi rst6 h of the simulations,Morakot(2009)wasa weak storm thathad amaximum Vtof25m s?1andminimum Vrof?10 m s?1,located at the radius of~170 km.Later, Vtincreased gradually w ith time and the radiusofmaximum tangentialw ind speed(RMW)became smaller.Themagnitude of Vrdoubled from?10 to?20m s?1in the next24 h and the radiusof the peak inflow also contracted(Figs.4d–f) in response to the intensification of the storm.Inward from the RMW,the radial velocity decelerated rapidly ata larger rate thanoutside theRMW.Afteranother16 h,thestorm centerwas close to Taiwan such that themagnitudes of both Vtand Vrdecreased.Among the three experiments,the peak Vtand Vrin TC1 were the largest,mainly because the ratio of CK/CDused in TC1was larger than that in CTL and TC2.

    Fig.4.Time–radius Hovm¨oller plots of azimuthally averaged tangentialw ind(units:m s?1)at2 km altitude from(a) CTL,(b)TC1 and(c)TC2,and radialw ind(units:m s?1)at250m altitude from(d)CTL,(e)TC1 and(f)TC2.The thick lines depict the RMW at2 km altitudeand themaximum inflow at0.25 km altitude.

    The time–height evolution of the mean temperature anomalies(referred to as warm-core anomalies hereafter) from 0000 UTC 6 August to 1200 UTC 7 August is plotted in Fig.5 for CTL,TC1 and TC2.Follow ing Liu etal.(1999) and Lietal.(2013),we fi rstcalculated themean temperature w ithin the region of 630×630 km from the typhoon’sm inimum surface pressure centerateach vertical level.Then,the temperature anomaly was obtained by subtracting themean temperature from the temperatureateach grid pointand each level.Thewarm-coreanomaly was then defined as the average value of the temperature anomaliesw ithin the region of 300×300 km from the storm centerateach level.Theheight of the peakmean temperature anomaly represents thewarm core height.It is evident from Fig.5 that the peak warm coreanomaly in CTL and TC2wassmaller than thatin TC1, indicating that the warm-core anomaly was correlated w ith the storm intensity.According to thehydrostatic balance,the lower theminimum sea level pressure the larger thewarmcoreanomaly(Zhang and Chen,2012).Thus,thewarm-core anomaly in TC2 was the smallest among the three experiments as the storm intensity in TC2 was the lowest.On the otherhand,thewarm-coreheightwasnotcorrelatedw ith the storm intensity,because thewarm-coreheightwas located at 8–12 km for the three experiments.

    4.3. Surface wind,surface flux and boundary layer heights

    The horizontal distributions of 10 m w ind speed valid at 1800 UTC 6 August 2009 are shown in Figs.6a–c for CTL,TC1 and TC2,respectively.It isevident that themaximum surfacew ind speed in TC1was larger than that in the other two experiments.The eyewall region,which covered themaximum w ind speed,was broader in TC1 than in CTL and TC2.Comparing TC1 and TC2,the resultsuggests that increasing CKalone would increase the storm intensity in termsof themaximum w ind speed in theeyewall region.The surfacew ind distribution in TC2 was sim ilar to that in CTL, but themaximum w ind speed in the right-rear quadrantwas slightly larger in TC2 than in CTL.Although CKin TC1was close to that in CTL,themaximum surface w ind speed in TC1 was larger than that in CTL,especially forw ind speeds of>25m s?1.Thisdifferencewasmainly due to the fact that the CDused in TC1wassmaller than that in CTL.

    Large values of latent fluxes were found in the eyewall and primary rain band regions where surface w ind speeds were also large(Figs.6d–f).Themaximum latentheat flux in TC2 wasmuch smaller than that in the other two experiments,mainly because the CKused in TC2 was smaller than that in CTL and TC1(Fig.2).The difference in the sensibleheat fluxesamong the threeexperiments(Figs.6g–i)was much smaller than the difference in the latent heat fluxes, although the same formulawasused for calculating both the latentand sensible heat fluxes.Itappears that the simulated maximum sensibleheatflux in TC1wasthe largestamong the three experiments,while that in TC2 was the smallest.The asymmetric distribution of the sensible heat fl ux in TC2 was closer to that in CTL than in TC1.Overall,the result(Fig. 6)implies that CKand CDalonehaveoppositeeffectson surface enthalpy flux,and CKinfluences the enthalpy fluxmore than CD.A larger CKinducesmoresensibleand latent fluxes, which supportmore energy for a storm to intensify.On the otherhand,a larger CDinduces largersurface friction,which reduces the surfacew ind speed and in turn reduces the sensi-ble and latent fluxes because these fluxes are also a function of thew ind speed.

    Fig.5.Time–height diagrams of temperature deviation(units: K)from the three experiments:(a)CTL;(b)TC1;(c)TC2. The average was computed w ithin the area of 300×300 km from the surfaceminimum pressure center for simulations.The anomalieswere obtained by subtracting the averaged temperaturew ithin the region atevery height level.

    Fig.6.Themodel-simulated 10m w ind speed(colorscale;units:m s?1)andwind vectors(arrows;units:m s?1)from (a)CTL,(b)TC1 and(c)TC2;latentheat flux(color scale;units:W m?2)and w ind vectors(arrows;units:m s?1) from(d)CTL,(e)TC1 and(f)TC2;and sensibleheat flux(colorscale;units:W m?2)and w ind vectors(arrows;units: m s?1)from(g)CTL,(h)TC1 and(i)TC2 at1800UTC 6 August2009.The large vector indicates themotion of the storm,and the thin crosses divide the storm into fourquadrants.

    Figure7 shows the radius–heightplotsof theazimuthally averaged tangentialand radialvelocities forall three experiments at 1800 UTC 6 August 2009.Themagnitudes of the tangentialand radial velocities in TC1 were generally larger than those in the other two experiments,consistentw ith the simulated storm intensities.However,the peak tangential w ind speed in TC2was smaller than that in CTL,which was notconsistentw ith the intensity differencebetween these two experiments in terms of themaximum surfacew ind,as TC2 had a largermaximum surfacewind speed.Thisdiscrepancy can be explained by the role of CDin regulating the boundary layer dynam ics.Montgomery et al.(2010)pointed out thatan increase in CDleads to an increase in storm intensity in terms ofmaximum tangentialw ind speed in the boundary layer,although it reduces the surfacew ind speed through surface friction.Our result is consistentw ith that of Montgomery etal.(2010),indicating that surface flux parameterization affects the vertical structure ofw ind velocitiesabove the surface layer.

    In many PBL schemes used in full-physics numerical models,one of the crucial elements is the boundary layer height,because it is coupled w ith the energy transport from the surface layer to the boundary layer and above(e.g.,Beljaars and Viterbo,1998;Noh et al.,2003).The boundary layer height is also a key variable that regulates the verticaldistribution of turbulent fluxesand helpsdeterm inewhere turbulent fluxes tend to become negligible(Stull,1988).Follow ing Zhang et al.(2011),the kinematic boundary layer height is defined by the heightofmaximum tangentialw ind speed(hvtm).Inflow layer depth(hinf),defined as the height where the inflow reduces to 10%of the peak value,also represents the kinematic boundary layer height.In all three experiments,hvtmdecreased with decreasing radius toward the storm center(Figs.7a–c).Thisbehavior isconsistentw ith the resultgiven by Zhang etal.(2011),who composited hundreds of dropsondes data collected from 13 hurricanes to study thecharacteristics of hurricane boundary layer heights.Within the radiusof125 km,hvtmin TC2wassmaller than thatin the other two experiments,andwascloser to observations(Zhang etal.,2011,Fig.5a).

    Fig.7.Azimuthally averaged radius–heightcrosssectionsof tangentialw ind(units:m s?1)from(a)CTL,(b)TC1 and (c)TC2;radialw ind(units:m s?1)from(d)CTL,(e)TC1 and(f)TC2;and the gradient force imbalance(FPEC,units: m s?1h?1)from(g)CTL,(h)TC1 and(i)TC2 at1800UTC 6 August2009.The thick lines in(a–c)depict the height of themaximum w ind speed varying w ith radius,and the thick dashed lines in(d–f)depict the inflow layer height, defined as the heightwhere the radialw ind speed is10%of the peak inflow.

    The difference in the inflow layer depth(hinf)among the three experiments wasmuch larger than that in hvtm(Figs. 7d–f).Itappears that hinfin TC1 was the highest.A ll three experiments captured the decrease of hinfw ith deceasing radius,consistent w ith observations.Furthermore,hvtmwas smaller than hinfin all three experiments,which was also consistentw ith observations.Overall,themagnitude of hinfin TC2 was closest to observations(Zhang etal.,2011,Fig. 5b).The above results indicate that boundary layer heights are tied to the surface flux parameterization.

    According to Zhang et al.(2001),themomentum equation of the radialw ind velocity in the cylindrical coordinates system can bew ritten as

    Fig.8.(a)U.S.Air Force DefenseMeteorological Satellite Program(DMSP)polar-orbiting satellite F-17microwave imagery of polarization corrected temperaturew ith horizontal polarization w ith 91 GHz at2112 UTC 6 August,and simulated composite radar reflectivity(units:dB Z)at2100UTC 6 August2009 from(b)CTL,(c)TC1 and(d)TC2. (e)The observed composite reflectivity and simulated composite reflectivity from experiment(f)CTL,(g)TC1 and(h) TC2 at1200 UTC 7 August2009.(i)The observed composite reflectivity and simulated composite reflectivity from experiment(j)CTL,(k)TC1 and(l)TC2 at0000UTC 8August2009.

    where w is verticalw ind velocity;Ωis the angular velocity andφis the latitude;r is the radius from the center;λis the azimuthal angle.Equation 11 states that the radial acceleration isdetermined by the radial pressuregradient force (FP;fi rst term on the right-hand side of the equation),the centrifugal force(FE;second term),the Coriolis force(FC; third and fourth terms),and diffusion(Ud;last term).Thedegreeofgradient force imbalanceornetagradient force(FPEC) is evaluated by adding FP,FE,and FCtogether(Figs.7g–i). Firstly,FPECwas larger in TC1 than in the other two exper-iments,supporting the fact that the simulated storm in TC1 wasstronger,because thestorm tended to spin-up fasterwhen FPECwas larger(Sm ith et al.,2009).A lthough the CK/CDwas alike in CTL and TC2,FPECwas larger in CTL than that in TC2.Thiswasmainly due to the different CDused in those two experiments.Follow ing the dynamicalexplanation of Montgomery et al.(2010),the agradient tendencies near thesurface caused the inflow ing ringsofboundary-layer air to converge farther inwards in the storm centerbefore rising outof the boundary layerand ascending into theeyewall updraught,resulting in enhanced maximum tangentialw ind speed.Our result is consistentw ith this argument(Fig.7). Italso suggests that intense positive supergradientacceleration occurs in the vicinity of themaximum tangentialw ind speed and isassociated w ith the outflow jetabove the boundary layer.

    4.4. Radarreflectivityand precipitation

    Next,we investigate the simulated radar reflectivity and precipitation in the three experiments.Figure 8a shows the observed U.S.Air Force Defense Meteorological Satellite Program polar-orbiting satellite F-17microwave imagery of polarization corrected temperature w ith a horizontal polarization at 91 GHz valid at 2112 UTC 6 August.The simulated composite radar reflectivity valid at 2100 UTC 6 August from CTL,TC1 and TC2 are shown in Figs.8b–d,respectively.Due to the interaction of the typhoon circulation w ith themonsoon flow and verticalw ind shear(Wang etal., 2012),the storm became asymmetric.All three experiments captured theasymmetric distribution of radar reflectivity and reproduced the unclosed eyewall in the southern part of the storm.The simulated reflectivity in TC2 was only slightly closer to observations than in the other two experiments,as it captured broaderhigh reflectivity area in the southern eyewall.Otherw ise,the overall rainfall structurewas sim ilar in allexperiments.

    Theobserved and simulated radar reflectivity composites at1200UTC 7Augustand 0000UTC 8Augustareshown in Figs.8e–h and Figs.8i–l,respectively.At1200 UTC 7 August,Typhoon Morakot(2009)was located on the east side of Taiwan before landfall.It is evident from the observation(Fig.8e)that thestrongest reflectivitywas located on the southwestside of the storm.Thisasymmetric rainfallpattern wascaptured by all theexperiments.The simulated reflectivity in the eyewall region was stronger in CTL and TC1 than in TC2(Figs.8f–h).The simulated reflectivity in TC2 was slightly closer to observations than in the other two experimentsbecause the unclosed eyewall,w ith the principal rainband located in the southern partof the storm and the strong echoover Taiwan,werecaptured in TC2.At0000UTC 8August,the observed precipitation pattern becamemore asymmetric than in earlier periods.The strong echoes were observed on the south side of Taiwan(Fig.8i).In CTL and TC1,the simulated reflectivity wasstillstronger than that in TC2(Figs.8j–l).Again,TC2 performed slightly better than in the other two experiments.

    Overall,the total precipitation was strongest in TC1 and weakest in CTL,especially after 1200 UTC 6 August(Fig. 9a).Note that theaccumulation period of precipitation is1 h. As TC1 simulated the strongest storm while CTL simulated theweakeststorm,this resultsuggests that the totalprecipitation iscorrelatedw ith thestorm intensity.The10m domain–averaged divergencesofmoisture flux(note thatnegativevalues representconvergence of themoisture flux)are shown in Fig.9b.It is evident that the convergence ofmoisture flux was correlated w ith the total precipitation.This result isnot surprising,as the low-levelmoisturewas themain source of the rainfall.Nonetheless,the resultsuggests thatsurface flux parameterizations have a substantial impacton precipitation simulations.

    Although the simulated rainfall over the ocean was strongly tied to surface flux parameterization,interestingly, we found that the rainfall over land(i.e.,Taiwan)wasmuch less sensitive to the surface flux parameterization.Figure 10 compares the 12 h accumulated precipitation from CTL, TC1 and TC2 togetherw ith the objective analysesof rainfall measured by automaticweatherstationsover Taiwan,valid at 0600UTC and 1800 UTC 7 August.Prior to the landfallof Morakot(2009)in Taiwan(1800UTC 6Augustto 0600UTC 7 August),the observation(Fig.10a)shows two regions of strong precipitation across the island:oneon thenorth sideofTaiwan and the other on the south side over highmountains. A ll three experiments simulated these two regions of strong rainfall,although the simulated rainfall wasmuch stronger than the observed valued.In particular,all three simulations over-predicted the precipitation from central to southern Taiwan.However,the difference in the precipitation distribution in the threeexperimentswasvery smallover thewhole island.

    Fig.9.Averaged(a)accumulated rainfall and(b)10 m divergence ofmoisture flux from experiment CTL,TC1 and TC2. The average was computed w ithin the area of 300×300 km from the surfacem inimum pressure center.

    Fig.10.12 h accumulated precipitation(units:mm)valid at 0600 UTC 7 August 2009 from(a)automatic weatherstation hourly observations,(b)CTL,(c)TC1 and(d)TC2 over Taiwan;and 12 h accumulated precipitation valid at1800 UTC 7 August2009 from(e)automaticweather station hourly observations,(f)CTL,(g) TC1 and(h)TC2 over Taiwan.

    The above result suggests that the precipitation was lesssensitive to the exchange coefficients over land than over the ocean in Morakot(2009).Over the ocean,the role of the exchange coefficient in precipitation wasan indirectone, through influencing the(horizontal)moisture flux convergence(and ultimately rainfall)and through affecting the intensity of thestorm.For theprecipitationover land,themoisture came from the ocean,even though the end precipitation fellover land.In Typhoon Morakot(2009),the effectof terrain(i.e.,forced uplift)played adom inantrole in thedistribution of theprecipitation over Taiwan(Halletal.,2013;Wang etal.,2013),which is likely themain reason for the sim ilar rainfallsimulationsamong the threeexperiments.

    5. Summary and discussion

    In this study,three numerical experiments were performed w ith the WRF-ARW model to study the impact of surface flux parameterizations on the structure and intensity of Typhoon Morakot(2009).The initial conditions of the threeexperimentswereall from the JMA RSManalysis field. The simulated track and intensity of Morakot(2009)were verified against thebest track.Different formulasofmomentum and heat roughness lengthswere tested in sensitivity experiments thatgoverned thebehaviorof the surfaceexchange coefficients formomentum and heat transfers.The results showed that the simulated track was notsensitive to the exchange coefficients,but the simulated intensity and structure were.

    Our results indicate that the surfaceexchange coefficients are key factors for the simulation of surfacew ind speed and fluxes.The effect of CKon the surface enthalpy flux is straightforward because of the linear relationship between these two parameters.On the other hand,the effect of CDon the enthalpy flux takes place via the surfacew ind speed. When the CDissmall,themaximum surfacew ind speed tends to be larger due to reduced surface friction.In turn,the enthalpy flux becomes larger because of the largerw ind speed. Overall,we found CKhad a larger impacton theenthalpy flux simulation than CD.

    Consistentw ith previous studies(e.g.,Emanuel,1995), the simulated storm intensity was found to bemore sensitive to the ratio of CK/CDthan to CKor CDalone.According to the idealized numerical simulation given by Montgomery et al.(2010),CK/CDshould havea criticalvalue for the intensification of storm.If CDis too large,the storm w illnot intensify.When CK/CDis larger,the simulated storm is stronger and vice versa.In the CTL experiment,the intensity simulation is comparable to that in TC2,because a sim ilar CK/CDwasused in these two experiments.

    Thepressure–w ind relationshipwasalso found to besensitive to CDand CK,consistentw ith recentnumericalstudies of Atlantic hurricanes(Bao et al.,2012;Green and Zhang, 2013).Overall,the simulated intensity and pressure–w ind relationship in TC2wasclosest to thebest track than those in CTL and TC1.This resultisencouraging because the CKand CDused in TC2 were close to recent field and wave tank observations.This resultisalso consistentw ith thatof Zhang et al.(2012),who showed thatobservation-based surface layer and boundary layer physics led to improvements in the operationalHWRFmodeland better intensity forecasts.

    Our resultsalso indicate thatsimulated structures,such as the surfacew ind distribution,boundary layerheights,warmcore anomaly and height,and precipitation are affected by CDand CK.Compared to the dropsonde observations from Zhang etal.(2011),the simulated kinematic boundary layer heights in TC2 are closer to observations than the other two experiments.Thewarm-coreanomaly is tied to the storm intensity but not the warm-core height,consistentw ith Stern and Nolan(2012).The difference in the rainfall over the ocean is consistent w ith the difference in storm intensity, which can be explained by the difference in the convergence ofmoisture flux in the boundary layer.Over land,the simulated rainfall ismuch less sensitive to CDand CKthan over the ocean,which we attribute to the dominance of the terrain effect on the precipitation in Typhoon Morakot(2009), aspointed outbyWang etal.(2013)and Halletal.(2013).

    Wealso conducted dynamicalanalyses to investigatewhy CDand CKaffect the vertical structure of w ind velocities in the boundary layer.Consistentw ith Montgomery et al. (2010),we found thata larger drag coefficient can lead to a largergradientw ind imbalance in theboundary layer(Fig.7). Asa resultof the largeragradient forcing,theboundary-layer air converged farther inward near the storm centerbefore rising outof the boundary layer and ascending into the eyewall updraft.The end resultwas enhanced maximum tangential w ind speed,despite the loss of absolute angularmomentum en route.

    In this study,we focused on investigating the sensitivity of the simulated intensity and structure of Typhoon Morakot (2009)to the surfaceexchange coefficientsonly,while keeping the restof themodelphysics the same.Wenote thatother partsof themodelphysics(e.g.,planetary boundary layerparameterization and radiation parameterization)may also be important for TC simulations.Futurework w illevaluate the impactof otheraspectsofmodelphysicson numericalsimulationsof TC structure and intensity change.

    Acknow ledgements.Jie MING was primarily supported by the National Fundamental Research 973 Program of China(Grant Nos.2015CB452801 and 2013CB430100),the National Natural Science Foundation of China(Grant No.41105035),and the Fundamental Research Funds for the Central Universities(Grant Nos.20620140054 and 20620140347).Jun ZHANG was supported by NOAA’s Hurricane Forecast and Improvement Project(HFIP), GrantNos.NA14NWS4680028 and NASA GrantNNX14AM69G. Wearegrateful to theHigh PerformanceComputing Centerof Nanjing University for carrying out the numerical calculations in this paperon its IBMBlade cluster system.

    REFERENCES

    Bao,J.-W.,S.G.Gopalakrishnan,S.A.Michelson,F.D.Marks, and M.T.Montgomery,2012:Impact of physics repre-sentations in the HWRFX on simulated hurricane structure and pressure-w ind relationships.Mon.Wea.Rev.,140,3278–3299.

    Beljaars,A.C.M.,and P.Viterbo,1998:Role of the boundary layer in a numerical weather prediction model.Clear and Cloudy Boundary Layers,A.A.M.Holtslag,and P.G. Duynkerke,Eds.,Royal Netherlands Academy of Arts and Sciences,Amsterdam,287–304.

    Bell,M.M.,M.T.Montgomery,and K.A.Emanuel,2012:Air-sea enthalpy and momentum exchange atmajor hurricane w ind speeds observed during CBLAST.J.Atmos.Sci.,69,3197–3222.

    Black,P.G.,and Coauthors,2007:Air-seaexchange in hurricanes: Synthesis of observations from the coupled boundary layer air-sea transferexperiment.Bull.Amer.Meteor.Soc.,88,357–374.

    Braun,S.A.,and W.-K.,Tao,2000:Sensitivity of high-resolution simulations of hurricane Bob(1991)to planetary boundary layer parameterizations.Mon.Wea.Rev.,128,3941–3961.

    Bryan,G.H.,2012:Effects of surface exchange coefficients and turbulence length scales on the intensity and structure of numerically simulated Hurricanes.Mon.Wea.Rev.,140,1125–1143.

    Carlson,T.N.,and F.E.Boland,1978:Analysis of urban-rural canopy using a surface heat flux/temperaturemodel.J.Appl. Meteor.,17,998–1013.

    Charnock,H.,1955:Wind stresson awatersurface.Quart.J.Roy. Meteor.Soc.,81,639–640.

    Chen,S.H.,andW.-Y.Sun,2002:A one-dimensional time dependentcloudmodel.J.Meteor.Soc.Japan,80(1),99–118.

    Chien,F.-C.,and H.-C.Kuo,2011:On the extreme rainfallof Typhoon Morakot(2009).J.Geophys.Res.,116,D05104,doi: 10.1029/2010JD015092.

    Davis,C.,and Coauthors,2008:Prediction of landfalling hurricanes w ith the advanced hurricane WRFmodel.Mon.Wea. Rev.,136,1990–2005.

    DeCosmo,J.,K.B.Katsaros,S.D.Smith,R.J.Anderson,W.A. Oost,K.Bumke,and H.Chadw ick,1996:Air-sea exchange of water vapor and sensible heat:The Humidity Exchange over the Sea(HEXOS)results.J.Geophys.Res.,101,12 001–12 016.

    Donelan,M.A.,B.K.Haus,N.Reul,W.J.Plant,M.Stiassnie,H. C.Graber,O.B.Brown,and E.S.Saltzman,2004:On the lim iting aerodynam ic roughness of the ocean in very strong w ind.Geophys.Res.Lett.,31,L18306,doi:10.1029/2004GL 019460.

    Drennan,W.M.,J.A.Zhang,J.R.French,C.McCormick,and P. G.Black,2007:Turbulent fluxes in the hurricane boundary layer.Part II:Latentheat flux.J.Atmos.Sci.,64,1103–1115.

    Dudhia,J.,1989:Numerical study of convection observed during the w intermonsoon experiment using amesoscale two dimensionalmodel.J.Atmos.Sci.,46,3077–3107.

    Dudhia,J.,and Coauthors,2008:Prediction of Atlantic tropical cyclonesw ith the Advanced HurricaneWRF(AHW)model. 28th Conf.on Hurricanes and Tropical Meteorology,Orlando,Florida,Preprint18A.2.

    Emanuel,K.A.,1986:An air-sea interaction theory for tropical cyclones.Part I:Steady-statemaintenance.J.Atmos.Sci.,43, 585–605.

    Emanuel,K.A.,1995:Sensitivity of tropical cyclones to surface exchange coefficients and a revised steady-statemodel incorporating eye dynamics.J.Atmos.Sci.,52,3969–3976.

    French,J.R.,W.M.Drennan,J.A.Zhang,and P.G.Black,2007: Turbulent fluxes in the hurricane boundary layer.Part I:Momentum flux.J.Atmos.Sci.,64,1089–1102.

    Garratt,J.R.,1992:The Atmospheric Boundary Layer.Cambridge University Press,316 pp.

    Geernaert,G.L.,K.L.Davidson,S.E.Larsen,and T.Mikkelsen, 1988:Wind stressmeasurements during the Tower Ocean Wave and Radar Dependence Experiment.J.Geophys.Res., 93,13 913–13 923.

    Green,B.W.,and F.Q.Zhang,2013:Impacts of air-sea flux parameterizations on the intensity and structure of tropical cyclones.Mon.Wea.Rev.,141,2308–2324.

    Hall,J.D.,M.Xue,L.K.Ran,L.M.Leslie,2013:High-resolution modeling of Typhoon Morakot(2009):Vortex Rossbywaves and their role in extreme precipitation over Taiwan.J.Atmos. Sci.,70,163–186.

    Haus,B.K.,D.Jeong,M.A.Donelan,J.A.Zhang,and I. Savelyev,2010:Relative rates of sea-air heat transfer and frictional drag in very high w inds.Geophys.Res.Lett.37, L07802,doi:10.1029/2009GL042206.

    Hong,S.Y.,Y.Noh,and J.Dudhia,2006:A new verticaldiffusion packagew ith an explicit treatmentof entrainmentprocesses. Mon.Wea.Rev.,134(9),2318–2341.

    Hosomi,T.,2005:Implementation of targeted moisture diffusion for the JMA Regional SpectralModel(RSM).CAS/JSC WGNE Research Activities in Atmospheric and Oceanic Modelling/WMO.,35,7–8.

    Huang,H.-L.,M.-J.Yang,and C.-H.Sui,2014:Waterbudgetand precipitationefficiency of TyphoonMorakot(2009).J.Atmos. Sci.,71,112–129.

    JMA,2007:Outline of the operational numericalweather prediction at the Japan Meteorological Agency.Appendix toWMO NumericalWeather Prediction Progress Report,Japan MeteorologicalAgency,Tokyo,194 pp.

    Large,W.G.,and S.Pond,1981:Open ocean momentum flux measurements in moderate to strong w inds.J.Phys. Oceanogr.,11,324–336.

    Li,X.,J.Ming,Y.Wang,K.Zhao,and M.Xue,2013:Assim ilation of T-TREC-retrieved w ind data w ith WRF 3DVAR for the short-term forecasting of typhoon Meranti(2010)near landfall.J.Geophys.Res.,118,10 361–10 375.

    Lin,Y.-L.,R.D.Farley,and H.D.Orville,1983:Bulk parameterization of the snow field in a cloud model.J.Climate Appl. Meteor.,22(6),1065–1092.

    Liu,Y.B.,D.-L.Zhang,and M.K.Yau,1999:Amultiscalenumerical study of Hurricane Andrew(1992).Part II:Kinematics and inner-core structures.Mon.Wea.Rev.,127,2597–2616.

    Malkus,J.S.,and H.Riehl,1960:On the dynam ics and energy transformations in steady-state hurricanes.Tellus,12,1–20.

    Ming,J.,Y.Q.Ni,and X.Y.Shen,2009:The dynam ical characteristicsand wave structureof typhoon Rananim(2004).Adv. Atmos.Sci.,26,523–542,doi:10.1007/s00376-009-0523-0.

    Ming,J.,J.J.Song,B.J.Chen,and K.F.Wang,2012:Boundary layer structure in typhoon Saomai(2006):Understanding the effects of exchange coefficient.J.Trop.Meteor.,18(2),195–206.

    Mlawer,E.J.,S.J.Taubman,P.D.Brown,M.J.Iacono,and S.A. Clough,1997:Radiative transfer for inhomogeneous atmosphere:RRTM,a validated correlated-k model for the longwave.J.Geophys.Res.,102,16 663–16 682.

    Montgomery,M.T.,R.K.Sm ith,and S.V.Nguyen,2010:Sensitivity of tropical-cyclonemodels to the surface drag coef-ficient.Quart.J.Roy.Meteor.Soc.,136,1945–1953,doi: 10.1002/qj.702.

    Noh,Y.,W.G.Cheon,S.Y.Hong,and S.Raasch,2003:Improvementof the K-profi lemodel for the planetary boundary layer based on large eddy simulation data.Bound.-LayerMeteor., 107(2),401–427.

    Nolan,D.S.,J.A.Zhang,and D.P.Stern,2009:Evaluation of planetary boundary layer parameterizations in tropical cyclonesby comparison of in situ dataand high-resolution simulationsofHurricane Isabel(2003).Part I:Initialization,maximum w inds,and the outer-core boundary layer structure. Mon.Wea.Rev.,137,3651–3674.

    Ooyama,K.,1969:Numericalsimulation of the life cycle of tropicalcyclones.J.Atmos.Sci.,26,3–40.

    Powell,M.D.,P.J.Vickery,and T.A.Reinhold,2003:Reduced drag coefficient forhighw ind speedsin tropicalcyclones.Nature,422,279–283.

    Rosenthal,S.L.,1971:The response of a tropical cyclonemodel to variations in boundary layer parameters,initial conditions, lateralboundary conditionsand domain size.Mon.Wea.Rev., 99,767–777.

    Rotunno,R.,and K.A.Emanuel,1987:An air-sea interaction theory for tropical cyclones.Part II:Evolutionary study using a nonhydrostaticaxisymmetricnumericalmodel.J.Atmos.Sci., 44,542–561.

    Schwartz,C.S.,Z.Q.Liu,Y.S.Chen,and X.Y.Huang,2012:Impactof assimilatingmicrowave radiancesw ith a limited-area ensemble data assim ilation system on forecasts of Typhoon Morakot.Wea.Forecasting,27,424–437.

    Skamarock,W.C.,and Coauthors,2008:Description of the advanced researchWRFversion 3,Rep.NCAR/TN-475++STR, Natl.Cent.for Atmos.Res.,Boulder,Colo.,113 pp.

    Sm ith,R.K.,M.T.Montgomery,and N.Van Sang,2009:Tropicalcyclonespin-up revisited.Quart.J.Roy.Meteor.Soc.,135, 1321–1335.

    Sm ith,S.D.,1980:Wind stress and heat flux over the ocean in gale forcew inds.J.Phys.Oceanogr.,10,709–726.

    Stern,D.P.,and D.S.Nolan,2012:On the height of the warm core in tropical cyclones.J.Atmos.Sci.,69,1657–1680.

    Stull,R.B.,1988:An Introduction to Boundary LayerMeteorology.K luwer Academ ic,666 pp.

    Wang,C.-C.,H.-C.Kuo,Y.H.Chen,H.-L.Huang,C.-H.Chung, and K.Tsuboki,2012:Effects of asymmetric latent heating on typhoonmovementcrossing Taiwan:The caseofMorakot (2009)w ith extreme rainfall.J.Atmos.Sci.,69,3172–3196.

    Wang,C.-C.,H.-C.Kuo,T.-C.Yeh,C.-H.Chung,Y.H.Chen, S.-Y.Huang,Y.W.Wang,and C.-H.Liu,2013:Highresolution quantitative precipitation forecastsand simulations by the Cloud-Resolving Storm Simulator(CReSS)for Typhoon Morakot(2009).J.Hydrol.,506,26–41.

    Xie,B.G.,and F.Q.Zhang,2012:Impacts of typhoon track and island topography on theheavy rainfalls in Taiwan associated w ith Morakot(2009).Mon.Wea.Rev.,140,3379–3394.

    Zhang,D.L.,Y.B.Liu,and M.K.Yau,2001:Amultiscalenumerical study of Hurricane Andrew(1992).Part IV:Unbalanced flows.Mon.Wea.Rev.,129,92–107.

    Zhang,D.L.,and H.Chen,2012:Importance of the upperlevelwarm core in the rapid intensification of a tropical cyclone.Geophys.Res.Lett.,39,L02806,doi:10.1029/2011GL 050578.

    Zhang,F.Q.,Y.H.Weng,Y.-H.Kuo,J.S.Whitaker,and B.G. Xie,2010:Predicting Typhoon Morakot’s catastrophic rainfallw ith a convection-perm ittingmesoscaleensemblesystem. Wea.Forecasting,25,1816–1825.

    Zhang,J.A.,P.G.Black,J.R.French,andW.M.Drennan,2008: First directmeasurements of enthalpy flux in the hurricane boundary layer:the CBLAST results.Geophys.Res.Lett., 35(14),L14813,doi:10.1029/2008GL034374.

    Zhang,J.A.,R.F.Rogers,D.S.Nolan,and F.D.Marks,2011: On the characteristic heightscalesof the hurricane boundary layer.Mon.Wea.Rev.,139,2523–2535.

    Zhang,J.A.,S.Gopalakrishnan,F.D.Marks,R.F.Rogers,and V. Tallapragada,2012:A developmental framework for improving hurricane model physical parameterizations using aircraftobservations.Trop.Cycl.Res.Rev.,1(4),419–429,doi: 10.6057/2012TCRR04.01.

    2April2015;revised 2 June2015;accepted 16 June2015)

    :Ming,J.,and J.A.Zhang,2015:Effects of surface flux parameterization on the numerically simulated intensity and structureof Typhoon Morakot(2009).Adv.Atmos.Sci.,33(1),58–72,

    10.1007/s00376-015-4202-z.

    ?Corresponding author:JieMING

    Email:jming@nju.edu.cn

    国产区一区二久久| 乱人伦中国视频| 久久久国产欧美日韩av| 琪琪午夜伦伦电影理论片6080| 丝袜人妻中文字幕| 精品久久久久久成人av| 不卡av一区二区三区| www.自偷自拍.com| 亚洲av电影在线进入| 亚洲人成电影观看| 国产单亲对白刺激| 欧美国产精品va在线观看不卡| 天天躁夜夜躁狠狠躁躁| 欧美日韩中文字幕国产精品一区二区三区 | 黑人猛操日本美女一级片| 黑人猛操日本美女一级片| 久久人妻av系列| 国产aⅴ精品一区二区三区波| 精品久久久久久成人av| 亚洲精品国产色婷婷电影| 夜夜躁狠狠躁天天躁| av国产精品久久久久影院| 午夜福利影视在线免费观看| 手机成人av网站| 超碰97精品在线观看| 波多野结衣高清无吗| 午夜福利欧美成人| 黄色女人牲交| 在线天堂中文资源库| 99精国产麻豆久久婷婷| 亚洲国产精品一区二区三区在线| 母亲3免费完整高清在线观看| 日韩高清综合在线| 国产av又大| 亚洲一区二区三区欧美精品| 日韩成人在线观看一区二区三区| 精品乱码久久久久久99久播| 国产一区二区在线av高清观看| 久久久精品欧美日韩精品| 日本a在线网址| 欧美日韩视频精品一区| 狂野欧美激情性xxxx| 校园春色视频在线观看| 黄片播放在线免费| 欧美午夜高清在线| 成年人免费黄色播放视频| 露出奶头的视频| 老鸭窝网址在线观看| bbb黄色大片| 精品久久久久久,| 婷婷六月久久综合丁香| 亚洲精品国产精品久久久不卡| 久久久水蜜桃国产精品网| 中文字幕人妻丝袜制服| 亚洲色图 男人天堂 中文字幕| 老司机福利观看| 欧美最黄视频在线播放免费 | 老熟妇仑乱视频hdxx| 搡老岳熟女国产| 岛国视频午夜一区免费看| 80岁老熟妇乱子伦牲交| 亚洲第一青青草原| 亚洲,欧美精品.| 国产xxxxx性猛交| 日本 av在线| www.999成人在线观看| 精品久久久久久久毛片微露脸| 天堂中文最新版在线下载| 亚洲黑人精品在线| 久久人妻av系列| 亚洲精品在线美女| 淫秽高清视频在线观看| 欧美日韩视频精品一区| 一级片'在线观看视频| 怎么达到女性高潮| 国产成+人综合+亚洲专区| 国产无遮挡羞羞视频在线观看| 在线观看一区二区三区激情| 亚洲精品久久成人aⅴ小说| 久久天堂一区二区三区四区| 女生性感内裤真人,穿戴方法视频| 黄色视频,在线免费观看| 最近最新免费中文字幕在线| 这个男人来自地球电影免费观看| 国产欧美日韩一区二区精品| 桃色一区二区三区在线观看| 欧美一级毛片孕妇| 丁香六月欧美| 一级a爱片免费观看的视频| 狠狠狠狠99中文字幕| 欧洲精品卡2卡3卡4卡5卡区| 欧美日韩黄片免| 国产黄a三级三级三级人| bbb黄色大片| 丁香六月欧美| 亚洲第一av免费看| 国产成人影院久久av| 丰满迷人的少妇在线观看| 妹子高潮喷水视频| 国产精品自产拍在线观看55亚洲| 高清av免费在线| 九色亚洲精品在线播放| 国产黄a三级三级三级人| 精品国产一区二区久久| 亚洲午夜精品一区,二区,三区| 黄色女人牲交| 99久久国产精品久久久| 欧美不卡视频在线免费观看 | 国产精品一区二区在线不卡| 韩国精品一区二区三区| av网站免费在线观看视频| 级片在线观看| ponron亚洲| 久久久精品欧美日韩精品| 国产成人欧美在线观看| 黄网站色视频无遮挡免费观看| 十八禁网站免费在线| 日韩有码中文字幕| 欧美激情久久久久久爽电影 | 又大又爽又粗| 99久久99久久久精品蜜桃| 国产欧美日韩一区二区三区在线| 乱人伦中国视频| 精品高清国产在线一区| 老司机亚洲免费影院| 中文亚洲av片在线观看爽| 韩国精品一区二区三区| 亚洲成人国产一区在线观看| 欧美日韩乱码在线| 这个男人来自地球电影免费观看| 欧美中文日本在线观看视频| 99精品在免费线老司机午夜| 老司机在亚洲福利影院| 日本五十路高清| 欧美黄色片欧美黄色片| 久久人妻熟女aⅴ| 大香蕉久久成人网| 最新在线观看一区二区三区| 久久人人爽av亚洲精品天堂| 成年版毛片免费区| 人人妻,人人澡人人爽秒播| 亚洲国产精品999在线| 搡老岳熟女国产| 国产一区二区在线av高清观看| 国产精品偷伦视频观看了| 黄色怎么调成土黄色| 亚洲中文字幕日韩| 国产高清激情床上av| 欧美日韩乱码在线| 女警被强在线播放| 级片在线观看| 757午夜福利合集在线观看| 久久久久久人人人人人| 日韩欧美一区二区三区在线观看| 久久狼人影院| 男女之事视频高清在线观看| 免费看十八禁软件| 免费少妇av软件| 日本黄色日本黄色录像| 国产单亲对白刺激| 日韩中文字幕欧美一区二区| 一夜夜www| 国产片内射在线| 亚洲国产欧美日韩在线播放| 中文字幕人妻丝袜制服| av国产精品久久久久影院| 亚洲自拍偷在线| 久久精品91无色码中文字幕| 久久精品国产综合久久久| 日韩精品免费视频一区二区三区| 一级作爱视频免费观看| 国产精品免费视频内射| 天天添夜夜摸| 最近最新中文字幕大全电影3 | 高清av免费在线| 一区二区三区激情视频| 国产伦人伦偷精品视频| 人人澡人人妻人| av中文乱码字幕在线| a级毛片在线看网站| 欧美最黄视频在线播放免费 | 在线观看免费午夜福利视频| 精品福利观看| 午夜福利在线免费观看网站| 制服人妻中文乱码| 日本五十路高清| 18美女黄网站色大片免费观看| 亚洲精品av麻豆狂野| 国产又爽黄色视频| 午夜老司机福利片| 黑人巨大精品欧美一区二区mp4| 女性生殖器流出的白浆| 国内久久婷婷六月综合欲色啪| 国产精品综合久久久久久久免费 | 如日韩欧美国产精品一区二区三区| 亚洲午夜理论影院| 日韩欧美三级三区| 久久精品国产99精品国产亚洲性色 | 人妻丰满熟妇av一区二区三区| av中文乱码字幕在线| 一边摸一边抽搐一进一小说| 中文字幕高清在线视频| 久热这里只有精品99| 国产国语露脸激情在线看| 欧美日韩精品网址| 午夜福利免费观看在线| 国产人伦9x9x在线观看| 国产一区在线观看成人免费| 热99re8久久精品国产| av免费在线观看网站| 亚洲精品一二三| 高清毛片免费观看视频网站 | 老司机在亚洲福利影院| 亚洲国产看品久久| 69av精品久久久久久| 操美女的视频在线观看| 窝窝影院91人妻| 黄网站色视频无遮挡免费观看| 91麻豆精品激情在线观看国产 | 岛国在线观看网站| 在线视频色国产色| 免费女性裸体啪啪无遮挡网站| 亚洲免费av在线视频| 黑人操中国人逼视频| 亚洲av成人不卡在线观看播放网| 黄色视频,在线免费观看| 国产高清videossex| 99久久人妻综合| a在线观看视频网站| 好看av亚洲va欧美ⅴa在| 美女高潮喷水抽搐中文字幕| 国产精品永久免费网站| 精品一区二区三卡| 国产精品亚洲av一区麻豆| 国产一区二区三区视频了| 男女下面插进去视频免费观看| 可以在线观看毛片的网站| 亚洲精品久久成人aⅴ小说| bbb黄色大片| 在线观看免费视频日本深夜| 日韩视频一区二区在线观看| 国产精品电影一区二区三区| 午夜福利,免费看| 婷婷丁香在线五月| 男女午夜视频在线观看| 视频区图区小说| 久久国产精品影院| cao死你这个sao货| 欧美 亚洲 国产 日韩一| 岛国视频午夜一区免费看| 婷婷丁香在线五月| 亚洲成人精品中文字幕电影 | 在线观看66精品国产| 国产精品av久久久久免费| 丰满迷人的少妇在线观看| av天堂在线播放| 国产伦人伦偷精品视频| 国产97色在线日韩免费| 国产亚洲精品综合一区在线观看 | 在线观看免费视频网站a站| 久久人妻福利社区极品人妻图片| 久久草成人影院| 淫妇啪啪啪对白视频| 一级黄色大片毛片| 精品人妻在线不人妻| 午夜亚洲福利在线播放| 欧美日韩福利视频一区二区| 韩国av一区二区三区四区| 欧美激情 高清一区二区三区| 亚洲成国产人片在线观看| 老熟妇仑乱视频hdxx| www.www免费av| 欧美国产精品va在线观看不卡| 欧美激情极品国产一区二区三区| 黄色a级毛片大全视频| 精品人妻在线不人妻| 亚洲一区高清亚洲精品| 午夜免费成人在线视频| 女警被强在线播放| 欧美一区二区精品小视频在线| 狠狠狠狠99中文字幕| 男女高潮啪啪啪动态图| 天堂俺去俺来也www色官网| 日本黄色日本黄色录像| 三级毛片av免费| 免费一级毛片在线播放高清视频 | 九色亚洲精品在线播放| 久久精品aⅴ一区二区三区四区| 大香蕉久久成人网| 黑人欧美特级aaaaaa片| 超碰成人久久| 99re在线观看精品视频| 美女高潮到喷水免费观看| 黄片小视频在线播放| 一级,二级,三级黄色视频| 制服人妻中文乱码| 韩国精品一区二区三区| 欧美日韩亚洲综合一区二区三区_| 成人亚洲精品av一区二区 | 久久人妻熟女aⅴ| 精品欧美一区二区三区在线| 97人妻天天添夜夜摸| 免费高清视频大片| 美女 人体艺术 gogo| 午夜91福利影院| 丁香欧美五月| 久久久久国内视频| 女警被强在线播放| 国产99久久九九免费精品| 乱人伦中国视频| 国产野战对白在线观看| 午夜精品久久久久久毛片777| 国产亚洲精品久久久久久毛片| 日日干狠狠操夜夜爽| 好看av亚洲va欧美ⅴa在| 90打野战视频偷拍视频| 真人做人爱边吃奶动态| 亚洲av日韩精品久久久久久密| 人妻丰满熟妇av一区二区三区| 18禁国产床啪视频网站| 久久国产精品人妻蜜桃| 免费在线观看完整版高清| 丰满人妻熟妇乱又伦精品不卡| 日韩 欧美 亚洲 中文字幕| 国产精品亚洲av一区麻豆| 国产精品影院久久| 两人在一起打扑克的视频| 天天躁夜夜躁狠狠躁躁| 精品国产一区二区久久| 精品人妻1区二区| 热99re8久久精品国产| 一进一出抽搐gif免费好疼 | 国产成人系列免费观看| 欧美中文日本在线观看视频| 女人高潮潮喷娇喘18禁视频| 国产麻豆69| 老熟妇仑乱视频hdxx| 一区在线观看完整版| 国产精品亚洲一级av第二区| 嫁个100分男人电影在线观看| 侵犯人妻中文字幕一二三四区| 亚洲一码二码三码区别大吗| 亚洲va日本ⅴa欧美va伊人久久| 露出奶头的视频| 精品高清国产在线一区| 黄色视频,在线免费观看| 精品卡一卡二卡四卡免费| 视频区欧美日本亚洲| 日日摸夜夜添夜夜添小说| 一二三四社区在线视频社区8| 搡老乐熟女国产| 国产欧美日韩一区二区三| 国产精品 欧美亚洲| 欧美黑人精品巨大| 中文字幕最新亚洲高清| 人人妻人人爽人人添夜夜欢视频| 人人澡人人妻人| 免费在线观看完整版高清| 欧美大码av| 久久天堂一区二区三区四区| 精品久久久久久久久久免费视频 | 亚洲欧美精品综合久久99| 亚洲国产精品999在线| 久久久国产成人精品二区 | 成人18禁高潮啪啪吃奶动态图| 男女做爰动态图高潮gif福利片 | 又大又爽又粗| 热re99久久国产66热| 国产真人三级小视频在线观看| 午夜精品久久久久久毛片777| 国产精品99久久99久久久不卡| 99精品欧美一区二区三区四区| 很黄的视频免费| 老熟妇乱子伦视频在线观看| 亚洲久久久国产精品| 婷婷丁香在线五月| 精品一区二区三区视频在线观看免费 | 欧美黑人欧美精品刺激| 亚洲一区高清亚洲精品| 999精品在线视频| 午夜亚洲福利在线播放| 免费av中文字幕在线| 12—13女人毛片做爰片一| 欧美在线一区亚洲| 一级黄色大片毛片| 午夜精品在线福利| 18禁观看日本| 性少妇av在线| 日本欧美视频一区| 熟女少妇亚洲综合色aaa.| av欧美777| 免费av中文字幕在线| 亚洲激情在线av| 久久热在线av| 国产一区二区激情短视频| 男男h啪啪无遮挡| 91成年电影在线观看| 午夜日韩欧美国产| 51午夜福利影视在线观看| 久久久久国产精品人妻aⅴ院| 亚洲一码二码三码区别大吗| 国产乱人伦免费视频| 18美女黄网站色大片免费观看| 久久香蕉精品热| 香蕉国产在线看| 成人av一区二区三区在线看| 性少妇av在线| 久久国产乱子伦精品免费另类| 香蕉国产在线看| 水蜜桃什么品种好| 国产黄a三级三级三级人| 99久久综合精品五月天人人| 亚洲五月天丁香| 中文字幕人妻熟女乱码| 香蕉久久夜色| 男女床上黄色一级片免费看| 欧美人与性动交α欧美精品济南到| 99久久国产精品久久久| 亚洲欧美日韩另类电影网站| 国产欧美日韩精品亚洲av| 99久久精品国产亚洲精品| 手机成人av网站| avwww免费| 91国产中文字幕| 国产99久久九九免费精品| 亚洲国产看品久久| 免费看十八禁软件| 男男h啪啪无遮挡| 一区二区日韩欧美中文字幕| 夫妻午夜视频| 88av欧美| 国产aⅴ精品一区二区三区波| 国产精品香港三级国产av潘金莲| 人妻久久中文字幕网| 久久狼人影院| 黑人猛操日本美女一级片| 在线观看一区二区三区| 亚洲精品国产一区二区精华液| 国产不卡一卡二| 欧美一级毛片孕妇| 久久精品aⅴ一区二区三区四区| 精品日产1卡2卡| 免费少妇av软件| 女人被狂操c到高潮| 涩涩av久久男人的天堂| 精品久久蜜臀av无| 精品一品国产午夜福利视频| 免费久久久久久久精品成人欧美视频| 另类亚洲欧美激情| 一级毛片女人18水好多| 午夜福利,免费看| 日本三级黄在线观看| 少妇被粗大的猛进出69影院| 午夜成年电影在线免费观看| 国产91精品成人一区二区三区| 午夜免费鲁丝| 不卡一级毛片| 在线观看免费日韩欧美大片| 老司机靠b影院| 真人一进一出gif抽搐免费| 欧美一区二区精品小视频在线| bbb黄色大片| 神马国产精品三级电影在线观看 | 老司机午夜福利在线观看视频| 两个人免费观看高清视频| 9色porny在线观看| 大陆偷拍与自拍| 亚洲色图av天堂| 欧美成狂野欧美在线观看| 久久欧美精品欧美久久欧美| 欧美日韩亚洲国产一区二区在线观看| 亚洲激情在线av| 黄色毛片三级朝国网站| 在线看a的网站| 天堂俺去俺来也www色官网| 热99re8久久精品国产| 国产一区在线观看成人免费| 人人澡人人妻人| 免费高清在线观看日韩| 国产精品1区2区在线观看.| 久9热在线精品视频| 国产av在哪里看| 欧美日韩亚洲国产一区二区在线观看| 欧美黄色片欧美黄色片| 国产熟女xx| 久久精品91无色码中文字幕| 亚洲专区中文字幕在线| 熟女少妇亚洲综合色aaa.| 亚洲久久久国产精品| 精品久久久精品久久久| 99riav亚洲国产免费| 嫩草影视91久久| 自拍欧美九色日韩亚洲蝌蚪91| 精品乱码久久久久久99久播| 欧美日韩国产mv在线观看视频| 日韩人妻精品一区2区三区| 丝袜在线中文字幕| 久久精品国产99精品国产亚洲性色 | 中文字幕av电影在线播放| 欧美黑人欧美精品刺激| 正在播放国产对白刺激| 国产欧美日韩一区二区三| 午夜久久久在线观看| 久久亚洲精品不卡| 日本三级黄在线观看| 精品人妻在线不人妻| 午夜久久久在线观看| 国产一区在线观看成人免费| 首页视频小说图片口味搜索| 国产精品影院久久| bbb黄色大片| 精品电影一区二区在线| 久久香蕉激情| 亚洲全国av大片| 人人妻人人澡人人看| 99久久久亚洲精品蜜臀av| 亚洲色图 男人天堂 中文字幕| 精品少妇一区二区三区视频日本电影| 午夜福利在线观看吧| 国产精品亚洲一级av第二区| 国产成人影院久久av| 日韩中文字幕欧美一区二区| 高清毛片免费观看视频网站 | 国内毛片毛片毛片毛片毛片| 亚洲五月色婷婷综合| x7x7x7水蜜桃| 国产精品国产av在线观看| 免费一级毛片在线播放高清视频 | 国产精品永久免费网站| 90打野战视频偷拍视频| 少妇 在线观看| 精品无人区乱码1区二区| 两个人免费观看高清视频| 一区二区三区国产精品乱码| 老鸭窝网址在线观看| 亚洲 国产 在线| 多毛熟女@视频| 久久久久久久午夜电影 | 中文字幕人妻熟女乱码| 男女下面进入的视频免费午夜 | av福利片在线| 一区二区日韩欧美中文字幕| 国产成人精品久久二区二区免费| 亚洲一区二区三区欧美精品| 亚洲少妇的诱惑av| 国产在线精品亚洲第一网站| 精品免费久久久久久久清纯| 国产一区二区三区视频了| 亚洲欧美精品综合久久99| 亚洲熟妇熟女久久| 国产乱人伦免费视频| 国产精品秋霞免费鲁丝片| 一区二区日韩欧美中文字幕| 亚洲av成人一区二区三| 国产黄色免费在线视频| 亚洲自拍偷在线| 人妻丰满熟妇av一区二区三区| 日日干狠狠操夜夜爽| 女人被躁到高潮嗷嗷叫费观| 日韩欧美在线二视频| 久久久久久大精品| 国产精品98久久久久久宅男小说| 日韩一卡2卡3卡4卡2021年| 亚洲,欧美精品.| 不卡av一区二区三区| 69精品国产乱码久久久| 亚洲精品久久成人aⅴ小说| 色播在线永久视频| 大香蕉久久成人网| 欧美在线黄色| 咕卡用的链子| 国产精品成人在线| 日本欧美视频一区| 在线观看免费高清a一片| 一级,二级,三级黄色视频| 新久久久久国产一级毛片| 黑人欧美特级aaaaaa片| 亚洲av熟女| 黑人巨大精品欧美一区二区蜜桃| 制服人妻中文乱码| 亚洲精品中文字幕一二三四区| ponron亚洲| 亚洲欧洲精品一区二区精品久久久| 国产有黄有色有爽视频| 悠悠久久av| 成人18禁在线播放| 少妇裸体淫交视频免费看高清 | 国产精品久久视频播放| 激情在线观看视频在线高清| www.www免费av| 久久久久久免费高清国产稀缺| 久久久国产成人免费| 韩国精品一区二区三区| 欧美日韩福利视频一区二区| 国产三级在线视频| 午夜久久久在线观看| 亚洲久久久国产精品| 午夜91福利影院| 一级片免费观看大全| 免费在线观看影片大全网站| 免费不卡黄色视频| 男女下面进入的视频免费午夜 | 精品久久久久久,| www.熟女人妻精品国产| av天堂久久9| a级片在线免费高清观看视频| 黄色怎么调成土黄色| 大陆偷拍与自拍| 亚洲一区中文字幕在线| 视频在线观看一区二区三区| 欧美亚洲日本最大视频资源| 久久这里只有精品19| 一区福利在线观看| xxxhd国产人妻xxx| 亚洲九九香蕉|