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

    Comparisons of Three-Dimensional Variational Data Assimilation and Model Output Statistics in Improving Atmospheric Chemistry Forecasts

    2018-05-19 06:00:38ChaoqunMATijianWANGZengliangZANGandZhijinLI
    Advances in Atmospheric Sciences 2018年7期

    Chaoqun MA,Tijian WANG?,Zengliang ZANG,and Zhijin LI

    1School of Atmospheric Sciences,Nanjing University,Nanjing 210023,China

    2CMA-NJU Joint Laboratory for Climate Prediction Studies,Nanjing 210023,China

    3Jiangsu Collaborative Innovation Center for Climate Change,Nanjing 210023,China

    4Institute of Meteorology and Oceanography,PLA University of Science and Technology,Nanjing 210000,China

    5Joint Institute for Regional Earth System Science and Engineering,University of California,Los Angeles 90001,California,USA

    1.Introduction

    In recent years an unexpected outbreak of severe air pollution events has engulfed China during the autumn and winter months of the year.These air pollution episodes have aroused deep concern and panic amongst the public and have subsequently become a top priority to address for the government.For example,the State Council released the Air Pollution Prevention and Control Action Plan in September 2013,aimed at reducing particulate pollution.In addition,targets to control sulfur dioxide and nitrogen oxides were listed in the 11th and 12th Five-Year Plan drawn up by the National Development and Reform Commission.However,to control atmospheric pollution efficiently,it is necessary to achieve more accurate forecasting of atmospheric chemical constituents.

    Since the beginning of the 21st century,air quality forecast systems,such as the Weather Research and Forecasting-Chemistry(WRF-Chem)model,have been gradually put into operation in key cities across China by many organizations and institutions.However,without any additional measures,these numerical forecast systems are not accurate enough to be applied in operational air quality forecasts due to the uncertainties within their parameterization schemes and input data(van Loon et al.,2007;Zhang et al.,2016).Accordingly,scientists have developed multiple pre-or post-processing techniques for model improvement in the operational prediction of meteorological and chemical fields.For example,data assimilation(DA)—a measure applied before the model run—is an effective approach in improving the model forecast skill of air pollution via reducing the uncertainty of chemical initial conditions(CICs)or other parameters.For instance,Barbu et al.(2009)achieved a better forecast by assimilating measurements of sulfur dioxide(SO2)and sulfate to adjust the emission and conversion rates of SO2in the model;and the research of Liu et al.(2011)and Yin et al.(2016)showed improved aerosol analysis and forecasting by assimilating the MODIS total aerosol optical depth retrieval products.Different models and observational data have been tested and their results show similar conclusions(Wang et al.,2014;Zhang et al.,2015;Mizzi et al.,2016;Tang et al.,2016)—readers who are interested should refer to Bocquet et al.(2015)for more details.Alternatively,the forecast error can also be corrected effectively by another approach,called model output statistics(MOS)(Glahn and Lowry,1972).MOS works through statistically relating the historical model output with corresponding observations,and then applying the relationship to the model forecast.This approach has been widely used in the post-processing of operational numerical weather prediction,and in doing so most forecast biases can be corrected(Wilson and Vall′ee,2003),especially for temperature(Taylor and Leslie,2005;Libonati et al.,2008)and humidity(Anadranistakis et al.,2004).

    Although it has been proven that both DA and MOS are effective in improving the forecast performance,little attention has been paid to comparing or combining the two methods,especially with respect to atmospheric chemistry modeling.There are two obstacles that make it difficult to compare the two methods fairly using the same observational dataset.First,most early(and even recent)research concerning 3DVar in atmospheric chemistry models has focused on assimilating observations from satellite-derived products to generate analyses that are skillful in improving the forecasts of variables like carbon monoxide(Barret et al.,2008),carbon dioxide,ozone(O3),nitrogen dioxide(NO2)(Inness et al.,2015;Wargan et al.,2015),methane(Alexe et al.,2015),and aerosols(Benedetti et al.,2009;Yumimoto et al.,2016).However,MOS works only with in-situ observations from surface stations.Second,unlike numerical weather prediction,MOS still remains in its infancy in terms of the operational numerical forecasting of atmospheric chemistry variables.Such works,if any,are usually based on regression approaches(Denby et al.,2008;Honore et al.,2008;Struzewska et al.,2016),which are effective in improving the air quality forecast for all analyzed species.However,methods based on these approaches are usually too dependent on local pollution conditions,which makes them inconvenient to be applied as widely as DA approaches.

    This study compares the potential of the two approaches in improving the atmospheric chemistry forecasts of the WRF-Chem modeling system in an operational context.To overcome the problems mentioned above,we firstly adapted a 3DVar DA system based on Li et al.(2013)(L13 henceforth)and Jiang et al.(2013)to the assimilation of observational data from surface stations.Then,we modified a MOS scheme from Galanis and Ana dranistakis(2002)(hereafter G02)for use in adjusting the meteorological forecast,to enable it to correct the chemistry output from the atmospheric model.Following this introduction,the model,setup and experimental design are described in section 2.Section 3 evaluates the model improvement with the two methods.A summary and conclusion are given in section 4.

    2.Method and data

    2.1.Model setup

    Fig.1.The nested domains of the WRF-Chem model:D01(blue frame)is the outer domain,which covers most of East Asia;D02(black frame)is the inner domain,in which Hebei Province(red boundary)is located.

    The WRF-Chem model is an online 3D,Eulerian chemical transport model that considers the complex physical and chemical processes in the troposphere(Grell et al.,2005).It has been applied in various research settings,especially those concerning feedbacks of air pollution to meteorological and chemical DA(Saide et al.,2012,2015;Makar et al.,2015;Mizzi et al.,2016).In this study,version 3.7.1 of WRF-Chem was used to simulate the air quality in Hebei Province,China.Two nested domains were set,as shown in Fig.1.The outer domain covered East Asia with a hori-zontal resolution of 75 km×75 km and 106(lon)×81(lat)grids,while the inner domain covered Hebei Province with a horizontal resolution of 15 km×15 km and 76(lon)×81(lat)grids.The vertical resolution of the model was 24 vertical levels(about 6 levels below 1 km and 20 levels below 10 km),with 100 hPa as the model top.The 0.5°×0.5°data from the NCEP’s Global Forecast System were used to provide the meteorological initial conditions and lateral boundary meteorological conditions every 12 hours.Atmospheric gaseous chemistry and aerosols were simulated using the Second Generation Regional Acid Deposition Model(RADM2)with the Modal Aerosol Dynamics Model for Europe(MADE)and the Secondary Organic Aerosol Model(SORGAM)(Stock well et al.,1990;Ackermann et al.,1998;Schell et al.,2001)scheme.The anthropogenic emissions inventory was the Multi-resolution Emissions Inventory for China in 2012(http://www.meicmodel.org).The con figuration of WRF-Chem is detailed in Table 1.

    2.2.MOS method

    What MOS does is to find statistic relationships from training samples that can then be applied in model forecast outputs.By doing so,the expectation is that model errors will be corrected and a forecast generated that will better fit the observation.In this study,a one-dimensional Kalman filter was chosen as the algorithm to realize the MOS process.The algorithm was formulated in a way that generally resembled G02;the only modi fications were as follows:

    Firstly,in this work,the measurementy(t),as well as the real valuex(t),could both be the difference and ratio between the forecast and observation,whereas in Eqs.(1)and(2)of G02 they only denoted the difference.Furthermore,hourly concentrations of five species from the three-day model output were split into 3×5×24 independent daily concentration series.Lastly,given Kalman filtering can only predict one time-step ahead(one day ahead,in this context),the correction could only work during forecast hours 0–24(24-h forecast hereafter),while leaving forecast hours 24–48 and 48–72(48-h and 72-h forecast hereafter)uncorrected.Therefore,to extend the algorithm further,the corrected results from the 24-h(48-h)forecast was used as a proxy or substitute observation at the corresponding time to correct the 48-h(72-h)model output.Appendix A describes the steps in more detail.

    Table 1.Con figuration of the physical and chemical schemes of WRF-Chem.

    2.3.DA con figuration

    In this study,3DVar DA was implemented to optimize the CICs for the inner model domain.The DA system and formulation used were based on L13 with the following modi fications:

    In addition to the fine particulate matter(PM2.5)assimilated in L13,particulate matter with diameters between 2.5μm and 10 μm(PM2.5–PM10)was also assimilated and the analysis increment was added to the corresponding model variables following L13.Gaseous species,including SO2,NO2and O3,were also assimilated to decrease the uncertainty of their concentrations in the model CICs.

    Following L13,the National Meteorological Center(NMC)method(Parrish and Derber,1992)was adopted to estimate the background root-mean-square error(RMSE)and the three Kronecker product members of the background error correlation matrix.The NMC method utilized the difference between the 12-and 24-h WRF-Chem forecasts valid at the same time of 1200 UTC for a whole month.No cross correlation between different species was assigned for background error.The domain-average RMSEs for five species are shown in Fig.2.The vertical distributions of RMSE for all species display a relatively rapid decrease with height—except for O3,which peaks at around 4 km above the ground.The vertical correlation matrices are displayed in Fig.3.The correlation between different height levels experiences a jump at the top of the boundary layer,which seems to be a common feature for all species.In addition,the band of high correlation along the diagonal seems wider in the middle troposphere than the upper or lower troposphere.

    Considering that all the stations were built and maintained under the same standards,no difference in the measurement and representativeness error between different stations was assumed.In addition,cross correlation between different species and stations was set to zero because of the lack of information.Therefore,observation error consisted merely of five values—one for each species.The measurement error was assigned as 1.0μg m?3,1.0μg m?3,1.0 ppb,1.0 ppb and 1.0 ppb for PM2.5,PM2.5–PM10,SO2,NO2and O3,respectively.Representativeness error was estimated following Elbern et al.(2007)and Schwartz et al.(2012)using the formula

    Fig.2.Vertical distribution of the root-mean-square ofthebackground errors,in mass concentration for particulate matter and ppb for gaseous species.

    Fig.3.Vertical correlations of the background errors for PM2.5,SO2,NO2,and O3.The plot for PM2.5–PM10is very similar to that of PM2.5and so is not shown here.Both the x-and y-axis are limited within 10 km,though the real model top is higher.

    where ε0and εrare the measurement error and representativeness error,γ is an adjustable parameter that accounts for the lifetime of the species(0.5 for PM2.5,PM2.5–PM10and O3;1 for SO2and 2 for NO2),Δxis the grid spacing(here,15 km),andLis the radius of in fluence determined according to the location of stations(here,4.0 km for suburban stations—assumed for all sites).If the total observation error is defined as the sum of measurement error and representativeness error,the standard deviation of observation error is 2.0μg m?3,2.0μg m?3,3.0 ppb,4.9 ppb and 2.0 ppb for PM2.5,PM2.5–PM10,SO2,NO2and O3,respectively.Though the observation error was determined fairly arbitrarily and empirically here,the uncertainty relating to it should not have a signi ficant in fluence on the conclusion.That is because the results are usually not very sensitive to the speci fication of error(Geer et al.,2006),and similar analysis fields were obtained from our experiments when the observation error was increased or reduced by a factor of two or three.

    2.4.Experimental design

    To compare the relative importance of MOS and DA,four parallel experiments were designed:Sim base,SimDA,Sim-MOS and Sim-DM.Simbase worked as the base simulation without applying DA or MOS;Sim-DA was an experiment with only DA employed to optimize the model CICs;Sim-MOS was the same as Sim-base but with the model output corrected by one-dimensional Kalman filtering;and Sim-DM used both the DA and MOS methods.

    To simulate the operational forecast scenes,as Fig.4 shows,all experiments initiated a new WRF-Chem forecast at 1200 UTC,every day,between 30 November 2014 and 31 December 2014.Each forecast was integrated for 84 h to generate 72-h forecasts for each day,with the earliest 12 h discarded as spin-up time.The CICs for each initiation came from the 24-h forecasts of the previous cycle,which would be the background fields to be assimilated with valid observations for experiments with DA before initializing WRFChem.The first CICs at 1200 UTC 30 November 2014 came from a two-days’spin-up of the climatological background chemistry pro file.For all experiments,the chemical boundary conditions came from the default climatological chemistry pro file for the outer domain and the interpolation of the outer domain for the inner domain.

    Fig.4.Time settings of the model in the four experiments.In each forecast cycle,the model was integrated 84 h in advance,with the first 12 h discarded(black dashed arrow)and the remaining 72 h divided into three parts:24-h forecast(black solid arrow),48-h forecast(blue solid arrow)and 72-h forecast(red solid arrow).The blue thick arrows mean using CICs directly for experiments without DA,while for others as the background fields of the 3DVar DA.

    2.5.Observational data

    Hourly concentrations of SO2,NO2,PM10,O3and PM2.5at surface level from 207 sites were provided by the Ministry of Environmental Protection of China.Data covered the whole month of December 2014 and had been subjected to routine quality control.As shown in Fig.5,only 155 stations were selected(randomly)from the 207 stations to be assimilated,and the data of the remaining 52 were used to verify the assimilation process.Because all stations were located at surface level,the adjustment of the CICs from the 3DVar DA was limited within several layers near the surface according to the vertical background error covariance.Furthermore,it should be noted that only those 155 sites that provided their data for the 3DVar DA participated in the MOS process.

    3.Results

    3.1.Model evaluation

    Table 2 presents the mean bias(MB),relative bias(RB),RMSE and correlation coefficient(Corr)for the 24-h,48-h and 72-h forecast of Sim base.In general,the base model simulation provides a fairly good result—especially for NO2,whose bias is small and correlation high.In terms of particulate matter and SO2,the model tends to systematically underestimate the concentration of SO2as well as that of PM2.5and PM10.Even so,the model reproduces the temporal variations of particulate matter well,with Corr values higher than 0.47 for PM10and 0.54 for PM2.5.For O3,the model encounters a problem—the surface O3simulated concentrations(5–45μg m?3from observation)are seriously overestimated by the model(20–80 μg m?3from simulation),which leads to positive bias(~ 40 μg m?3)and lower Corr(0.44)than for other species.Fortunately,when viewing the RMSE of O3from the aspect of MB,it is apparent that MB contributes the largest portion of RMSE,and therefore the model is still able to reproduce the variation of O3.The biases mentioned above can usually be attributed to the uncertainties from the emissions inventory,meteorological forecasting(Tang et al.,2011)and model schemes(Yerramilli et al.,2010).Although the 24-h forecast performs the best for all species,the 48-h and 72-h forecasts are also good enough to yield fairly reliable results,which is critical to the success of MOS in the whole 72 hours’forecast.In short,the model shows forecast skill that is sufficient to be competent for the success of the DA and MOS process.

    Fig.5.The terrain height of the Hebei Province,with monitoring sites plotted as filled dots.Red dots are sites that participated in the 3DVar and MOS process,while the blue ones are those used only in the validation of the DA effect.

    Table 2.Site-averaged mean bias(MB),relative bias(RB),root mean square error(RMSE)and correlation coefficient(Corr)for the 24-h,48-h and 72-h forecasts of five species in Sim base.Statistics were calculated according to the hourly concentration from each of the 207 sites before being averaged.Units areμg m?3for MB and RMSE.

    Table 3.Site-averaged MB,RMSE and Corr for the 24-h,48-h and 72-h forecasts of five species in Sim MOS.Statistics were calculated according to the hourly concentration from each of the 207 sites before being averaged.Units areμg m?3for MB and RMSE.

    3.2.Validation of MOS

    Figure 6 depicts the site-averaged hourly concentration simulated by Sim MOS,plotted against ground observations.Note that,although the hourly concentrations were averaged over 155 stations,the MB and RMSE values in Table 3 were generated by first calculating the individual errors of the 155 stations,before averaging.

    From Fig.6 it can be concluded that the forecast from Sim MOS fits the observation fairly closely—especially for SO2,NO2and O3.However,when it comes to PM2.5and PM10,the points locate within a wider space,and those extremely high observations are hard for MOS to forecast.Even so,when comparing Table 3 with Table 2,PM2.5and PM10,together with the other three species,demonstrate that a clear correction can be obtained for all forecast times.Excluding the 48-h forecast of NO2,MOS can reduce the MB to a large extent,meaning this method can remove the majority of the model systematic bias.Because of the reduction in MB,the RMSE also decreases for all cases except the 24-h forecast of SO2.In addition,the effect on reducing the error is unlikely to become poorer as the forecast time advances.For the 72-h and 48-h forecasts,the effect MOS has on the former may rival or even exceed that on the latter,e.g.,the RMSE reduction of PM2.5is even larger for the 72-h than 48-h forecast.For the 48-h and 24-h forecasts,the same effect can be found.Among all five species,O3seems to bene fit the most from the MOS process.This is because O3usually follows a very regular daily variation,which makes the hourly-split but dailylinked concentration series almost perfect for the assumptions of one-dimensional Kalman filtering.

    Admittedly,MOS degrades the forecast in a few cases(e.g.,the 48-h forecast of NO2and 24-h forecast of SO2,as mentioned above).Such increases in error,however,will usually not be of concern to users,and may well be accepted,as they are extremely small and only appear at times when the model outputs to be corrected are already fairly close to the observation.Nonetheless,when viewed from the correlation perspective,such degradation becomes moreobvious.Except for NO2and O3,the correlations all experience a reduction by 0.1–0.2.Thus,the MOS approach tends to reduce the bias and error at the expense of correlation.

    3.3.Validation of DA

    Figures 7 and 8 show the change in RMSE and Corr over the integration time from?12 h(right after the DA)to 10 h(already integrated for 22 h),respectively,for experiments Sim base and SimDA.To keep the veri fication independent from the observations assimilated,the RMSE and Corr were only averaged over the 52 stations that did not provide their observational data in the 3DVar DA process.

    Fig.6.Hourly concentrations simulated by Sim MOS,plotted against ground station observations(obs)averaged over 155 stations.

    Fig.7.RMSE change over the integration time from the?12 h forecast(right after the DA)to the 10 h forecast(already integrated for 22 h)for the five species.All RMSE values were averaged over the 52 stations that did not provide observations in the DA.

    From the?12 h forecast of Fig.7 and Fig.8,DA leads to better initial conditions for the simulation—especially for NO2,PM10,PM2.5and O3,whose RMSEs decrease substantially at almost all sites.For example,PM10and PM2.5shown an RMSE reduction of about 50–100 μg m?3,which is about half the RMSE of Simbase.Such results are as good as those obtained by L13 and Jiang et al.(2013),who also worked on assimilating ground observations using 3DVar.For SO2,the reduction in RMSE is less apparent(although the change in RMSE is still negative when 52 sites are averaged),but the correlation after DA is still obviously larger than before.The marginal RMSE reduction for SO2may be acceptable considering the increase in correlation is rather obvious and the data representativeness of some stations is dubious(Zhang et al.,2016).

    Fig.8.Similar to Fig.7 but for Corr.

    However,as expected,the effect of DA slowly diminishes as the integration continues,which has also been observed in other works(Jiang et al.,2013;Li et al.,2013).After the model has been integrated for more than 14 h,the RMSE after DA minus that before DA(RMSE change henceforth)is still negative,but their absolute values are apparently smaller when compared to the earlier ones.The effect of DA remains for longer in the case of O3,PM10and PM2.5(RMSE change remains negative for the 14 h),benefiting mainly from their relatively long lifetimes.For example,PM2.5and PM10still maintain an RMSE reduction of about 10–20 μg m?3,which is even better than the results reported by L13 and Jiang et al.(2013).However,for NO2,whose lifetime is short,the two experiments show almost no difference in RMSE after four hours’integration.Because the initial improvement from DA is relatively small,the forecast of SO2soon loses its improvement from DA and shows little RMSE change almost immediately after the run of the model.For SO2and NO2,the RMSE change is positive in some cases,but the increase in RMSE is usually very small compared with the original RMSE,and therefore unlikely to be a problem.Conclusions from the viewpoint of correlation are similar to those from the RMSE,except the effect of DA seems more obvious and long-lasting.

    Overall,for most cases,DA successfully produces better CICs for the model and may help to improve the forecast skill in the following half to one day.

    3.4.Effect of MOS and DA

    From Fig.9,which plots the RMSE of the four experiments and five species at different forecast hours from ?12 h to 72 h,we can see that the forecast error varies with forecast hours.Considering the RMSE was calculated from the statistics of one month,Fig.9 is a good representation of the general features of the four experiments.

    By comparing the experiments using MOS(solid lines)with those without MOS(dotted lines),it is clear that all species show a large reduction in the 72 hours forecast span,and that this reduction is much larger than can be provided by DA(solid lines are below the dotted lines by a larger distance than the blue lines are below the red).For example,when compared to Sim base,the average SO2RMSE decreases by 4.34 μg m?3for Sim MOS(average taken over all stations throughout the whole 72-h forecast),while the decrease is only 0.48μg m?3for Sim DA.Worse still,when the forecast runs to its second or third day,the effect of DA inevitably diminishes(as evidenced by the overlapping red and blue dotted lines after 24 h),while MOS can still work during this period(solid lines do not overlap the dotted lines,even after 24 h).

    The blue dotted lines represent the simulation RMSE corrected only with MOS,while the red dotted lines are the results processed by both MOS and DA.Overlapping of the two lines can be seen at almost all times from all species,meaning the DA of the initial conditions provides little help to the effect of MOS,although does provide better initial conditions for the model to generate a better forecast.However,sometimes the two lines do not overlap,and show some differences,which is common for all species but most obvious for NO2and most unobvious for PM10and PM2.5.When DA can still improve the forecast,or the red dotted line is below the blue dotted line,the red solid line could be either above(1-h forecast for SO2)or below(10-h forecast for O3)the blue line,which means a better forecast from DA may either improve or degrade the MOS effect.Because in this work MOS corrects one day’s forecast using the correction results and forecast from previous days,it is not a surprise that Sim DM and SimMOS show discrepancies when SimDA and Simbase coincide after 24 h.However,as mentioned,such a discrepancy could be either an improvement or degradation.

    3.5.Discussion

    Fig.9.RMSE of the four experiments and five species at different forecast hours from ?12 h to 72 h.

    This section provides additional discussion around two facts.First,that MOS may improve forecasts far more than the 3DVar DA of CICs.This is reasonable because MOS can remain effective throughout the whole forecast period,whereas the effect of 3DVar via optimized initial conditions usually diminishes after 24 h of model integration.The loss of bene fit from 3DVar DA is unavoidable because atmospheric chemistry is less sensitive to CICs than other driving factors like meteorological conditions and emissions(Henze et al.,2009;Semane et al.,2009;Tang et al.,2011).Worse still,the forecast during the earliest 12 h,which bene fits the most from 3DVar DA,usually makes no sense in a real operational forecasting environment and is excluded from evaluation as spin-up.In fact,when compared with Sim base,Sim DA can account for 43.85%of the O3RMSE reduction in the first 12 h after initialization(forecast hour?12 to?1),closely rivalling the contribution of Sim MOS(55.94%)in its first 12 h(forecast hour 0 to 11).However,when discussed within the same forecast period,e.g.,12 h to 24 h,Sim DA can produce only 3.93%of the O3RMSE reduction,which is far less than Sim MOS(61.26%),despite the following hours during which DA has no effect at all.

    The second fact to be explained is that Sim DM does not al way so ut per form Sim MOS.This result is somehow against the experience that,when corrected with the same MOS algorithm,better input should lead to better or at least not worse output.However,for example,assume the difference between the forecast and observation remains constant before 3DVar is applied.(This condition may never happen in reality,but this does not matter for the purposes of demonstration).Then,no matter how large the constant is,the MOS method might work perfectly to eliminate it.However,after DA is applied,it is possible that the difference may reduce though no longer remain constant.In this case,although the forecast becomes better before MOS,it is now more difficult for MOS to eliminate the error.So,it should be noted that the error’s temporal consistency,rather than its magnitude,determines the effects of MOS on the model outputs.When reducing the magnitude of the error of model outputs,the 3DVar DA process may at the same time violate or increase its consistency to degrade or improve the effects of MOS from case to case.Therefore,such a phenomenon is uncorrelated with the inherent or necessary nature of the model,DA and MOS processes,and will be changed randomly whenever the three change their setup.The assumption is supported by the evidence that the results are very different when the whole experiment is replicated but with the spatial resolution of the original anthropogenic emissions changed from 0.1°×0.1°to 0.25°×0.25°[see Fig.10(PM2.5and PM10not plotted given the problem for them was not obvious in Fig.9)].For example,in Fig.9,SO2is predicted better by Sim DA than Sim base at the second forecast hour,but Sim DM is beaten by Sim MOS.However,in Fig.10,the same comparison leads to an inverse result,which shows Sim DM performs better than Sim MOS.Given the fact that the experiments cover a period of only 1 month,it is possible that the forecast ability of Sim MOS is slightly worse than, or almost the same as,Sim DM,if the experiment is carried out over a longer time.However,results from short-term experiments,which contain random error just like everyday forecast,still demonstrate that using MOS and DA together does not guarantee a better output than MOS only,which should be carefully considered by researchers and forecasters working in this field.

    4.Conclusion

    A comparison between the effect of MOS and DA in improving the forecast skill of an atmospheric chemistry model was performed in a near real operational context.The evaluation based on observations showed that both 3DVar DA and MOS based on one-dimensional Kalman filtering are effective measures to reduce the errors in the model forecast.

    Fig.10.RMSE of three species at different forecast hours from ?12 h to 72 h from the same four experiments with solely emissions perturbed.

    The forecast solely with MOS(SimMOS)performed better than that solely with 3DVar DA of CICs(SimDA).Such superiority of MOS could also be seen in all five species.That is,the implementation of MOS rather than 3DVar DA on CICs is more suitable for the aim of improving the operational forecast ability.

    Considering the randomness of the influence of DA on error consistency,it is not impossible that combined use of both techniques will sometimes yield a worse forecast than MOS only.The potential degradation,which is expected to be mitigated by long-term averaging,should be considered,but is unlikely to concern forecasters because of the relatively limited difference yielded.However,considering the complexity of the two approaches,the feature of MOS,rather than DA of CICs,is recommended first when planning to improve your model forecast.

    Given that refinement of the model grid is promising for additional forecast skill(Elbern et al.,2007),future work will concentrate on testing finer model grids in the vertical and horizontal dimensions,as well as different 3DVar observation error setups,in order to improve the effect of DA on the forecast.Of particular interest are species like SO2and NO2,which are unsatisfactorily forecasted by models with 3DVar.Also,to see if the conclusions will be different,we intend to try different DA methods,such as 4DVar and inverse modeling,which are able to adjust model parameters and emissions,which work as control parameters in a successful forecast(Schmidt and Martin,2003;Dubovik et al.,2008).Finally,another line of research would be extending the experiment time to explore whether any long-term statistical properties exist when using the DA and MOS techniques together.

    Acknowledgements.This work was supported by the State Key Research and Development Program(Grant Nos.2017YFC0209803,2016YFC0208504,2016YFC0203303and 2017YFC0210106)and the National Natural Science Foundation of China(Grant Nos.91544230,41575145,41621005 and 41275128).The Hebei Meteorological Bureau of China is also thanked for making available the pollution monitoring data.

    APPENDIX

    Assume we have a systemxthat varies with timet,andxat timetis almost the same asxat timet?1 expect for a random errorw.Then,this system can be described by a system equation such as

    For the systemx,assume we have some kind of measurementythat measuresxwith a random errorv.Then,we can write the measurement equation as

    Ifwandvare both white noise,with their covariance already known asWandVand independent of each other,then we can obtain the best estimate ofxusing the following equations:

    In these equations,?x(t)stands for the best estimate of the real systemx(t),which is usually unknown;andK(t)andP(t)are both intermediate variables whose meanings are not important.To complete a single one-dimensional Kalman filter iteration,the inputs of the equations areP(t?1)and?x(t?1)from the previous calculation,as well as a newly obtained measurementy(t).Then,the equations will yield newP(t)and?x(t),which will be saved for the next iteration.A detailed deduction of Eq.(A3)is beyond the scope of this work;interested readers should instead refer to Kalman(1960).

    Here,three sequences of O3concentrations at 1200 UTC from three days’model output will be taken as an example to show the detailed steps of completing one iteration of MOS.It is assumed that on dayt(here,trepresents “today”),the following variables had been prepared:observation sequenceOi:i=1,2,···,tfrom the ground station;model 24-h forecast sequence and its correctionfi,24:i=1,2,···,t+1 and ?fi,24:i=1,2,···,t;model 48-h forecast sequence and its correctionandmodel 72-h forecast sequence and its correctionand:i=1,2,···,t+2(model forecast sequence can be obtained by interpolating the 3D field of the model output to the position of the site);as well as all thePand?xfrom the previous iteration.What we want is to generatefrom all variables above.

    To obtaintwo different approaches will be applied simultaneously.In one approach,which is the same as in G02,x(i)=Oi?fi,24:i=1,2,···,tand1,2,···,tin Eqs.(A2)and(A3),so this approach is called the “difference approach”hereafter.v(i)in Eq.(A3)is set to zero to assume no measurement error,which meansVin Eq.(A3b)equals zero and the estimation ofWwould be unnecessary here.The calculation of Eq.(A3a)to Eq.(A3b)gives thethen,we assumexandThe subscript“d”means the corrected result from the difference approach.Additionally,we add another approach called the “ratio approach”.In this approach,all formulas are the same except(i)=Oi/fi,24:i=1,2,···,t,y(i)=Oi/fi,24:i=1,2,···,tandwhere the subscript“r”means the ratio approach.The final resultwill be chosen fromaccording to the method described in the final paragraph of this appendix.

    Withobtained,it is now possible to calculatewe still carry out the calculation via two approaches:the difference and ratio approach.Now,andy(i)So,to some extent,works as a “measurement”of the real observationOi,or a“proxy”of observation as in section 2.2.At this time,v(i)cannot be set to zero in Eq.(A3)and the evaluation ofVandWbecomes necessary in completing the iteration.The estimation ofWandVresembles G02,in which

    The variablenin Eqs.(A4)and(A5)is the number of sequence members participating in the statistics.A value from 7 to 9 is enough fornto generate a fairly good MOS effect,which means historical data from the past 7–9 days’observations and model output works as the virtual training sample for MOS.After completing the calculation of Eq.(A3a)–(A3d),?x(t+1)will be obtained and then?x(t+2)andcan be produced by replicating the process from the paragraph above.Another ratio approach is also performed to generateis chosen fromin the same way described in the final paragraph.

    Fig.A1.Steps in applying the one-dimensional Kalman filter to correction over the whole 72 hours forecast period.The circles from left to right stand for values at the same hour from different days.Black circles are those existing before today’s calculation,while red ones are those generated from the calculation of today.

    Fromto obtainingwhat has to be done is almost the same as fromThe differences mainly lie in some of the formulations.Taking the difference approach as an example,x(i)=Oi?fi,72:i=1,2,···,t,and the estimation ofVandWis changed to

    Figure A1 is a schematic diagram showing how the algorithm described above is carried out,step by step.

    The final output should be chosen from the difference and ratio approach according to their reasonability.From our experience,the difference approach tends to occasionally yield unreasonably low values,while the ratio approach sometimes gives results that are too high.Fortunately,the two conditions never happen simultaneously,and therefore the final output is from the ratio approach when its result is not too high(lower than the yearly averaged value of the species,for example).In cases where the ratio approach appears too high,the difference approach should be used instead.

    REFERENCES

    Ackermann,I.J.,H.Hass,M.Memmesheimer,A.Ebel,F.S.Binkowski,and U.Shankar,1998:Modal aerosol dynamics model for Europe:Development and first applications.Atmos.Environ.,32,2981–2999,http://dx.doi.org/10.1016/S1352-2310(98)00006-5.

    Alexe,M.,and Coauthors,2015:Inverse modelling of CH4emissions for 2010–2011 using different satellite retrieval products from GOSAT and SCIAMACHY.Atmos.Chem.Phys.,15,113–133,http://dx.doi.org/10.5194/acp-15-113-2015.

    Anadranistakis,M.,K.Lagouvardos,V.Kotroni,and H.Elefteriadis,2004:Correcting temperature and humidity forecasts using Kalman filtering:Potential for agricultural protection in Northern Greece.Atmos.Res.,71,115–125,http://dx.doi.org/10.1016/j.atmosres.2004.03.007.

    Barbu,A.L.,A.J.Segers,M.Schaap,A.W.Heemink,and P.J.H.Builtjes,2009:A multi-component data assimilation experiment directed to sulphur dioxide and sulphate over Europe.Atmos.Environ.,43,1622–1631,http://dx.doi.org/10.1016/j.atmosenv.2008.12.005.

    Barret,B.,and Coauthors,2008:Transport pathways of CO in the African upper troposphere during the monsoon season:A study based upon the assimilation of spaceborne observations.Atmos.Chem.Phys.,8,3231–3246,http://dx.doi.org/10.5194/acp-8-3231-2008.

    Benedetti,A.,and Coauthors,2009:Aerosol analysis and forecast in the European Centre for Medium-Range Weather Forecasts Integrated Forecast System:2.Data assimilation.J.Geophys.Res.,114,D13205,http://dx.doi.org/10.1029/2008JD011115.

    Bocquet,M.,and Coauthors,2015:Data assimilation in atmospheric chemistry models:Current status and future prospects for coupled chemistry meteorology models.Atmos.Chem.Phys.,15,5325–5358,http://dx.doi.org/10.5194/acp-15-5325-2015.

    Denby,B.,M.Schaap,A.Segers,P.Builtjes,and J.Hor′alek,2008:Comparison of two data assimilation methods for assessing PM10exceedances on the European scale.Atmos.Environ.,42,7122–7134,http://dx.doi.org/10.1016/j.atmosenv.2008.05.058.

    Dubovik,O.,T.Lapyonok,Y.J.Kaufman,M.Chin,P.Ginoux,R.A.Kahn,and A.Sinyuk,2008:Retrieving global aerosol sources from satellites using inverse modeling.Atmos.Chem.Phys.,8,209–250,http://dx.doi.org/10.5194/acp-8-209-2008.

    Elbern,H.,A.Strunk,H.Schmidt,and O.Talagrand,2007:Emission rate and chemical state estimation by 4-dimensional variational inversion.Atmos.Chem.Phys.,7,3749–3769,http://dx.doi.org/10.5194/acp-7-3749-2007.

    Galanis,G.,and M.Anadranistakis,2002:A one-dimensional Kalman filter for the correction of near surface temperature forecasts.Meteorological Applications,9,437–441,http://dx.doi.org/10.1017/S1350482702004061.

    Geer,A.J.,and Coauthors,2006:The ASSET in tercomparison of ozoneanalyses:Method and first results.Atmos.Chem.Phys.,6,5445–5474,http://dx.doi.org/10.5194/acp-6-5445-2006.

    Glahn,H.R.,and D.A.Lowry,1972:The use of model output statistics(MOS)in objective weather forecasting.J.Appl.Meteor.,11,1203–1211,http://dx.doi.org/10.1175/1520-0450(1972)011<1203:TUOMOS>2.0.CO;2.

    Grell,G.A.,S.E.Peckham,R.Schmitz,S.A.McKeen,G.Frost,W.C.Skamarock,and B.Eder,2005:Fully coupled “online”chemistry within the WRF model.Atmos.Environ.,39,6957–6975,http://dx.doi.org/10.1016/j.atmosenv.2005.04.027.

    Henze,D.K.,J.H.Seinfeld,and D.T.Shindell,2009:Inverse modeling and mapping US air quality in fluences of inorganic PM2.5 precursor emissions using the adjoint of GEOSChem.Atmos.Chem.Phys.,9,5877–5903,http://dx.doi.org/10.5194/acp-9-5877-2009.

    Honore,C.,and Coauthors,2008:Predictability of European air quality:Assessment of 3 years of operational forecasts and analyses by the PREV’AIR system.J.Geophys.Res.,113,http://dx.doi.org/10.1029/2007JD008761.

    Inness,A.,and Coauthors,2015:Data assimilation of satelliteretrieved ozone,carbon monoxide and nitrogen dioxide with ECMWF’s Composition-IFS.Atmos.Chem.Phys.,15,5275–5303,http://dx.doi.org/10.5194/acp-15-5275-2015.

    Jiang,Z.Q.,Z.Q.Liu,T.J.Wang,C.S.Schwartz,H.C.Lin,and F.Jiang,2013:Probing into the impact of 3DVAR assimilation of surface PM10 observations over China using process analysis.J.Geophys.Res.,118,6738–6749,http://dx.doi.org/10.1002/jgrd.50495.

    Kalman,R.E.,1960:A new approach to linear filtering and prediction problems.Journal of Basic Engineering,82,35–45,http://dx.doi.org/10.1115/1.3662552.

    Li,Z.,Z.Zang,Q.B.Li,Y.Chao,D.Chen,Z.Ye,Y.Liu,and K.N.Liou,2013:A three-dimensional variational data assimila-tion system for multiple aerosol species with WRF/Chem and an application to PM2.5 prediction.Atmos.Chem.Phys.,13,4265–4278,http://dx.doi.org/10.5194/acp-13-4265-2013.

    Libonati,R.,I.Trigo,and C.C.Dacamara,2008:Correction of 2m-temperature forecasts using Kalman Filtering technique.Atmos.Res.,87,183–197,http://dx.doi.org/10.1016/j.atmosres.2007.08.006.

    Liu,Z.Q.,Q.H.Liu,H.C.Lin,C.S.Schwartz,Y.H.Lee,and T.J.Wang,2011:Three-dimensional variational assimilation of MODIS aerosol optical depth:Implementation and application to a dust storm over East Asia.J.Geophys.Res.,116,D23206,http://dx.doi.org/10.1029/2011JD016159.

    Makar,P.A.,and Coauthors,2015:Feedbacks between air pollution and weather,Part 1:Effects on weather.Atmos.Environ.,115,442–469,http://dx.doi.org/10.1016/j.atmosenv.2014.12.003.

    Mizzi,A.P.,A.F.Arellano Jr.,D.P.Edwards,J.L.Anderson,and G.G.P fister,2016:Assimilating compact phase space retrievals of atmospheric composition with WRF-Chem/DART:Aregional chemical transport/ensembleKalman filterdataassimilation system.Geoscientific Model Development,9,965–978,http://dx.doi.org/10.5194/gmd-9-965-2016.

    Parrish,D.F.,and J.C.Derber,1992:The National Meteorological Center’s spectral statistical-interpolation analysis system.Mon.Wea.Rev.,120,1747–1763,https://doi.org/10.1175/1520-0493(1992)120<1747:TNMCSS>2.0.CO;2.

    Saide,P.E.,G.R.Carmichael,S.N.Spak,P.Minnis,and J.K.Ayers,2012:Improving aerosol distributions below clouds by assimilating satellite-retrieved cloud droplet number.Proceedings of the National Academy of Sciences of the United States of America,109,11 939–11 943,http://dx.doi.org/10.1073/pnas.1205877109.

    Saide,P.E.,and Coauthors,2015:Central American biomass burning smoke can increase tornado severity in the U.S.Geophys.Res.Lett.,42,956–965,http://dx.doi.org/10.1002/2014GL062826.

    Schell,B.,I.J.Ackermann,H.Hass,F.S.Binkowski,and A.Ebel,2001:Modeling the formation of secondary organic aerosol within a comprehensive air quality model system.J.Geophys.Res.,106,28 275–28 293,http://dx.doi.org/10.1029/2001JD000384.

    Schmidt,H.,and D.Martin,2003:Adjoint sensitivity of episodic ozone in the Paris area to emissions on the continental scale.J.Geophys.Res.,108,8561,http://dx.doi.org/10.1029/2001JD001583.

    Schwartz,C.S.,Z.Q.Liu,H.C.Lin,and S.A.McKeen,2012:Simultaneous three-dimensional variational assimilation of surface fine particulate matter and MODIS aerosol optical depth.J.Geophys.Res.,117,D13202,http://dx.doi.org/10.1029/2011JD017383.

    Semane,N.,and Coauthors,2009:On the extraction of wind information from the assimilation of ozone pro files in M′et′eo-France 4-D-Var operational NWP suite.Atmos.Chem.Phys.,9,4855–4867,http://dx.doi.org/10.5194/acp-9-4855-2009.

    Stock well,W.R.,P.Middleton,J.S.Chang,and X.Y.Tang,1990:The second generation regional acid deposition model chemical mechanism for regional air quality modeling.J.Geophys.Res.,95,16 343–16 367,http://onlinelibrary.wiley.com/doi/10.1029/JD095iD10p16343/full.

    Struzewska,J.,J.W.Kaminski,and M.Je fimow,2016:Application of model output statistics to the GEM-AQ high resolution air quality forecast.Atmos.Res.,181,186–199,http://dx.doi.org/10.1016/j.atmosres.2016.06.012.

    Tang,X.,J.Zhu,Z.F.Wang,and A.Gbaguidi,2011:Improvement of ozone forecast over Beijing based on ensemble Kalman filter with simultaneous adjustment of initial conditions and emissions.Atmos.Chem.Phys.,11,12 901–12 916,http://dx.doi.org/10.5194/acp-11-12901-2011.

    Tang,X.,J.Zhu,Z.F.Wang,A.Gbaguidi,C.Y.Lin,J.Y.Xin,T.Song,and B.Hu,2016:Limitations of ozone data assimilation with adjustment of NOxemissions:Mixed effects on NO2forecasts over Beijing and surrounding areas.Atmos.Chem.Phys.,16,6395–6405,http://dx.doi.org/10.5194/acp-16-6395-2016.

    Taylor,A.A.,and L.M.Leslie,2005:A single-station approach to model output statistics temperature forecast error assessment.Wea.Forecasting,20,1006–1020,http://dx.doi.org/10.1175/WAF893.1.

    van Loon,M.,and Coauthors,2007:Evaluation of long-term ozone simulations from seven regional air quality models and their ensemble.Atmos.Environ.,41,2083–2097,http://dx.doi.org/10.1016/j.atmosenv.2006.10.073.

    Wang,Y.,K.N.Sartelet,M.Bocquet,and P.Chazette,2014:Modelling and assimilation of lidar signals over Greater Paris during the MEGAPOLI summer campaign.Atmos.Chem.Phys.,14,3511–3532,http://dx.doi.org/10.5194/acp-14-3511-2014.

    Wargan,K.,S.Pawson,M.A.Olsen,J.C.Witte,A.R.Douglass,J.R.Ziemke,S.E.Strahan,and J.E.Nielsen,2015:The global structure of upper troposphere-lower stratosphere ozone in GEOS-5:A multiyear assimilation of EOS Aura data.J.Geophys.Res.,120,2013–2036,http://dx.doi.org/10.1002/2014JD022493.

    Wilson,L.J.,and M.Vall′ee,2003:The Canadian Updateable Model Output Statistics(UMOS)system:Validation against perfect prog.Wea.Forecasting,18,288–302,http://dx.doi.org/10.1175/1520-0434(2003)018<0288:TCUMOS>2.0.CO;2.

    Yerramilli,A.,and Coauthors,2010:Simulation of surface ozone pollution in the central gulf coast region using WRF/Chem model:Sensitivity to PBL and land surface physics.Advances in Meteorology,2010,Article ID 319138,http://dx.doi.org/10.1155/2010/319138.

    Yin, X. M., T. Dai, N. A. J. Schutgens, D. Goto, T. Nakajima,and G. Y. Shi, 2016: Effects of data assimilation on the global aerosol key optical properties simulations.Atmos.Res.,178–179,175–186,https://doi.org/10.1016/j.atmosres.2016.03.016.

    Yumimoto,K.,and Coauthors,2016:Aerosol data assimilation using data from Himawari-8,a next-generation geostationary meteorological satellite.Geophys.Res.Lett.,43,5886–5894,http://dx.doi.org/10.1002/2016GL069298.

    Zhang,L.,and Coauthors,2015:Source attribution of particulate matter pollution over North China with the adjoint method.Environmental Research Letters,10,084011,https://doi.org/10.1088/1748-9326/10/8/084011.

    Zhang,L.,and Coauthors,2016:Sources and processes affecting fine particulate matter pollution over North China:An adjoint analysis of the Beijing APEC Period.Environ.Sci.Technol.,50,8731–8740,http://dx.doi.org/10.1021/acs.est.6b03010.

    久久久国产成人免费| 久久人妻熟女aⅴ| 性色av乱码一区二区三区2| 亚洲人成电影免费在线| 精品久久蜜臀av无| 脱女人内裤的视频| 国产高清视频在线播放一区| 精品少妇黑人巨大在线播放| 久久久久国内视频| 法律面前人人平等表现在哪些方面| 在线观看舔阴道视频| 国产精品久久久久久精品古装| 97在线人人人人妻| 成年动漫av网址| 国内毛片毛片毛片毛片毛片| 久久中文字幕一级| 久久午夜亚洲精品久久| 亚洲第一欧美日韩一区二区三区 | 一级毛片精品| 老司机亚洲免费影院| 国产精品影院久久| 黄色视频,在线免费观看| 久久久国产成人免费| 国产精品 欧美亚洲| 大码成人一级视频| 美女扒开内裤让男人捅视频| 亚洲欧美一区二区三区黑人| 欧美日本中文国产一区发布| 亚洲av美国av| 露出奶头的视频| 女性生殖器流出的白浆| 人人妻人人澡人人爽人人夜夜| 国产一区二区三区在线臀色熟女 | 日本欧美视频一区| 亚洲精品在线美女| 欧美精品一区二区免费开放| 精品免费久久久久久久清纯 | 人人妻人人澡人人爽人人夜夜| 亚洲精品一二三| 女同久久另类99精品国产91| 欧美一级毛片孕妇| 国产欧美日韩综合在线一区二区| 美女高潮喷水抽搐中文字幕| 国产精品麻豆人妻色哟哟久久| 久久久久久久国产电影| 国产亚洲午夜精品一区二区久久| 午夜久久久在线观看| 大陆偷拍与自拍| 国产精品香港三级国产av潘金莲| 欧美精品啪啪一区二区三区| 亚洲国产欧美在线一区| 视频在线观看一区二区三区| 日韩欧美免费精品| 大香蕉久久网| 蜜桃在线观看..| 国产精品麻豆人妻色哟哟久久| 日韩欧美国产一区二区入口| 侵犯人妻中文字幕一二三四区| 欧美精品啪啪一区二区三区| 精品久久久精品久久久| 丁香六月欧美| 亚洲综合色网址| 最新在线观看一区二区三区| 国产精品久久久久久人妻精品电影 | 亚洲av片天天在线观看| 夜夜夜夜夜久久久久| 岛国毛片在线播放| 十八禁网站网址无遮挡| 在线观看免费视频网站a站| 欧美+亚洲+日韩+国产| 蜜桃国产av成人99| 国产精品久久久久久精品古装| 美国免费a级毛片| 91精品国产国语对白视频| 久久久精品94久久精品| 大码成人一级视频| 精品久久蜜臀av无| 国产色视频综合| 无人区码免费观看不卡 | 亚洲国产看品久久| 中文字幕高清在线视频| 日韩免费高清中文字幕av| 丁香六月欧美| 国产成人av激情在线播放| 亚洲国产av影院在线观看| 亚洲免费av在线视频| 久久香蕉激情| 欧美成狂野欧美在线观看| 国产成人精品久久二区二区91| 亚洲成av片中文字幕在线观看| 99re在线观看精品视频| 人人妻,人人澡人人爽秒播| 国产成人欧美| 国产欧美日韩一区二区三区在线| 亚洲欧美日韩另类电影网站| 王馨瑶露胸无遮挡在线观看| 国产色视频综合| 汤姆久久久久久久影院中文字幕| 高清黄色对白视频在线免费看| 天堂动漫精品| 亚洲男人天堂网一区| 日韩精品免费视频一区二区三区| 亚洲精品在线观看二区| 国产深夜福利视频在线观看| 悠悠久久av| 欧美 亚洲 国产 日韩一| 欧美一级毛片孕妇| 国产精品久久久久久人妻精品电影 | 女人被躁到高潮嗷嗷叫费观| 老司机在亚洲福利影院| 50天的宝宝边吃奶边哭怎么回事| 老司机午夜十八禁免费视频| 亚洲黑人精品在线| 1024香蕉在线观看| 99re在线观看精品视频| av电影中文网址| svipshipincom国产片| aaaaa片日本免费| 国内毛片毛片毛片毛片毛片| 国产免费现黄频在线看| 1024香蕉在线观看| 好男人电影高清在线观看| 岛国在线观看网站| 亚洲成人免费电影在线观看| 老熟妇乱子伦视频在线观看| 丰满少妇做爰视频| 人人澡人人妻人| 国产成人精品无人区| 自线自在国产av| 日韩中文字幕欧美一区二区| 天堂俺去俺来也www色官网| 欧美黄色淫秽网站| 伦理电影免费视频| 午夜精品久久久久久毛片777| 久久久久久久久久久久大奶| 91麻豆精品激情在线观看国产 | 欧美精品高潮呻吟av久久| 国产亚洲精品第一综合不卡| 亚洲av片天天在线观看| 少妇猛男粗大的猛烈进出视频| 欧美国产精品一级二级三级| 亚洲免费av在线视频| 日本撒尿小便嘘嘘汇集6| 伦理电影免费视频| 亚洲美女黄片视频| 国产区一区二久久| 日日爽夜夜爽网站| 一级毛片精品| 黄色成人免费大全| 亚洲美女黄片视频| 久久精品人人爽人人爽视色| 成年人黄色毛片网站| 日韩大片免费观看网站| 日本av手机在线免费观看| 久久久久网色| 老汉色∧v一级毛片| 精品人妻1区二区| 国产精品免费一区二区三区在线 | 亚洲全国av大片| 在线观看一区二区三区激情| videosex国产| 久久九九热精品免费| 交换朋友夫妻互换小说| 桃花免费在线播放| 国产成人啪精品午夜网站| 国产成人精品久久二区二区免费| 动漫黄色视频在线观看| 国产av一区二区精品久久| 蜜桃在线观看..| 大型av网站在线播放| 最黄视频免费看| 色视频在线一区二区三区| 午夜精品国产一区二区电影| 天天影视国产精品| 国产福利在线免费观看视频| 操美女的视频在线观看| 淫妇啪啪啪对白视频| 国产视频一区二区在线看| 日韩中文字幕欧美一区二区| 午夜福利欧美成人| 日本黄色日本黄色录像| 亚洲国产欧美一区二区综合| 亚洲人成电影免费在线| 国产男女内射视频| 老司机福利观看| 女性生殖器流出的白浆| 99精品久久久久人妻精品| 国产精品一区二区在线不卡| 色视频在线一区二区三区| 97人妻天天添夜夜摸| 99九九在线精品视频| 成人精品一区二区免费| 日韩中文字幕欧美一区二区| 日韩 欧美 亚洲 中文字幕| kizo精华| 丝袜喷水一区| 99国产精品免费福利视频| 波多野结衣一区麻豆| 亚洲成人手机| 另类精品久久| 亚洲精品美女久久av网站| 91国产中文字幕| 国产精品欧美亚洲77777| 下体分泌物呈黄色| 人妻 亚洲 视频| 在线天堂中文资源库| av福利片在线| 日日摸夜夜添夜夜添小说| 亚洲全国av大片| 黄色视频,在线免费观看| e午夜精品久久久久久久| 丰满少妇做爰视频| 亚洲av成人不卡在线观看播放网| 亚洲av美国av| 18禁黄网站禁片午夜丰满| 美女福利国产在线| 日韩免费av在线播放| www.999成人在线观看| bbb黄色大片| 久久午夜综合久久蜜桃| 欧美日韩av久久| 国产野战对白在线观看| 午夜福利免费观看在线| 大香蕉久久成人网| 国产一区二区 视频在线| 水蜜桃什么品种好| 日韩中文字幕欧美一区二区| 国产精品麻豆人妻色哟哟久久| 国产91精品成人一区二区三区 | 久久国产精品大桥未久av| 桃红色精品国产亚洲av| 国产成人免费观看mmmm| 国产aⅴ精品一区二区三区波| 国产在视频线精品| 日韩免费av在线播放| 岛国毛片在线播放| 国产黄频视频在线观看| 精品午夜福利视频在线观看一区 | 99久久人妻综合| 欧美日韩亚洲综合一区二区三区_| 一级片免费观看大全| 久久精品国产99精品国产亚洲性色 | 午夜精品久久久久久毛片777| av福利片在线| xxxhd国产人妻xxx| 97人妻天天添夜夜摸| 亚洲专区中文字幕在线| 精品亚洲成国产av| 国产一卡二卡三卡精品| 男人操女人黄网站| 久久国产亚洲av麻豆专区| 欧美亚洲日本最大视频资源| 国产精品免费视频内射| 狠狠精品人妻久久久久久综合| 在线亚洲精品国产二区图片欧美| av天堂在线播放| 国产伦理片在线播放av一区| 黄色视频,在线免费观看| av视频免费观看在线观看| 成人特级黄色片久久久久久久 | 五月天丁香电影| 我的亚洲天堂| 欧美精品av麻豆av| 激情在线观看视频在线高清 | 美国免费a级毛片| 又紧又爽又黄一区二区| 天天躁狠狠躁夜夜躁狠狠躁| 天天影视国产精品| 日韩欧美三级三区| 国产av一区二区精品久久| 亚洲久久久国产精品| 国产在线视频一区二区| 视频区欧美日本亚洲| 性高湖久久久久久久久免费观看| 老司机影院毛片| 亚洲欧美一区二区三区黑人| 亚洲成国产人片在线观看| 老司机深夜福利视频在线观看| 精品亚洲成a人片在线观看| 男女之事视频高清在线观看| 欧美日韩亚洲国产一区二区在线观看 | 男女下面插进去视频免费观看| 一级片'在线观看视频| 亚洲专区国产一区二区| 高清视频免费观看一区二区| tocl精华| 久久国产亚洲av麻豆专区| 黄网站色视频无遮挡免费观看| 岛国毛片在线播放| 欧美人与性动交α欧美精品济南到| 午夜福利影视在线免费观看| 无遮挡黄片免费观看| 母亲3免费完整高清在线观看| svipshipincom国产片| 一二三四在线观看免费中文在| 精品亚洲成a人片在线观看| 最近最新免费中文字幕在线| av天堂久久9| 国产日韩欧美视频二区| 国产精品久久久人人做人人爽| 国产黄频视频在线观看| 天堂8中文在线网| 国产无遮挡羞羞视频在线观看| 日本黄色日本黄色录像| 无人区码免费观看不卡 | 久久精品人人爽人人爽视色| 日本av免费视频播放| 婷婷丁香在线五月| av电影中文网址| 热99久久久久精品小说推荐| 大香蕉久久成人网| 欧美一级毛片孕妇| 曰老女人黄片| avwww免费| 亚洲国产av影院在线观看| 欧美成人午夜精品| 男女床上黄色一级片免费看| 手机成人av网站| 人人妻,人人澡人人爽秒播| 精品久久久久久久毛片微露脸| 麻豆国产av国片精品| 黄色丝袜av网址大全| 夜夜骑夜夜射夜夜干| 天天躁夜夜躁狠狠躁躁| 香蕉丝袜av| 久久精品国产99精品国产亚洲性色 | 亚洲第一青青草原| 国产福利在线免费观看视频| 成人18禁在线播放| 久热爱精品视频在线9| 国产男靠女视频免费网站| 久久人妻福利社区极品人妻图片| 国产亚洲欧美精品永久| 变态另类成人亚洲欧美熟女 | 亚洲精品国产精品久久久不卡| 国产人伦9x9x在线观看| 久久精品国产亚洲av香蕉五月 | 岛国在线观看网站| 精品免费久久久久久久清纯 | 久久久久久久国产电影| 视频在线观看一区二区三区| 亚洲五月色婷婷综合| 免费观看av网站的网址| 在线观看免费高清a一片| 在线观看免费视频网站a站| 精品久久蜜臀av无| 99re在线观看精品视频| 亚洲精品国产一区二区精华液| 国产精品久久久av美女十八| 欧美一级毛片孕妇| 国产又色又爽无遮挡免费看| 国产精品秋霞免费鲁丝片| 精品一区二区三区视频在线观看免费 | 热99国产精品久久久久久7| 亚洲九九香蕉| 亚洲av片天天在线观看| 国产成人精品在线电影| 精品少妇内射三级| 午夜福利免费观看在线| netflix在线观看网站| 亚洲熟妇熟女久久| 日本欧美视频一区| 18在线观看网站| 久久国产精品影院| 嫁个100分男人电影在线观看| 国产精品一区二区精品视频观看| 夜夜爽天天搞| 亚洲欧美日韩另类电影网站| 国产av精品麻豆| 美女福利国产在线| 国产亚洲欧美在线一区二区| 欧美精品av麻豆av| 美女高潮喷水抽搐中文字幕| 久久ye,这里只有精品| 亚洲专区中文字幕在线| 亚洲国产av影院在线观看| 欧美+亚洲+日韩+国产| 国产精品美女特级片免费视频播放器 | 国产精品免费一区二区三区在线 | 啦啦啦视频在线资源免费观看| 天堂动漫精品| 啦啦啦视频在线资源免费观看| 欧美激情极品国产一区二区三区| 丝袜美足系列| 大陆偷拍与自拍| 国产亚洲精品一区二区www | 老司机午夜福利在线观看视频 | 无人区码免费观看不卡 | 国产激情久久老熟女| 国产又爽黄色视频| 在线观看免费午夜福利视频| 两性午夜刺激爽爽歪歪视频在线观看 | 日韩免费高清中文字幕av| 亚洲 国产 在线| 满18在线观看网站| 国精品久久久久久国模美| 国产真人三级小视频在线观看| 久久午夜亚洲精品久久| 成年人午夜在线观看视频| 日韩三级视频一区二区三区| 91麻豆av在线| 18禁观看日本| 久久中文字幕一级| 国产成人av激情在线播放| 国产在线观看jvid| av福利片在线| 变态另类成人亚洲欧美熟女 | 无人区码免费观看不卡 | 大香蕉久久成人网| 激情在线观看视频在线高清 | 久久精品成人免费网站| 一级毛片女人18水好多| 国产精品久久久av美女十八| 欧美中文综合在线视频| 黄片小视频在线播放| 丰满迷人的少妇在线观看| 搡老岳熟女国产| 一区在线观看完整版| 黄色视频在线播放观看不卡| 99国产综合亚洲精品| 国产aⅴ精品一区二区三区波| 欧美av亚洲av综合av国产av| 久久久精品国产亚洲av高清涩受| 成人特级黄色片久久久久久久 | 日本av免费视频播放| 激情在线观看视频在线高清 | 18禁美女被吸乳视频| 精品国产乱子伦一区二区三区| 日韩精品免费视频一区二区三区| 在线av久久热| 日韩欧美国产一区二区入口| 午夜精品久久久久久毛片777| 亚洲欧美精品综合一区二区三区| 天堂中文最新版在线下载| 91老司机精品| 18禁国产床啪视频网站| 精品欧美一区二区三区在线| 91成人精品电影| 超色免费av| 成年版毛片免费区| 国产福利在线免费观看视频| 亚洲自偷自拍图片 自拍| 国产精品.久久久| 天堂8中文在线网| 午夜久久久在线观看| 涩涩av久久男人的天堂| 欧美成人免费av一区二区三区 | 老汉色av国产亚洲站长工具| 亚洲中文日韩欧美视频| 狠狠婷婷综合久久久久久88av| 国产精品国产高清国产av | av一本久久久久| 精品少妇黑人巨大在线播放| 在线观看舔阴道视频| 性色av乱码一区二区三区2| 国产精品98久久久久久宅男小说| 亚洲天堂av无毛| 十八禁高潮呻吟视频| 一级片'在线观看视频| 国产欧美日韩一区二区三| 人人澡人人妻人| 国产精品九九99| 午夜福利影视在线免费观看| 一级,二级,三级黄色视频| 欧美激情久久久久久爽电影 | av线在线观看网站| 亚洲avbb在线观看| 欧美黄色片欧美黄色片| 91精品国产国语对白视频| 免费女性裸体啪啪无遮挡网站| 欧美国产精品va在线观看不卡| 纵有疾风起免费观看全集完整版| 999久久久精品免费观看国产| 亚洲人成电影免费在线| 免费黄频网站在线观看国产| √禁漫天堂资源中文www| 免费在线观看完整版高清| 美女福利国产在线| 久久久久久久精品吃奶| 国产精品久久久久久精品古装| 精品久久久久久电影网| 日韩欧美一区视频在线观看| 午夜免费鲁丝| 欧美日本中文国产一区发布| 午夜日韩欧美国产| 老鸭窝网址在线观看| 精品国产亚洲在线| 国产国语露脸激情在线看| 天天影视国产精品| 丰满迷人的少妇在线观看| 午夜福利免费观看在线| 欧美日韩av久久| 亚洲色图 男人天堂 中文字幕| 女人高潮潮喷娇喘18禁视频| 久久久久精品国产欧美久久久| 欧美激情极品国产一区二区三区| 国产精品久久久久成人av| 成在线人永久免费视频| 亚洲av片天天在线观看| 大陆偷拍与自拍| 咕卡用的链子| 精品国产一区二区三区久久久樱花| 王馨瑶露胸无遮挡在线观看| 99精国产麻豆久久婷婷| 美女福利国产在线| 99在线人妻在线中文字幕 | 大片免费播放器 马上看| 美女高潮喷水抽搐中文字幕| 丰满迷人的少妇在线观看| 另类亚洲欧美激情| 搡老乐熟女国产| 1024香蕉在线观看| 久久国产精品人妻蜜桃| a级片在线免费高清观看视频| 多毛熟女@视频| 国产深夜福利视频在线观看| 日韩一卡2卡3卡4卡2021年| 麻豆乱淫一区二区| 精品国产乱码久久久久久男人| 一区二区三区乱码不卡18| 在线天堂中文资源库| 黄色a级毛片大全视频| 亚洲,欧美精品.| 国产片内射在线| 蜜桃在线观看..| 国产精品 欧美亚洲| 欧美+亚洲+日韩+国产| 天天操日日干夜夜撸| 啦啦啦在线免费观看视频4| av视频免费观看在线观看| 久久久精品94久久精品| 欧美老熟妇乱子伦牲交| 丰满迷人的少妇在线观看| 两个人免费观看高清视频| 久久久久国产一级毛片高清牌| 国产精品麻豆人妻色哟哟久久| 亚洲国产欧美在线一区| 国产淫语在线视频| 自拍欧美九色日韩亚洲蝌蚪91| 久久久国产欧美日韩av| 久久精品成人免费网站| 国产日韩欧美视频二区| 香蕉国产在线看| 亚洲精品中文字幕在线视频| 精品国产乱码久久久久久小说| 丝袜喷水一区| 色尼玛亚洲综合影院| 国产欧美日韩一区二区三| 国产精品一区二区在线不卡| 涩涩av久久男人的天堂| 国产成人欧美在线观看 | 国产不卡一卡二| 亚洲avbb在线观看| 欧美变态另类bdsm刘玥| 男女高潮啪啪啪动态图| 国产av一区二区精品久久| 一级片'在线观看视频| av欧美777| 欧美另类亚洲清纯唯美| 丁香六月天网| 在线观看免费午夜福利视频| 美女视频免费永久观看网站| tocl精华| 蜜桃国产av成人99| 女警被强在线播放| 人人妻,人人澡人人爽秒播| 首页视频小说图片口味搜索| 国产精品久久久av美女十八| 国产精品免费大片| 美女高潮喷水抽搐中文字幕| 久久国产亚洲av麻豆专区| 大陆偷拍与自拍| 99国产精品一区二区三区| 亚洲国产精品一区二区三区在线| 成年版毛片免费区| 日韩大码丰满熟妇| 丰满饥渴人妻一区二区三| 99香蕉大伊视频| 午夜免费成人在线视频| 99热网站在线观看| 免费在线观看影片大全网站| 亚洲人成伊人成综合网2020| 在线观看人妻少妇| 国产精品免费大片| 免费在线观看完整版高清| 丁香六月欧美| 欧美日韩中文字幕国产精品一区二区三区 | 国产高清国产精品国产三级| 少妇裸体淫交视频免费看高清 | 久久久久久久久免费视频了| 亚洲av日韩在线播放| 每晚都被弄得嗷嗷叫到高潮| 午夜成年电影在线免费观看| 欧美日韩中文字幕国产精品一区二区三区 | 久久亚洲真实| av线在线观看网站| 国产一区二区 视频在线| 18禁美女被吸乳视频| 日韩免费av在线播放| 伦理电影免费视频| 免费观看av网站的网址| 老熟妇仑乱视频hdxx| 国产精品 欧美亚洲| 丝袜美腿诱惑在线| 丝袜在线中文字幕| 国产精品 欧美亚洲| 欧美日韩亚洲综合一区二区三区_| 久久婷婷成人综合色麻豆| 夫妻午夜视频| 午夜福利免费观看在线| 国产亚洲一区二区精品| 大片电影免费在线观看免费| 亚洲av电影在线进入| 国产无遮挡羞羞视频在线观看| 免费久久久久久久精品成人欧美视频|