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

    Study on Multi-Scale Blending Initial Condition Perturbations for a Regional Ensemble Prediction System

    2015-06-09 21:37:22ZHANGHanbinCHENJingZHIXiefeiWANGYiandWANGYanan
    Advances in Atmospheric Sciences 2015年8期

    ZHANG Hanbin,CHEN Jing,ZHI Xiefei,WANG Yi,and WANG Yanan

    1College of Atmospheric Science,Nanjing University of Information Science and Technology,Nanjing 210044

    2Numerical Weather Prediction Center,China Meteorological Administration,Beijing 100081

    3Key Laboratory of Meteorological Disaster,Ministry of Education, Nanjing University of Information Science and Technology,Nanjing 210044

    4Earth and Atmospheric Sciences,University of Nebraska Lincoln,Lincoln,Nebraska,68588,USA5Center of Meteorological Service of Zhejiang,Hangzhou 310017

    Study on Multi-Scale Blending Initial Condition Perturbations for a Regional Ensemble Prediction System

    ZHANG Hanbin1,CHEN Jing?2,ZHI Xiefei3,WANG Yi4,and WANG Yanan5

    1College of Atmospheric Science,Nanjing University of Information Science and Technology,Nanjing 210044

    2Numerical Weather Prediction Center,China Meteorological Administration,Beijing 100081

    3Key Laboratory of Meteorological Disaster,Ministry of Education, Nanjing University of Information Science and Technology,Nanjing 210044

    4Earth and Atmospheric Sciences,University of Nebraska Lincoln,Lincoln,Nebraska,68588,USA5Center of Meteorological Service of Zhejiang,Hangzhou 310017

    An initial conditions(ICs)perturbation method was developed with the aim to improve an operational regional ensemble prediction system(REPS).Three issues were identified and investigated:(1)the impacts of perturbation scale on the ensemble spread and forecast skill of the REPS;(2)the scale characteristic of the IC perturbations of the REPS;and(3)whether the REPS’s skill could be improved by adding large-scale information to the IC perturbations.Numerical experiments were conducted to reveal the impact of perturbation scale on the ensemble spread and forecast skill.The scales of IC perturbations from the REPS and an operational global ensemble prediction system(GEPS)were analyzed.A“multi-scale blending”(MSB)IC perturbation scheme was developed,and the main findings can be summarized as follows:The growth rates of the ensemble spread of the REPS are sensitive to the scale of the IC perturbations;the ensemble forecast skills can benefit from large-scale perturbations;the global ensemble IC perturbations exhibit more power at larger scales,while the regional ensemble IC perturbations contain more power at smaller scales;the MSB method can generate IC perturbations by combining the small-scale component from the REPS and the large-scale component from the GEPS;the energy norm growth of the MSB-generated perturbations can be appropriate at all forecast lead times;and the MSB-based REPS shows higher skill than the original system,as determined by ensemble forecast verification.

    regional ensemble prediction system,spectral analysis,multi-scale blending,initial condition perturbations

    1.Introduction

    The errors related to initial conditions(ICs),numerical models,as well as the chaotic nature of the atmosphere(Lorenz,1969),lead to great uncertainty in numerical weather prediction.It is desirable for numerical weather prediction to describe the evolution of the model atmosphere as a probability density function,rather than provide a deterministic result.It was in this context that ensemble prediction theory(Leith,1974)was proposed.Ensemble prediction systems were originally implemented in the early 1990s at the National Centers for Environmental Prediction(NCEP;Toth and Kalnay,1993)and at the European Center for Medium-Range Weather Forecasts(ECMWF;Molteni et al.,1996). Various ensemble perturbation schemes have since been developed to address the uncertainties in global ensemble prediction systems(GEPSs)related to ICs,such as singular vectors(SVs;Buizza and Palmer,1995;Molteni et al.,1996), breeding growing modes(BGMs;Toth and Kalnay,1993; 1997)and perturbed observation(Houtekamer et al.,1996; Buizza et al.,2005).Furthermore,with an increased emphasis on representing the analysis probability density function in the initial analysis state,the ensemble transform Kalman filter(ETKF;Wang and Bishop,2003;Bowler et al.,2008; Ma et al.,2008),ensemble transform,and ensemble transform with rescaling(Wei et al.,2008),have been frequently used in recent years.

    The development of regional ensemble prediction lags behind that of global ensemble prediction.Nevertheless, the need for regional ensemble prediction systems(REPSs) to provide mesoscale severe weather prediction is clear,so increasing numbers of related studies have been conducted (e.g.,Du and Tracton,2001;Grimit and Mass,2002;Duan etal.,2012).As the probability distributions of various sources of errors are complicated,it is difficult to address the uncertainties in the ICs,model physics and lateral boundary conditions(LBCs)of regional numerical weather prediction(Chen et al.,2004,2005a).An important aspect of REPS research is exploring the effective methods to generate IC perturbations. One possibility in this regard is to apply dynamical downscaling to a GEPS(Marsigli et al.,2005;Frogner et al.,2006; Bowler et al.,2009).However,despite this method being attractive for its simplicity and good performance,small-scale uncertainties cannot be explicitly represented through downscaling(Wang et al.,2011).Other studies have generated IC perturbations for REPSs by using regional versions of traditional methods(i.e.,BGMs,SVs,ETKF etc.)(Stensrud et al.,1999;Du et al.,2003;Chen et al.,2005b;Li et al.,2008; Wang et al.,2011;Zhang et al.,2014a),and demonstrated that these methods can generate some ensemble spread and benefit the forecast skill of the REPS.While it remains uncertain whether or not these regional IC perturbation generators are superior to dynamical downscaling(Bowler and Mylne, 2009;Saito et al.,2011),there is no doubt that they can produce more information on small/mesoscale uncertainties,and this information is particularly useful for the forecasting of high-impact weather,quantitative precipitation,and local severe convective weather(Chen et al.,2005b;Stensrud and Yussouf,2007).However,when using these regionalversions of traditional IC perturbation generators,it was found that the mismatches in the LBCs can be a significant source of spurious perturbations,and thus the use of a“blending”scheme was proposed(Wang et al.,2011;Caron,2013;Du et al., 2014).This method can combine the small-scale component from the IC perturbations of a regional ensemble with the large-scale component from the IC perturbations of a global ensemble,and has proven to be effective in ameliorating the problem of mismatch in the LBCs(Caron,2013).

    The Numerical Weather Prediction Center of the China Meteorological Administration(CMA)has been routinely running a REPS since 2014.The system calculates IC perturbations by using a regional version of ETKF.Although this REPS has proven to be feasible in providing probabilistic forecasts of mesoscale phenomena,the system still faces several obstacles,such as a lack of ensemble spread,obvious systematic bias,and imperfect prediction of mesoscale uncertainty(Zhang et al.,2014b).Therefore,further improvement is needed.Bowler et al.(2009)reported that using downscaling,rather than the regional ETKF approach,to produce IC perturbations for REPSs can achieve better ensemble spread and forecast skill.As downscaling can create IC perturbations with more information at large scales,while regional perturbation generators can create ICperturbations with more information at small scales(Bowler and Mylne,2009),questions arise as to whether the perturbation scale is an important factor,or to what extent different scale perturbations can affect the REPS’spread and forecast skill.Addressing these issues could lead to improvements in current REPSs.

    In this study,we examine the impacts of perturbation scale on the ensemble spread and forecast skill of a REPS. Meanwhile,the scale characteristics of the IC perturbations generated by an operational ETKF are studied.We also explore the benefits of adding more large-scale information into the REPS IC perturbations,leading to the development of a multi-scale blending(MSB)perturbation method.The paper is structured as follows:section 2 introduces the model and data;section 3 compares the results of different numerical experiments that used different initial perturbation scale settings;and section 4 demonstrates the positive impacts of the MSB perturbation method on the ensemble spread and forecast skill of the REPS.A conclusion and discussion are provided in section 5.

    2.Model and data

    2.1.The regional ensemble prediction system

    2.1.1.Operational configuration

    2.1.2.Operational scheme

    GRAPES-REPS generates IC perturbations via the ETKF method,while the model uncertainty is represented by multiphysics.The LBCs of GRAPES-REPS are perturbed through coupling with T639-GEPS.

    2.2.Experimental data

    2.2.1.The T639-based global ensemble prediction system

    The background states and LBCs of GRAPES-REPS are provided by T639-GEPS.The resolution of this GEPS is T639L60(spectral triangular T639 with 60 vertical levels, corresponding to 30 km resolution),with ICs perturbed using BGMs.The BGM cycle interval is 12 h,initiated at 0000 and 0012 UTC each day with a forecast lead time of 15 days. A single control forecast is also initiated at 0006 and 0018 UTC to meet the requirements of GRAPES-REPS in terms of the background states and LBCs at these two initiation times. Model uncertainties are simulated using a stochastic physical parameterization scheme(Buizza et al.,1999).

    2.2.2.Data for verification

    The T639-GEPS analysis states corresponding to each forecast lead time are interpolated to a common regular 0.15?×0.15?grid to verify the variables at different pressure levels,while the precipitation isverified againstobservational accumulative precipitation from 2507 meteorologicalstations in China.

    3.Numerical experiments with multi-scale IC perturbations

    To investigate whether or not the perturbation scale is an important factor affecting the perturbation growth and forecast skill of the REPS,the results from a series of numerical experiments that used different IC perturbation scales are studied.

    3.1.Construction of IC perturbations with different scales

    A“scale selective perturbation generator”,which can construct IC perturbations with a particular scale,was developed with the following steps:

    Step 1:Generate a series of uniform random numbers ribetween-1.5 and 1.5 using a random number generator (Buizza et al.,1999).

    Step 2:Divide the regular grid domain into small domains(D1,D2,...,Di)of the same size(the size of the small domains corresponds to the required scale in a test),and then assign the random number riderived from Step 1 to the conjunctional grid points between the small domains as their values.Meanwhile,the values of grid points within the small domains are obtained through bilinear interpolation, and thereby a random number grid state characterized by a particular scale is formed,with values of all grids distributed smoothly.

    Step 3:Generate seven(half of the number of perturbed members)two-dimensionalrandom grid statesvia Step 2,and then multiply the random grid states by the statistical analysis error of the five variables(u,v,θ,πand q)at different model levels to obtain seven groups of perturbations for all variables at all levels.

    Step 4:Add the seven groups of perturbations to(or subtract from)the IC of the control,and then the IC perturbations of the first(or last)seven of the 14 ensemble members are obtained.Meanwhile,the IC perturbations for the 14 perturbed members are positive–negative paired.

    By using the“scale selective perturbation generator”, three tests with different scale settings(namely,R1,R2 and R3)were conducted.The scales for the three tests are listed in Table 1.For all the tests in this section,the model configuration was the same as the operational run.A test period of 10 consecutive days(5–15 August 2012)was conducted. The ensembles were initiated at 1200 UTC on each day,with a forecast length of 72 h.

    3.2.Analysis of results

    The characteristics of different types of perturbations areinvestigated.Here,we employ a total energy norm,which is appropriate for weather forecasting and data assimilation (Palmer et al.,1998).For one grid located at(i,j,k),the energy of the perturbation is computed from winds and temperature using the approximate energy norm defined as

    Table 1.Configuration of the IC perturbations for the three regional ensemble forecast tests.

    where u′,v′and T′are wind and temperature perturbations,cpis the specific heat and Tris the reference temperature(Wang and Bishop,2003;Wei et al.,2006).

    Figure 2 showsthe horizontaldistribution ofthe vertically averaged IC perturbation energy norm for the first ensemble member.It can be seen that,although the amplitudes of the three types of perturbations are similar,there are remarkable differences in the horizontal distribution pattern.The perturbation for R1 is structured with very small scale,while the perturbation scale for the R2 scheme is increased and the numberofperturbation centersisreduced,and the R3 scheme presents perturbations with the largest scale and fewest centers.

    To examine the vertical distributions and evolutions of the three types of perturbations,we averaged the energy norm of all grid points at each level.Figure 3 shows the resulting vertical distribution profiles of the perturbation energy norm for the R1,R2 and R3 ensembles.It is clear that the perturbation growth is quite different for the three schemes.The energy norm at all levels in the R1 ensemble is larger than that of the other two at the initial time(00 h).For the R1 scheme, the energy norm at 200 hPa can reach 2.3 J kg-1,while for R2 and R3 the corresponding values are about 2.15 J kg-1. For the 6 h forecast,the R1 scheme perturbations at most levels are sharply reduced to less than 1 J kg-1,whereas for the R2 and R3 ensembles the perturbations exhibit growth.As the forecast lead time increases,the perturbations of the three ensembles grow gradually.The R3 perturbations show the largest growth,as the 36 h energy magnitude at 250 hPa is 3.6 J kg-1.For the R2 ensemble,the perturbation growth is relatively smaller;the 36 h energy value at 250 hPa is 2.9 J kg-1.The perturbations for R1 are smallest,as the energy at all levels is even smaller than that at the initial time.For the total energy norm of all three schemes demonstrated in Fig. 3,it seems that the differences at the upper levels are dominated by kinetic energy,while the differences at the lower levels are dominated by internal energy(not shown).This is because the kinetic energy has larger magnitude at the upper levels and the internal energy has larger magnitude at the lower levels.

    The above result indicates that perturbations with larger scale exhibit more adequate perturbation growth.The reason for this behavior might be that larger scale perturbations can easily form“organized structure”,and perturbations with this“organized structure”can easily evolve into“flow dependent”perturbations that can develop well with the atmospheric flow (Toth and Kalnay,1993).Meanwhile,smaller scale perturbations are more like inertial gravity waves and are generally decaying(e.g.Lacarra and Talagrand,1988),and thus the perturbation growth would be constrained.This result indicates that the spread growth of the regional ensemble is sensitive to the perturbation scale,with large-scale perturbations more prone to growth.

    Verification of the three ensembles was performed(not shown).Results averaged over 10 days showed remarkabledifferences in the spread and root-mean-square error (RMSE).R3 maintained the largest spread at all forecast lead times,followed by R2 and then R1.Thisindicatesthatperturbations with largerscale favorthe growth ofensemble spread. As the current REPS is known for its lack of spread(Zhang et al.,2014a),the amplified spread due to the enlarged perturbation scale could also lead to RMSE reduction,thus enabling large-scale perturbations to produce a much better ensemble forecast.The improved probability verification scores due to amplified perturbation scale also supported this conclusion.

    4.Multi-scale blending experiment

    4.1.Power spectra analysis of IC perturbations derived from T639-GEPS and GRAPES-REPS

    Because of the effect of perturbation scale on ensemble forecast spread and forecastskill,asreported in section 3,itis necessary to investigate the scale characteristics of GRAPEREPS IC perturbations generated by the ETKF.The scale characteristics of T639-GEPS IC perturbations are also studied as comparison.The power spectra of IC perturbations were calculated for both T639-GEPS and GRAPE-REPS using a 2-dimemsional discrete cosine transform(2D-DCT,Denis et al.,2002),which is suitable for spectral analysis and spectral filtering of data in a limited area.

    Fora two-dimensionalfield f(i,j)of Niby Njgrid points, the direct and inverse 2D-DCT are defined as

    respectively,with

    Here,f(i,j)is the field value at grid point(i,j),and F(m,n) is the spectral coefficient corresponding to the(m,n)dimensional wave numbers.A 2D-DCT applied to a physical field f(i,j)of Niby Njvalues can produce an Niby Njarray of F(m,n)real spectral coefficients.

    In this section,the 2D-DCT is used to compute the power spectra from two-dimensional perturbation fields.In two dimensions,the perturbation state p(i,j)with Niby Njgrid points,for a variable at one level of a member,is defined as

    where a(i,j)is an ensemble member’s IC state and actl(i,j) is the control analysis state.According to Eq.(3)the twodimensional perturbation fields p(i,j)can be decomposed into a spectral field forming an Niby Njtwo-dimensional array,P(m,n).Since our goal in generating spectra is to evaluate the perturbation of two-dimensional fields as a function of spatial scales,each two-dimensional wave number pair(m,n) needs to be associated with a single-scale parameter;namely, a wavelengthλ.For a rectangular domain of Niby Nj,wehave

    where Δ is the grid point spacing.

    Figure 4a shows all-member averaged power spectra of 500 hPa temperature IC perturbations as a function of wavelength,with T639-GEPS downscaled to the same domain and resolution as GRAPES-REPS.We can see that the power of GRAPES-REPS perturbations is greater than that of T639-GEPS perturbations at wavelengths less than 1100 km.In particular,for scales less than 60 km(around two grid lengths of T639 global model),there is no power for T639-GEPS perturbations,and this is mainly because these scales cannot be resolved by the global model.Whereas,for scales over 1100 km,more power can be found in the global ensemble, as the maximum powercan reach 40 K2(corresponding to the wavelength of5000 km),while the maximum value forthe regional ensemble can only reach 20 K2(corresponding to the wavelength of 3200 km).According to the analysis in section 3,the more power at large scales may lead to a larger spread of the global ensemble,and thus if the large-scale component of global ensemble perturbations is introduced into the regional ensemble,we believe that this may be beneficial to the spread and forecast skill of the REPS.

    4.2.Construction methodology for the multi-scale blending IC perturbations

    Based on the above analysis,it is highly desirable to obtain the IC perturbations that contain the optimal components of the regional ensemble and global ensemble.Thus,we explore an MSB perturbation approach.The MSB-generated IC perturbations can be partitioned into two parts:a small-scale component and a large-scale component,in which the terms“small”and“l(fā)arge”are reflective of the relative relationship between the two components,and are not defined by exact values.The small-scale component is provided by the regional ensemble IC perturbations computed from the ETKF, while the large-scale component is provided by the global ensemble IC perturbations.A digital filter was used to conduct the scale selection.

    4.2.1.Introduction to the filter

    A low-pass filter for separating horizontal meteorological fields into different scales can be easily designed based on the aforementioned 2D-DCT(Denis et al.,2002).The direct application of the 2D-DCT to a physical field produces an array of spectral coefficients in which the spatial scales are related to the two-dimensional wave numbers(or wavelength).Lowpass filtering can easily be performed by applying a transfer function onto the two-dimensional spectral coefficients. This is done by multiplying the spectral coefficient array by a transfer function array with values between 0 and 1.In this study,a gradually varying transfer function with a soft cutoff was used to avoid Gibbs phenomena(Wang et al.,2014).The transfer function used for the cutoff transition zone follows a squared cosine.Figure 5 shows an example of such a transfer function,also commonly called the amplitude response function of the filter.For this low-pass filter,all scales shorter than“w1”km were removed and all scales larger than“w2”km were preserved.Thereafter,an inverse transform[Eq.(4)] was applied to rebuild the filtered physical field.

    4.2.2.Filter scheme for IC perturbations of T639-GEPS and GRAPES-REPS

    To conduct the scale filter for T639-GEPS and GRAPESREPS,we needed to choose appropriate critical wavelengths,“w1”and“w2”,for the transition zone.We first divided thescales that can be represented by both the GEPS and REPS into severalintervals:60–120 km,120–240 km,240–480 km, 480–960 km,960–1920 km,and 1920–3840 km,with the smallest wavelength(60 km)able to be represented by T639-GEPS and the longest wavelength(3840 km)by GRAPESREPS(7530 km in the latitudinal direction and 4950 km in the longitudinal direction).By comparing the one-month averaged power-scale curves(similar to Fig.4a)of the global ensemble and the regional ensemble,we can identify the intersection scale(at scales larger than the intersection scale, the global ensemble contains more power,while at scales smaller than the intersection scale,the regional ensemble contains more power).The intersection scales for different variables and different levels are listed in Table 2.In order to obtain the optimal combination that contains the relatively more powerful part of the two ensembles,the critical wavelengths(scales)for different variables and levels were determined according to the corresponding intersection scales listed in Table 2.Hence,for a particular variable at a given level,the transition zone was determined by the aforementioned scale interval within which the intersection scale lays.

    After the transition zones for all variables at all levels were determined,the 2D-DCT low-pass filter was performed directly to the downscaled T639-GEPS IC perturbations.This enabled the scales larger than w2 to remain untouched,while the scalessmallerthan w1 to be removed completely.

    In order to obtain the small-scale component from the regional ensemble IC perturbations,we also applied the 2DDCT low-pass filter to the GRAPES-REPS IC perturbations, with the critical wavelengths the same as those used in the T639-GEPS filtering.We then subtracted the filtered largescale component from the original IC perturbation fields to obtain the small-scale perturbation components.

    4.2.3.Method for blending the filtered perturbations

    The large-scale component of the T639-GEPS IC perturbations(denoted as T639-LS)and the small-scale component of the GRAPES-REPS IC perturbations(denoted as ETKFSS,because the IC generation method is the ETKF)were linearly combined with equal weight,such that the MSB IC perturbations were obtained by

    Table 2.Intersection scales of the IC perturbation power spectra for GRAPES-REPS and T639-GEPS.

    where IPMSBdenotes the MSB IC perturbations.

    Figure 6 shows 850 hPa u wind perturbation for T639-LS, ETKF-SS and MSB.It is clear that,after the filtering process,only the large-scale component remains in the T639-LS perturbation states,while the small-scale component remains in the ETKF-SS perturbation states.As a combination of the two,the MSB perturbation states include information at all scales.

    4.3.Experiment results

    The MSB scheme was tested in the operational GRAPESREPS environment and compared with the operational ETKF scheme.The experimental model configurations for both schemes were the same,and the experiment was carried out over a period of one month(August 2012).

    4.3.1.Spectral analysis of MSB perturbations

    The power spectrum of the MSB perturbations is evaluated.From Fig.4b,we can see that the MSB perturbations of 500 hPa temperature present significant power not only at large scales but also at small scales.Compared with Fig.4a, the MSB perturbations contain the same amount of power as GRAPES-REPS at small scales,while exhibiting the same characteristics as T639-GEPS at large scales,with the maximum power value of 43 K2.

    4.3.2.Comparison of MSB perturbations with ETKF perturbations

    The operational test results of GRAPES-REPS posed an issue that the ETKF generated perturbations exhibited slow growth.Indeed,the problem of“under-dispersion”of an ensemble system is often encountered in today’s REPS studies. This section presents the characteristics of the MSB perturbations and perturbation evolution.Perturbations derived from the ETKF are also presented as a comparison.

    Figure 7 shows the energy norm profile at all levels for 0–36 h forecast lead times.For the initial time(00 h),the ETKF ensemble perturbation energy norm has two maximum values:one located at the upper levels of 150 hPa,with a value of 2.2 J kg-1,and the other located at the bottom level, with a value of 2.17 J kg-1.We speculate that this maximum value at the bottom mainly lies underground,and that these are spurious perturbations thatcome from interpolation,especially in the plateau area.Thus,these perturbations attenuate rapidly within shortforecastlead times,and then startto grow gradually.For the MSB,the 00 h energy norm distribution at the upper levels is different from that of the ETKF,with a maximum value of 3 J kg-1at 250 hPa.This is mainly due to the introduction of the large-scale component from GEPS. For levels above 850 hPa,the energy norm is significantly larger than that of the ETKF ensemble,while for levels lower than 850 hPa,the MSB energy norm profile is similar to that of the ETKF.

    We also note that the most remarkable growth level of the ETKF perturbations is between 200 and 300 hPa,as the 36 h energy norm at 250 hPa can reach 4 J kg-1.For MSB,the most obvious energy norm growth can also be found at 250 hPa,where the 36 h energy norm is 4.4 J kg-1,which is more than that of the ETKF perturbations.The growth of MSB perturbations at other levels can also exceed that of the ETKF perturbations.Since a rapid amplification of analysis error can lead to large forecast errors,it is desirable for an ensemble to contain perturbations representative of likely analysis errors that can grow quickly(Wang and Bishop,2003).For the ETKF perturbations in Fig.7a,the maximum energy norm level does not correspond well to the fastest growth level.In other words,the most powerful perturbation within the ETKF ensemble subspace does not point to the most rapid growth direction,whereas the MSB method can adjust the energy norm distribution through adding more perturbations to the levels where energy can grow quickly,and thus higher growth rates can be found within the MSB ensemble perturbation subspace.

    4.3.3.Statistical verification

    An attempt was made to verify the performance of the MSB-based REPS and ETKF-based REPS using various score measures,including the RMSE of ensemble mean,ensemble spread,continuous rank probability score(CRPS), and percentage of outliers.

    Figures 8a–c show the one month averaged RMSE of the ensemble mean and the ensemble spread of 500 hPa temperature(T500),500 hPa wind speed(WS500)and 10m u wind (U10m),for ETKF and MSB,respectively.For all variables presented,the MSB shows an improved RMSE–spread relationship compared to that of the ETKF,given that the RMSE is smaller while the spread is larger for the MSB ensemble. Similar results can also be observed for other variables at different levels,indicating that the large-scale component fromthe GEPS can benefit the growth of ensemble spread and improve the ensemble mean forecast quality.

    Figures 8d–f show the CRPSs of T500,W500 and U10m, respectively.The CRPS has a negative orientation,the smaller the better(Hersbach,2000).It is clear that MSB performs better than the ETKF,and such an advantage of MSB over the ETKF increases with forecast lead time.CRPS verification on other variables produced similar results(not shown).

    Another measure of statistical reliability is the percentage of outliers.This is a statistic of the number of cases when the verifying analysis at any grid point lies outside the whole ensemble.A more reliable ensemble should have a smaller percentage of outliers(Wang et al.,2011).It is evident that MSB has fewer outliers than ETKF for all the variables(Figs. 8g–i),indicating a greater tendency for the observation to lie inside the MSB ensemble.

    From the above analyses,we find an overall improvement of MSB over ETKF,either for upper-air or near-surface variables.While the improvement of MSB compared to the ETKF is not overly substantial,in the context that enhancing the performance of a REPS through the initial perturbations is recognized a challenging issue,even a slight improvement is desirable.Therefore,these one-month averaged results do indeed indicate the superiority of the MSB method.

    4.3.4.A typical heavy rainfall case

    A typical heavy rainfall case in summer 2012 was studied.Both ensembles were initiated at 1200 UTC 19 August 2012,with all model configurations the same as those in the operational setting.

    Figure 9 shows the observation,the heavy precipitation probability and the spread for 24 h accumulated precipitation from 0000 UTC 20 to 0000 UTC 21 August 2012.As shown in the observation(Fig.9a),this case was characterized by a large precipitation area over central China,with the rainfallband exhibiting a west–east pattern.Figures 9b and c present the probability of24 h accumulated precipitation over50 mm. All of the probability magnitudes are represented by shading,with magnitudes over 40%highlighted by contours.For the ETKF ensemble(Fig.9b),the area with magnitudes of greater than 25%covers the northwest part of Shandong,and the high probability region(contoured)exhibits a northeast shift relative to the observed precipitation center.Also apparent is an area with probability greater than 25%located over the north coastal region of Shandong,whereas the observed precipitation was less than 5 mm in this region.Furthermore, an area with probability less than 25%can be seen over the center of Henan,but the observation shows there was a precipitation center at the same location,with a maximum valueover 90 mm.For MSB(Fig.9c),the locations and ranges of high probability areas are closer to the observation.For example,the range of probability over 40%near the border of Shandong and Henan is enlarged,and an area with probability over 25%emerged at the center of Henan,corresponding well with the observation.

    The spread of the ETKF ensemble(Fig.9d)places great emphasis on the northern part of Shandong,where the observed rain band was not located,and thus such a large precipitation spread atthatlocation ismeaningless.Whereas,the spread of the MSB ensemble(Fig.9e)has large values over the west of Shandong,which is close to the observed rainfall band.This larger spread with better locations effectively increases the chances of MSB ensemble members to capture the location of heavy rainfall accurately.

    The performances of the two ensembles with respect to their precipitation forecasts were examined by computing the one month averaged Brier Score(BS)and Area of Relative Operating Characteristic(AROC),which are appropriate for measuring the probability forecast skill of an ensemble in terms of the quantitative precipitation forecast.The BS measures the mean squared difference between the predicted probability and the observed occurrence of an event,producing a value of between zero and one,with a smaller value indicating better performance.The AROC denotes the relative relationship between the hit rate and the false alarm rate for a threshold of an event.A higher AROC means a higher hit rate and better probability forecast,and vice versa.Figure 10a presents the BS for both the ETKF and MSB,with a forecast lead time of 36 h.The MSB ensemble shows better performance than the ETKF ensemble at all precipitation magnitudes.For heavy rain greater than 50 mm,the BS of the MSB ensemble(0.013)exhibits 50%improvement over that of the ETKF(0.028).Figure 10b presents the AROC for both ensembles with a forecast lead time of 36 h.For this AROC index,the MSB ensemble has a higher score not only for light and moderate rain but also for intense rain.Therefore,this result supports the conclusion obtained from the BS analyses.

    Since the spatialdistribution ofheavy precipitation is sensitive to the atmospheric flow,more information on MSB perturbations at larger scales can help to improve the representation of atmospheric circulation factors in the model, thereby improving the spatial distribution forecast of heavy rainfall.This case study clearly demonstrates that combining the small-scale component of the regional ensemble and the large-scale component of the global ensemble is an effective way to improve the quality of probabilistic precipitation forecasts.

    5.Conclusions and discussion

    To improve the operational regional ensemble prediction system GRAPES-REPS,three issues were investigated in this study:(1)whether the scale of perturbations is an important factor affecting the ensemble spread and ensemble forecast skill;(2)what the scale characteristics of GRAPES-REPS are;and(3)whether the REPS’s skill can be improved by adding more large-scale information to the IC perturbations. Experiments were conducted to test the impact of perturbation scale on ensemble spread growth and forecast skill; the scale characteristics of GRAPES-REPS were investigated and compared with that of T639-GEPS;and an MSB IC perturbation scheme was tested.The key findings of the study can be summarized as follows:

    A“scale selective perturbation generator”was designed to generate IC perturbations characterized by a particular scale.Three sets of experiments were carried out with different perturbation scales,and the results showed that the spread growth rate of GRAPES-REPS is sensitive to the scale of IC perturbations.Overall,perturbations with larger scale exhibit a higher rate of spread growth and better skill.

    The power spectra of IC perturbations from T639-GEPS and GRAPES-REPS were analyzed and compared using a 2D-DCT.Results indicated that perturbations derived from the global ensemble show more power at larger scales,while the regional ensemble perturbations contain more power at smaller scales.The fact that the ETKF generated regional ensemble perturbations lacking large-scale information led us to consider if adding more large-scale information from a GEPS would improve the REPS’s performance.

    An MSB IC perturbation scheme was developed.This method takes advantage of both the large-scale component of global ensemble perturbations and the small-scale compo-nent of ETKF-generated IC perturbations,with the two components obtained using a 2D-DCT filter and linearly combined with equal weight.The MSB and ETKF schemes were compared through a series of experiments in the operational environment.Results showed that,compared to the ETKF, the MSB scheme can generate IC perturbations with more power at larger scales,and the power at smaller scales is also preserved.The total energy norm of MSB perturbations can maintain appropriate growth at all forecast lead times.Verification showed a higher skill of the MSB ensemble than the ETKF ensemble.A comparison between MSB and ETKF in terms of precipitation forecasting also showed better performance of the MSB ensemble.All of the results indicated that introducing a large-scale component into the regional ensemble IC perturbations can enhance the REPS’s spread and forecast skill,basically because it is this introduced large-scale information that the regional ensemble otherwise lacks.The findings of this study are useful for directing future upgrades of the current REPS.

    The effect of the MSB method on the mismatch in the LBCs of GRAPES-REPS was not investigated in the present study because the domain for GRAPES-REPS covers a wide area including the whole of China,and so the spurious gravity waves caused by the mismatch in the LBCs take little account of the IC perturbations within the whole domain.Further experiments that take account of the mismatch in the LBCs in a small-domain ensemble will be carried out in the future.

    As the MSB method is better at representing large-scale uncertainties relative to the regional ETKF,we speculate that the advantages of MSB over ETKF might be greater in cases with stronger synoptic forcing,but less so in cases with weaker synoptic forcing.We plan to investigate more cases in the future to verify this hypothesis.

    Acknowledgements.This study was supported by the National Natural Science Foundation of China(Grant No.91437113), the Special Fund for Meteorological Scienti fic Research in the Public Interest(Grant Nos.GYHY201506007 and GYHY201006015), the National 973 Program of China(Grant Nos.2012CB417204 and 2012CB955200)and the Scienti fic Research&Innovation Projects for Academic Degree Students of Ordinary Universities of Jiangsu (Grant No.KYLX 0827).The authors would like to thank Dr.DU Jun at the National Centers for Environmental Prediction and Dr. WANG Yong at the Central Institute for Meteorology and Geodynamics of Austria for their constructive comments and suggestions.

    REFERENCES

    Bowler,N.E.,and K.R.Mylne,2009:Ensemble transform Kalman filter perturbations for a regional ensemble prediction system.Quart.J.Roy.Meteor.Soc.,135,757–766.

    Bowler,N.E.,A.Arribas,K.R.Mylne,K.B.Robertson and S.E. Beare,2008:The MOGREPS short-range ensemble prediction system.Quart.J.Roy.Meteor.Soc.,134,703–722.

    Bowler,N.E.,A.Arribas,S.E.Beare,K.R.Mylne,and G.J. Shutts,2009:The local ETKF and SKEB:Upgrades to the MOGREPS short-range ensemble prediction system.Quart. J.Roy.Meteor.Soc.,135,767–776.

    Buizza,R.,and T.N.Palmer,1995:The singular-vector structure of the atmospheric global circulation.J.Atmos.Sci.,52, 1434–1456.

    Buizza,R.,M.Milleer and T.N.Palmer,1999:Stochastic representation of model uncertainties in the ECMWF ensemble prediction system.Quart.J.Roy.Meteor.Soc.,125,2887–2908.

    Buizza,R.,P.L.Houtekamer,Z.Toth,P.Pellerin,M.Z.Wei,and Y.J.Zhu,2005:A comparison of the ECMWF,MSC,and NCEP global ensemble prediction systems.Mon.Wea.Rev.,133,1076–1097.

    Caron,J.F.,2013:Mismatching perturbations at the lateral boundaries in limited-area ensemble forecasting:A case study. Mon.Wea.Rev.,141,356–374.

    Chen,D.H.,and X.S.Shen,2006:Recent progress on GRAPES research and application.Journal of Applied Meteorological Science,17,773–777.(in Chinese)

    Chen,J.,and J.S.Xue,2009:Heavy rainfall ensemble prediction: Initial condition perturbation vs multi-physics perturbation. Acta Meteorologica Sinica,23(1),53–67.

    Chen,J.,J.Xue,H.Yan,2004:The impact of non-adiabatic physics on mesoscale short-range heavy rainfall prediction. Acta Meteorological Sinica,18(1),51–72.

    Chen,J.,J.Xue,H.Yan,2005a:The uncertainty of heavy rainfall of South China and experiments of ensemble prediction.Acta Meteorological Sinica,19(1),1–18.

    Chen,J.,J.S.Xue,and H.Yan,2005b:A new initial perturbation method of ensemble mesoscale heavy rain prediction.Chinese Journal of Atmospheric Sciences,29,717–726.(in Chinese)

    Denis,B.,J.C?ot′e,and R.Laprise,2002:Spectraldecomposition of two-dimensional atmospheric fields on limited-area domains using the discrete cosine transform(DCT).Mon.Wea.Rev.,130,1812–1829.

    Du,J.,and M.S.Tracton,2001:Implementation of a real-time short range ensemble forecasting system at NCEP:An update.Preprints,9th Conference on Mesoscale Processes,Ft. Lauderdale,Florida,Amer.Meteor.Soc.,355–356.

    Du,J.,G.DiMego,M.S.Tracton,and B.Zhou,2003:NCEP short-range ensemble forecasting(SREF)system:Multi-IC, multi-model and multi-physics approach.Research Activities in Atmospheric and Oceanic Modeling,J.C?ot′e,Ed.,Rep. 33,CAS/JSC Working Group Numerical Experimentation (WGNE),WMO/TD-No.1161,5.09–5.10.

    Du,J.,G.DiMego,B.B.Zhou,D.Jovic,B.Ferrier,B.Yang, S.Benjamin,2014:NCEP regional ensembles:Evolving toward hourly-updated convection-allowing scale and stormscale predictions within a unified regional modeling system. 26th Conf.on Weather Analysis and Forecasting and 22nd Conf.on Numerical Weather Prediction,Atlanta,GA,Amer. Meteor.Soc.,1–6 Feb.2014,paper J1.4.

    Duan,Y.H.,and Coauthors,2012:An overview of the Beijing 2008 Olympics Research and Development Project (B08RDP).Bull.Amer.Meteor.Soc.,93,381–403.

    Frogner,I.-L.,H.Haakenstad,T.Iversen,2006:Limited-area ensemble predictions atthe Norwegian MeteorologicalInstitute. Quart.J.Roy.Meteor.Soc.,132,2785–2808.

    Grimit,E.P.,and C.F.Mass,2002:Initial results of a mesoscale short-range ensemble forecasting system over the Pacific Northwest.Wea.Forecasting,17,192–205.

    Hersbach,H.,2000:Decomposition of the continuous rankedprobability score for ensemble prediction systems.Wea.Forecasting,15,559–570.

    Houtekamer,P.L.,L.Lefaivre,J.Derome,H.Ritchie,and H.L. Mitchell,1996:A system simulation approach to ensemble prediction.Mon.Wea.Rev.,124,1225–1242.

    Lacarra,J.F.,and O.Talagrand,1988:Short-range evolution of small perturbations in a barotropic model.Tellus A,40A,81–95.

    Leith,C.E.,1974:Theoretical skill of Monte Carlo forecasts. Mon.Wea.Rev.,102,409–418.

    Li,X.L.,M.Charron,L.Spacek,and G.Candille,2008:A regional ensemble prediction system based on moist targeted singular vectors and stochastic parameter perturbations.Mon. Wea.Rev.,136,443–462.

    Lorenz,E.N.,1969:Predictability of a flow which possesses many scales of motion.Tellus,21,289–307.

    Ma,X.L.,J.S.Xue,and W.S.Lu,2008:Preliminary study on ensemble transform Kalman filter-based initial perturbation scheme in GRAPES global ensemble prediction.Acta Meteorologica Sinica,66(4),526–536.(in Chinese)

    Marsigli,C.,F.Boccanera,A.Montani,and T.Paccagnella,2005: The COSMO-LEPS mesoscale ensemble system:Validation of the methodology and verification.Nonlinear Processes in Geophysics,12,527–536.

    Molteni,F.,R.Buizza,T.N.Palmer,and T.Petroliagis,1996:The ECMWF ensemble prediction system:Methodology and validation.Quart.J.Roy.Meteor.Soc.,122,73–119.

    Palmer,T.N.,R.Gelaro,J.Barkmeijer and R.Buizza,1998:Singular vectors,metrics,and adaptive observations.J.Atmos. Sci.,55,633–653.

    Saito,K.,M.Hara,H.Seko,M.Kunii,and M.Yamaguchi,2011: Comparison of initial perturbation methods for the mesoscale ensemble prediction system of the Meteorological Research Institute for the WWRP Beijing 2008 Olympics Research and Development Project(B08RDP).Tellus A,63,445–467.

    Stensrud,D.J.,and N.Yussouf,2007:Reliable probabilistic quantitative precipitation forecasts from a short-range ensemble forecasting system.Wea.Forecasting,22,3–17.

    Stensrud,D.J.,H.E.Brooks,J.Du,M.S.Tracton,and E.Rogers, 1999:Using ensembles for short-range forecasting.Mon. Wea.Rev.,127,433–446.

    Toth,Z.,and E.Kalnay,1993:Ensemble forecasting at NMC: The generation of perturbations.Bull.Amer.Meteor.Soc.,74, 2317–2330.

    Toth,Z.,and E.Kalnay,1997:Ensemble forecasting at NCEP and the breeding method.Mon.Wea.Rev.,125,3297–3319.

    Wang,X.G.,and C.H.Bishop,2003:A comparison of breeding and ensemble transform Kalman filter ensemble forecast schemes.J.Atmos.Sci.,60,1140–1158.

    Wang,X.G.,C.H.Bishop,and S.J.Julier,2004:Which is better, an ensemble of positive–negative pairs or a centered spherical simplex ensemble?Mon.Wea.Rev.,132,1590–1605.

    Wang,Y.,M.Bellus,J.F.Geleyn,X.L.Ma,W.H.Tian,and F.Weidle,2014:A new method for generating initial condition perturbations in a regional ensemble prediction system: Blending.Mon.Wea.Rev.,142,2043–2059.

    Wang,Y.,and Coauthors,2011:The Central European limitedarea ensemble forecasting system:ALADIN-LAEF.Quart. J.Roy.Meteor.Soc.,137,483–502.

    Wei,M.Z.,Z.Toth,R.Wobus,and Y.J.Zhu,2008:Initial perturbations based on the ensemble transform(ET)technique in the NCEP global operational forecast system.Tellus A,60, 62–79.

    Wei,M.A.,Z.Toth,R.Wobus,Y.J.Zhu,C.H.Bishop and X.G. Wang,2006:Ensemble Transform Kalman Filter-based ensemble perturbations in an operational global prediction system at NCEP.Tellus A,58,28–44.

    Zhang,H.B.,J.Chen,X.F.Zhi,K.J.Long and Y.N.Wang, 2014a:Design and comparison of perturbation schemes for GRAPES Meso based ensemble forecast.Transactions of Atmospheric Sciences,37,276–284.(in Chinese)

    Zhang,H.B.,J.Chen,X.F.Zhi,Y.L.Li and Y.Sun,2014b: Study on the application of GRAPES regional ensemble prediction system.Meteorological Monthly,40,1076–1087.(in Chinese)

    :Zhang,H.B.,J.Chen,X.F.Zhi,Y.Wang,and Y.N.Wang,2015:Study on multi-scale blending initial condition perturbationsfor a regionalensemble prediction system.Adv.Atmos.Sci.,32(8),1143–1155,

    10.1007/s00376-015-4232-6.

    27 October 2014;revised 16 January 2015;accepted 22 January 2015)

    ?Corresponding author:CHEN Jing

    Email:chenj@cma.gov.cn

    最新中文字幕久久久久| 天美传媒精品一区二区| 免费av观看视频| 婷婷色av中文字幕| 亚洲婷婷狠狠爱综合网| 欧美不卡视频在线免费观看| 欧美3d第一页| 日韩 亚洲 欧美在线| 成年免费大片在线观看| 成人亚洲精品av一区二区| 亚洲美女搞黄在线观看| 欧美最黄视频在线播放免费| 成人一区二区视频在线观看| 伊人久久精品亚洲午夜| 一级黄色大片毛片| 久久久久久久久久久免费av| 日本成人三级电影网站| 不卡视频在线观看欧美| 91精品一卡2卡3卡4卡| 国产乱人偷精品视频| 免费黄网站久久成人精品| 十八禁国产超污无遮挡网站| av在线老鸭窝| 少妇高潮的动态图| 亚洲人成网站在线播放欧美日韩| 欧美极品一区二区三区四区| 99久国产av精品| 女人被狂操c到高潮| 在线播放国产精品三级| 人妻久久中文字幕网| 国产av一区在线观看免费| 国产单亲对白刺激| 亚洲国产精品国产精品| 91狼人影院| 看黄色毛片网站| 99久久成人亚洲精品观看| 舔av片在线| 两个人视频免费观看高清| 18禁裸乳无遮挡免费网站照片| 欧美精品国产亚洲| 少妇的逼水好多| 国产av麻豆久久久久久久| 色5月婷婷丁香| 久久久精品大字幕| 嘟嘟电影网在线观看| 能在线免费观看的黄片| 少妇熟女欧美另类| 一进一出抽搐动态| 久久国内精品自在自线图片| 国模一区二区三区四区视频| 性色avwww在线观看| 毛片一级片免费看久久久久| 毛片一级片免费看久久久久| 亚洲高清免费不卡视频| 一夜夜www| 国产成人福利小说| 色播亚洲综合网| 99riav亚洲国产免费| 日韩欧美精品免费久久| 亚洲四区av| 久久久国产成人免费| 三级毛片av免费| 国产黄色视频一区二区在线观看 | 久久精品国产自在天天线| 日韩av不卡免费在线播放| 啦啦啦韩国在线观看视频| 十八禁国产超污无遮挡网站| 韩国av在线不卡| 国产高潮美女av| 黄色一级大片看看| 亚洲内射少妇av| av免费在线看不卡| 一区福利在线观看| 嫩草影院精品99| 少妇裸体淫交视频免费看高清| 色播亚洲综合网| 给我免费播放毛片高清在线观看| 一本一本综合久久| 高清毛片免费观看视频网站| 毛片一级片免费看久久久久| 亚洲av第一区精品v没综合| 免费看美女性在线毛片视频| 午夜老司机福利剧场| 国产黄色视频一区二区在线观看 | 国产精品av视频在线免费观看| 夜夜看夜夜爽夜夜摸| 少妇的逼水好多| 两性午夜刺激爽爽歪歪视频在线观看| 亚洲av电影不卡..在线观看| 中国国产av一级| а√天堂www在线а√下载| 亚洲精品国产成人久久av| 亚洲国产日韩欧美精品在线观看| 特级一级黄色大片| 欧美色视频一区免费| 男人的好看免费观看在线视频| 亚洲欧美清纯卡通| 在线免费观看不下载黄p国产| 99精品在免费线老司机午夜| 亚洲人成网站在线播放欧美日韩| 国产亚洲5aaaaa淫片| 长腿黑丝高跟| 97超视频在线观看视频| 狠狠狠狠99中文字幕| 免费人成视频x8x8入口观看| 国产精华一区二区三区| 欧美成人a在线观看| 日韩亚洲欧美综合| 精品日产1卡2卡| 黄色欧美视频在线观看| 亚洲在线观看片| 中国美白少妇内射xxxbb| 日韩,欧美,国产一区二区三区 | 熟女电影av网| 热99re8久久精品国产| 久久精品国产99精品国产亚洲性色| 久久这里有精品视频免费| 婷婷精品国产亚洲av| 神马国产精品三级电影在线观看| 男人舔女人下体高潮全视频| 26uuu在线亚洲综合色| 国产精品久久久久久精品电影| 中文精品一卡2卡3卡4更新| 亚洲国产欧美在线一区| 亚洲熟妇中文字幕五十中出| eeuss影院久久| 国产麻豆成人av免费视频| 午夜亚洲福利在线播放| www.色视频.com| 99热网站在线观看| 国内少妇人妻偷人精品xxx网站| 国产乱人视频| 日本-黄色视频高清免费观看| 亚洲精品久久久久久婷婷小说 | 级片在线观看| 亚洲五月天丁香| 中国美女看黄片| 日韩精品青青久久久久久| 网址你懂的国产日韩在线| 国语自产精品视频在线第100页| 人妻制服诱惑在线中文字幕| 麻豆成人av视频| 三级经典国产精品| 成人亚洲欧美一区二区av| 黄色配什么色好看| 色哟哟哟哟哟哟| 亚洲av.av天堂| 国产一区二区在线av高清观看| 国产一区二区三区在线臀色熟女| 免费搜索国产男女视频| 精品无人区乱码1区二区| 熟女人妻精品中文字幕| 国产精品伦人一区二区| 亚洲欧美成人精品一区二区| 国产精品久久久久久久电影| 在线观看一区二区三区| 国产成年人精品一区二区| а√天堂www在线а√下载| 亚洲欧美中文字幕日韩二区| 亚洲av成人精品一区久久| 天堂网av新在线| 国产精品.久久久| 九九在线视频观看精品| 99热精品在线国产| 精品不卡国产一区二区三区| 精品久久久久久久久久免费视频| 亚洲人成网站在线观看播放| 亚洲av二区三区四区| 日韩欧美精品v在线| 国产伦精品一区二区三区视频9| 热99在线观看视频| 国产黄色小视频在线观看| 一级毛片我不卡| 偷拍熟女少妇极品色| 日本一二三区视频观看| 99热这里只有是精品在线观看| 国产精品一二三区在线看| 天堂中文最新版在线下载 | 久久人人爽人人片av| 一进一出抽搐动态| 亚洲av电影不卡..在线观看| 亚洲人成网站在线观看播放| 一个人看的www免费观看视频| 麻豆久久精品国产亚洲av| 色综合亚洲欧美另类图片| 亚洲aⅴ乱码一区二区在线播放| 可以在线观看的亚洲视频| 如何舔出高潮| 赤兔流量卡办理| 亚洲经典国产精华液单| 欧美又色又爽又黄视频| 97人妻精品一区二区三区麻豆| 色噜噜av男人的天堂激情| 黄片wwwwww| 变态另类丝袜制服| 一个人观看的视频www高清免费观看| 久久99热6这里只有精品| 啦啦啦韩国在线观看视频| 成年版毛片免费区| 日韩成人伦理影院| 午夜福利在线观看免费完整高清在 | 亚洲在线观看片| 人妻制服诱惑在线中文字幕| 国产av在哪里看| 免费不卡的大黄色大毛片视频在线观看 | 99热只有精品国产| 亚洲熟妇中文字幕五十中出| 插阴视频在线观看视频| 最近2019中文字幕mv第一页| 非洲黑人性xxxx精品又粗又长| 老女人水多毛片| 婷婷亚洲欧美| 少妇丰满av| 亚洲精品国产av成人精品| 亚洲天堂国产精品一区在线| 亚洲在久久综合| 深夜a级毛片| 久久综合国产亚洲精品| 久久精品久久久久久久性| 可以在线观看毛片的网站| 身体一侧抽搐| 最近视频中文字幕2019在线8| 久久精品久久久久久久性| 91狼人影院| 欧美又色又爽又黄视频| 精品久久久久久久久av| 精品欧美国产一区二区三| 女同久久另类99精品国产91| 国产免费一级a男人的天堂| 国国产精品蜜臀av免费| 精品人妻熟女av久视频| 国产伦在线观看视频一区| avwww免费| 毛片一级片免费看久久久久| 午夜久久久久精精品| 在线a可以看的网站| 18禁黄网站禁片免费观看直播| 观看美女的网站| 人妻系列 视频| 国产麻豆成人av免费视频| 亚洲精品色激情综合| 性插视频无遮挡在线免费观看| 久久久久性生活片| 最近最新中文字幕大全电影3| 国产精品不卡视频一区二区| 久久久久久久久久成人| 身体一侧抽搐| 成人三级黄色视频| 国产欧美日韩精品一区二区| 中文字幕久久专区| 床上黄色一级片| 亚洲av免费在线观看| 亚洲国产精品成人综合色| 精品人妻熟女av久视频| 99热这里只有精品一区| 免费观看精品视频网站| 国产伦一二天堂av在线观看| 最近中文字幕高清免费大全6| 色播亚洲综合网| 老司机福利观看| 联通29元200g的流量卡| 国产伦理片在线播放av一区 | 别揉我奶头 嗯啊视频| 黄色欧美视频在线观看| 精品不卡国产一区二区三区| 亚洲精品色激情综合| 少妇人妻精品综合一区二区 | 婷婷亚洲欧美| 国产亚洲5aaaaa淫片| av在线天堂中文字幕| 欧美3d第一页| 国产精品一区二区性色av| av卡一久久| 欧美+日韩+精品| 亚洲国产欧美在线一区| 国产精品永久免费网站| 亚洲精品日韩av片在线观看| 青春草国产在线视频 | 欧美xxxx黑人xx丫x性爽| 中国美白少妇内射xxxbb| 高清日韩中文字幕在线| 亚洲av第一区精品v没综合| 午夜爱爱视频在线播放| 高清在线视频一区二区三区 | 精品久久久久久久人妻蜜臀av| av福利片在线观看| 亚洲精品成人久久久久久| 久久午夜福利片| 精品久久久噜噜| 精品一区二区三区人妻视频| 最近视频中文字幕2019在线8| 国产中年淑女户外野战色| 久久久欧美国产精品| 成年女人看的毛片在线观看| 国产精品免费一区二区三区在线| 亚洲av第一区精品v没综合| 日本撒尿小便嘘嘘汇集6| 免费av不卡在线播放| 一夜夜www| 亚洲欧美成人精品一区二区| 国产视频首页在线观看| 啦啦啦韩国在线观看视频| 亚洲成人精品中文字幕电影| 欧美成人a在线观看| 亚洲图色成人| 午夜福利在线观看免费完整高清在 | 亚洲va在线va天堂va国产| 内地一区二区视频在线| 不卡一级毛片| 国产精品爽爽va在线观看网站| 日本撒尿小便嘘嘘汇集6| 婷婷亚洲欧美| 哪里可以看免费的av片| 久久久久久久久久久丰满| 欧美+亚洲+日韩+国产| 国产精品久久视频播放| 人人妻人人看人人澡| 美女黄网站色视频| 非洲黑人性xxxx精品又粗又长| 国产精品国产高清国产av| 变态另类丝袜制服| 夜夜夜夜夜久久久久| 欧美不卡视频在线免费观看| 99热网站在线观看| 给我免费播放毛片高清在线观看| 真实男女啪啪啪动态图| 我的老师免费观看完整版| 乱人视频在线观看| 男女那种视频在线观看| 人妻夜夜爽99麻豆av| 日本爱情动作片www.在线观看| 黄色欧美视频在线观看| 看非洲黑人一级黄片| 美女 人体艺术 gogo| 99久久九九国产精品国产免费| 九九热线精品视视频播放| 国产 一区 欧美 日韩| 中文字幕av成人在线电影| 一个人观看的视频www高清免费观看| 国产在线男女| 色吧在线观看| a级毛色黄片| 久久鲁丝午夜福利片| 内射极品少妇av片p| 亚洲欧美日韩高清在线视频| 国内精品美女久久久久久| 精品久久国产蜜桃| 美女高潮的动态| 久久精品91蜜桃| 欧美性猛交黑人性爽| 春色校园在线视频观看| 哪个播放器可以免费观看大片| 日韩成人av中文字幕在线观看| 色尼玛亚洲综合影院| 在线天堂最新版资源| 久久久精品94久久精品| 高清毛片免费看| 午夜福利高清视频| 日日干狠狠操夜夜爽| 蜜桃久久精品国产亚洲av| 精品一区二区三区人妻视频| 日本av手机在线免费观看| 国产高清三级在线| 日韩精品青青久久久久久| 久久亚洲精品不卡| 波多野结衣高清无吗| 日韩成人av中文字幕在线观看| 欧美高清性xxxxhd video| 国产老妇女一区| av免费在线看不卡| 麻豆乱淫一区二区| 免费观看在线日韩| 搡老妇女老女人老熟妇| 日韩视频在线欧美| 国产午夜精品久久久久久一区二区三区| 在线观看午夜福利视频| 精品国内亚洲2022精品成人| 日本黄大片高清| 国产男人的电影天堂91| 毛片女人毛片| 久久久久久国产a免费观看| 少妇人妻精品综合一区二区 | 成年女人看的毛片在线观看| av卡一久久| 午夜免费男女啪啪视频观看| 波多野结衣高清作品| 精品日产1卡2卡| 又粗又爽又猛毛片免费看| 免费大片18禁| 亚洲在线自拍视频| 国产精品av视频在线免费观看| 三级毛片av免费| 国产女主播在线喷水免费视频网站 | 如何舔出高潮| 男的添女的下面高潮视频| 两个人视频免费观看高清| 高清在线视频一区二区三区 | 小蜜桃在线观看免费完整版高清| 免费一级毛片在线播放高清视频| 久久精品国产鲁丝片午夜精品| 你懂的网址亚洲精品在线观看 | 最好的美女福利视频网| 国产成人a∨麻豆精品| 天堂av国产一区二区熟女人妻| 国产私拍福利视频在线观看| 在线免费十八禁| 在线播放无遮挡| 日韩欧美精品v在线| 欧美激情在线99| 午夜精品一区二区三区免费看| 久久久午夜欧美精品| 午夜老司机福利剧场| 少妇熟女aⅴ在线视频| 国产av不卡久久| 又爽又黄a免费视频| 少妇被粗大猛烈的视频| 午夜激情福利司机影院| 亚洲中文字幕一区二区三区有码在线看| 两个人视频免费观看高清| 黄色配什么色好看| 午夜老司机福利剧场| 国产成人精品久久久久久| 性色avwww在线观看| 亚洲性久久影院| 亚洲av.av天堂| 成熟少妇高潮喷水视频| 人人妻人人澡欧美一区二区| 99久国产av精品国产电影| 久久久久网色| 老熟妇乱子伦视频在线观看| 狠狠狠狠99中文字幕| 国产亚洲精品av在线| 久久久久久国产a免费观看| 大又大粗又爽又黄少妇毛片口| 我的女老师完整版在线观看| 国产精品久久久久久久久免| 亚洲自拍偷在线| 色视频www国产| av天堂在线播放| 中文字幕精品亚洲无线码一区| 日韩欧美 国产精品| 一区二区三区四区激情视频 | 精华霜和精华液先用哪个| 可以在线观看的亚洲视频| 国产精品1区2区在线观看.| 亚洲av熟女| 婷婷亚洲欧美| 国产一区二区三区在线臀色熟女| 久久久久久久久中文| 亚洲欧洲日产国产| 中文字幕久久专区| 亚州av有码| 亚洲不卡免费看| 国产乱人视频| 国产成人aa在线观看| 人妻夜夜爽99麻豆av| 久久九九热精品免费| 日本黄色视频三级网站网址| 久久国内精品自在自线图片| 成人特级黄色片久久久久久久| 能在线免费看毛片的网站| 国产午夜精品一二区理论片| 深夜精品福利| 国产精品麻豆人妻色哟哟久久 | 亚洲国产欧美在线一区| 欧美色欧美亚洲另类二区| 波多野结衣高清无吗| 午夜激情欧美在线| 精品人妻偷拍中文字幕| 欧美成人一区二区免费高清观看| 久久精品国产亚洲av涩爱 | 免费观看精品视频网站| 99热只有精品国产| 老女人水多毛片| 国产精品一及| 免费黄网站久久成人精品| 美女大奶头视频| 啦啦啦啦在线视频资源| av天堂在线播放| 中文字幕精品亚洲无线码一区| 亚洲精品乱码久久久久久按摩| 免费观看的影片在线观看| 变态另类成人亚洲欧美熟女| 日韩av在线大香蕉| 国产精品不卡视频一区二区| 成人特级av手机在线观看| 欧美一级a爱片免费观看看| 日韩欧美三级三区| 国产av一区在线观看免费| 少妇的逼水好多| 久久久久久久午夜电影| 99久久精品一区二区三区| 色哟哟·www| 小蜜桃在线观看免费完整版高清| 国产一区二区在线av高清观看| 国产高清三级在线| 免费av毛片视频| 色尼玛亚洲综合影院| 日韩欧美一区二区三区在线观看| 午夜福利成人在线免费观看| 自拍偷自拍亚洲精品老妇| 51国产日韩欧美| 99热6这里只有精品| 天堂网av新在线| 国产高清激情床上av| 欧美成人免费av一区二区三区| 国产亚洲5aaaaa淫片| 热99re8久久精品国产| 午夜福利高清视频| 久久久久网色| 成年女人看的毛片在线观看| 一边摸一边抽搐一进一小说| 在线观看66精品国产| 嫩草影院入口| 少妇熟女aⅴ在线视频| 久久综合国产亚洲精品| 久久久欧美国产精品| 中文字幕熟女人妻在线| 免费人成在线观看视频色| 国产男人的电影天堂91| 青春草视频在线免费观看| 国产男人的电影天堂91| 国产精品电影一区二区三区| 色哟哟·www| 欧美高清成人免费视频www| 国产精品精品国产色婷婷| 色5月婷婷丁香| 91在线精品国自产拍蜜月| 欧美日韩综合久久久久久| 一本久久精品| 1000部很黄的大片| 能在线免费观看的黄片| 亚洲av中文av极速乱| 日本撒尿小便嘘嘘汇集6| 村上凉子中文字幕在线| 91av网一区二区| 久久久精品欧美日韩精品| 国产午夜精品论理片| 超碰av人人做人人爽久久| 国产成人一区二区在线| 久久久精品欧美日韩精品| 久久久久免费精品人妻一区二区| 一级毛片久久久久久久久女| 亚洲国产精品合色在线| 日韩欧美精品免费久久| 麻豆久久精品国产亚洲av| a级毛色黄片| 男女做爰动态图高潮gif福利片| 久久久精品大字幕| 青春草亚洲视频在线观看| 大香蕉久久网| 亚洲精品国产av成人精品| 国产高潮美女av| 国产单亲对白刺激| 五月伊人婷婷丁香| 丝袜喷水一区| 淫秽高清视频在线观看| 国产精品一区二区三区四区免费观看| 99riav亚洲国产免费| 又粗又爽又猛毛片免费看| 中出人妻视频一区二区| 久久久久久伊人网av| 99久久精品一区二区三区| 日本爱情动作片www.在线观看| 亚洲经典国产精华液单| 在线免费十八禁| 国产精品福利在线免费观看| 内地一区二区视频在线| 99热只有精品国产| 精品久久久久久久人妻蜜臀av| 一级二级三级毛片免费看| 国产精品国产高清国产av| 在线观看午夜福利视频| 亚洲成人久久性| a级毛色黄片| 色噜噜av男人的天堂激情| 国产综合懂色| 成人特级av手机在线观看| 看免费成人av毛片| 一本久久中文字幕| 国产乱人视频| 成年av动漫网址| 又粗又硬又长又爽又黄的视频 | 99久久久亚洲精品蜜臀av| 人妻系列 视频| 色噜噜av男人的天堂激情| 久久人妻av系列| 亚洲18禁久久av| 日韩一区二区视频免费看| 激情 狠狠 欧美| 国产精品国产高清国产av| 亚洲五月天丁香| av.在线天堂| 亚洲人与动物交配视频| 黄片wwwwww| 日韩大尺度精品在线看网址| 欧美日韩乱码在线| 成人亚洲欧美一区二区av| 国产人妻一区二区三区在| 村上凉子中文字幕在线| 老女人水多毛片| 国产精品精品国产色婷婷| 国产高清不卡午夜福利| 可以在线观看的亚洲视频| 亚洲成人中文字幕在线播放| 日本一本二区三区精品| 国产精品女同一区二区软件| 女同久久另类99精品国产91| 国产黄色视频一区二区在线观看 | 国产高清视频在线观看网站| 色5月婷婷丁香| 18禁裸乳无遮挡免费网站照片| 亚洲精品亚洲一区二区| 久久久久久久久久久丰满| a级毛片a级免费在线|