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

    Assimilating All-sky Infrared Radiances from Himawari-8 Using the 3DVar Method for the Prediction of a Severe Storm over North China

    2021-04-07 10:21:20DongmeiXUZhiquanLIUShuiyongFANMinCHENandFeifeiSHEN
    Advances in Atmospheric Sciences 2021年4期

    Dongmei XU,Zhiquan LIU,Shuiyong FAN,Min CHEN,and Feifei SHEN

    1The Key Laboratory of Meteorological Disaster, Ministry of Education (KLME)/Joint International Research Laboratory of Climate and Environment Change (ILCEC)/Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters (CIC-FEMD), Nanjing University of Information Science &Technology, Nanjing 210044, China

    2Heavy Rain and Drought-Flood Disasters in Plateau and Basin Key Laboratory of Sichuan Province,Chengdu 610225, China

    3National Center for Atmos.Res., Boulder, Colorado, 80301, USA

    4Institute of Urban Meteorology, China Meteorological Administration, Beijing 100089, China

    ABSTRACT Although radar observations capture storm structures with high spatiotemporal resolutions,they are limited within the storm region after the precipitation formed.Geostationary satellites data cover the gaps in the radar network prior to the formation of the precipitation for the storms and their environment.The study explores the effects of assimilating the water vapor channel radiances from Himawari-8 data with Weather Research and Forecasting model data assimilation system(WRFDA) for a severe storm case over north China.A fast cloud detection scheme for Advanced Himawari imager (AHI)radiance is enhanced in the framework of the WRFDA system initially in this study.The bias corrections,the cloud detection for the clear-sky AHI radiance,and the observation error modeling for cloudy radiance are conducted before the data assimilation.All AHI radiance observations are fully applied without any quality control for all-sky AHI radiance data assimilation.Results show that the simulated all-sky AHI radiance fits the observations better by using the cloud dependent observation error model,further improving the cloud heights.The all-sky AHI radiance assimilation adjusts all types of hydrometeor variables,especially cloud water and precipitation snow.It is proven that assimilating all-sky AHI data improves hydrometeor specifications when verified against the radar reflectivity.Consequently,the assimilation of AHI observations under the all-sky condition has an overall improved impact on both the precipitation locations and intensity compared to the experiment with only conventional and AHI clear-sky radiance data.

    Key words:Himawari-8,data assimilation,all-sky,storm case

    1.Introduction

    Remote sensing observations are important in terms of improving the numerical weather prediction (NWP),especially for domains lacking conventional observations(McNally et al.,2000;Xu et al.,2016) through data assimilation (DA).It has been shown that advanced geostationary satellites are able to provide nearly continuous four-dimensional information about the atmosphere,similarly to radar observations.Geostationary satellites are able to monitor the severe weather with higher spatial and temporal resolutions.In recent years,new generation geostationary weather satellites were launched that are superior to their predecessors in terms of temporal and spatial resolutions for more spectral bands (Schmit et al.,2005;Bessho et al.,2016;Yang et al.,2017b).Most of them,such as the FY-4 series of the China Meteorological Administration (Yang et al.,2017b),Himawari-8/9 of the Japan Meteorological Agency(Bessho et al.,2016),and the Geostationary Operational Environmental Satellite (GOES)-R series of the National Aeronautics and Space Administration (NASA) (Schmit et al.,2005) have been recently investigated in both global and regional NWP models.

    In the global NWP models,the use of geostationary satellite observations is well demonstrated.For example,based on the four dimensional variational (4DVar) system,the European Centre for Medium-Range Weather Forecasts(ECMWF) assimilates Meteosat visible and infrared imager(MVIRI) onboard Meteosat-7 (K?pken et al.,2004) and the Spinning Enhanced Visible and InfraRed Imager (SEVIRI)on board Meteosat-8 (Schmetz et al.,2002).In the National Centers for Environmental Prediction (NCEP) global data assimilation system,the GOES imager observations are assimilated using three dimensional variational assimilation(3DVar).For regional NWP models,Stengel et al.(2009)explored the impact of SEVIRI assimilation with the 3DVar system (10 km) and 4DVar system (22 km) with the France/Aldin model respectively.Zou et al.(2011) and Qin et al.(2013) assimilated GOES imager with the NCEP Gridpoint Statistical Interpolation (GSI) 3DVar scheme using 9 and 10 km resolutions,respectively,for improving the track and intensity forecast for Hurricane Debby (2012) as well as the precipitation skills over the coastal areas of the Gulf of Mexico.Wang et al.(2018a) assimilated the layer precipitable water (LWP) retrieved from the Advanced Himawari-8 imager (AHI) for local storms.The impact of the cumulus parameterization schemes on assimilating the AHI LWP data are also studied in Lu et al.(2019).In the Weather Research and Forecasting Data Assimilation (WRFDA) system (Barker et al.,2012),Yang et al.(2017a) assimilated GOES imager radiances for convection-permitting (4 km)forecasts over Mexico.Wang et al.(2018b) investigated the added value of clear-sky AHI radiances assimilation on the analyses and forecasts for a severe storm using the WRFDA.

    The above studies prove that NWP can be improved by assimilating the Infrared (IR) Radiances data from geostationary satellites measurements that are sensitive to the mass and thermal state,such as temperature,pressure and wind.However,less appreciated is the benefit of satellite data measuring atmospheric moisture,clouds and precipitation (Geer et al.,2017,2019),since it is not possible to improve the model initial conditions by holding cloud constant in DA.It is also very desirable to fully utilize the cloud information from the IR channels,since the cloud development is observed earlier by the IR channels than precipitation by the radar observations.Thus,the problem of under-utilizing humidity and precipitation information from IR channels needs to be addressed (Chevallier et al.,2001;Chevallier and Kelly,2002) through all-sky IR radiance assimilation.Because of the uncertainties and complexities of the cloudy radiance simulation in the radiative transfer models and the analysis techniques with clouds considered,the all-sky assimilation of IR channels is among one of the most challenging issues in DA (Li et al.,2016b).Generally,the difficulties exist in the zero-gradient problems for the clear background,the complex statistics of the situation dependent background error,the observation errors and the simulation accuracy in the radiative transfer models.Recently,Minamide and Zhang (2017) investigated an adaptive observation error inflation method to limit the erroneous analysis increments due to large representativeness errors,as is often the case for cloud-affected radiances for GOES-16.Harnisch et al.(2016) investigated a unified method to assimilate the two humidity channels at 6.2 and 7.3 μm from Meteosat Second Generation (MSG) considering the nonlinear impact of clouds on both observation and model forecast as well as the systematic deficiency of model and observation operator.Based on a local ensemble transform Kalman filter(LETKF,Hunt et al.,2007) method,Okamoto (2017) evaluated the all-sky Himawari-8/AHI IR radiance assimilation in a mesoscale NWP model with 5-km resolution using LETKF after studying the observation error models.Okamoto et al.(2019) further compared the assimilation clear-sky and all-sky IR radiances for Himawari-8/AHI in a mesoscale system using a single humidity band.Improvements for temperature,water vapor,wind,and rainfall are found in their study because of the universal coverage of Himawari-8/AHI.They also proposed adaptive bias correction schemes by setting different predictors for all-sky IR radiance assimilation.Similar studies include that of Sawada et al.(2019),who assimilated Himawari-8/AHI radiances every 10-minutes to improve convective predictability.Zhang et al.(2018) assimilated all-sky IR radiances from GOES-16 ABI using EnKF for 1 km model convection-permitting severe thunderstorms prediction for a tornado case with an adaptive observation and first guess error inflation methods.They found that by removing spurious clouds,good estimation of ice particle and liquid particle fields could be analyzed.Honda et al.(2018a) conducted Himawari-8/AHI radiances assimilation every 10 minutes to improve convective predictability.The assimilation of allsky Himawari-8 AHI radiance for the moisture and hydrometeor initialization for TC predictions were also studied with promising results (e.g.,Honda et al.,2018b,2019;Minamide and Zhang,2018).

    This study aims to evaluate the effects of assimilating the clear-sky and all-sky IR radiances from Himawari-8 for a severe storm over north China in the framework of WRFDA using the 3DVar.The differences between the clear-sky and all-sky AHI radiance DA analyses and forecasts are explored to provide some guidelines for the usage of AHI radiance data with high spatial and temporal resolutions.

    This paper is organized as follows.The introduction of the WRFDA-3DVar system,the AHI IR radiance assimilation methods,and the background error modeling method are provided in section 2.Section 3 describes the experiment design and DA configurations.The results are illustrated in Section 4 before concluding in section 5.

    2.Methodologies

    2.1.3DVar method in WRFDA

    The cost function is the core of the 3DVar data assimilation technique.For the traditional 3DVar method of the WRFDA (Barker et al.,2012),the prescribed cost function is designed to achieve the optimal analysis through minimization iteratively as:

    where yis the observation vector,

    x

    and

    x

    represent the analysis and first guess (background) states,respectively,

    B

    is the background error covariance assigning the forecast error,the observation error is described by the observation error covariance matrix

    R

    ,and

    H

    stands for the observation operator mapping the variables in model space to the variables in observation space.

    2.2.AHI radiance Observation Operator,Quality Control and Cloud Detection

    The Community Radiative Transfer Model (CRTM:Liu and Weng,2006) is the forward operator applied in our DA experiments.CRTM is enabled to simulate both the clear-sky and all-sky radiance.The temperature and moisture profiles from the model are the essential inputs for the clear-sky radiance simulation,while extra hydrometeor inputs (such as the cloud rainwater mixing ratio) are also needed for the all-sky radiance simulation.To conduct satellite DA,the corresponding tangent linear model and adjoint model are required in addition to the forward operator.CRTM has proven to be quite efficient with the Jacobian model included for calculating the Jacobians.

    As is well known,it is necessary to correct the systematic biases before DA (Dee and Uppala,2009) by modifying the observation operator with the bias considered (Liu et al.,2012;Xu et al.,2013;Otkin et al.,2018).The modification includes four potentially state-dependent predictors(200—50 hPa,1000—300 hPa layer thicknesses,surface skin temperature,and total column water vapor) along with their coefficients,as well as a constant component of total bias.The variational bias correction [VarBC:(Auligné et al.,2007)]scheme is employed in this study.The bias information for the water vapor channels is estimated based upon the departure (observation minus simulated clear-sky radiance from the model background) information of the two week-long period from 0000 UTC 25 June 2018 to 0000 UTC 9 July 2018.Then,the bias correction is performed with VarBC,so that the correction coefficients of the predictors are updated along with the variational analysis during the minimization procedure.As is known,the OMB (observation minus background) errors under the all-sky condition are from multiple unconfirmed sources (Sawada et al.,2019).It is rather complicated to consider and classify the error sources specifically for all sky pixels.Thus,the bias correction coefficients applied for the all-sky AHI IR radiance assimilation in this study are only estimated based on the clear sky pixels excluding the cloud sky pixels.Therefore,it is reliable to directly ingest the bias correction coefficients under clear-sky conditions to represent those under all-sky conditions.

    For the all-sky AHI IR radiance DA in our current study,no quality control procedure is conducted,since the situation-dependent observation error model is applied.It should be pointed out that the observation error model may not be able to include all the error sources from the observations.Nevertheless,it is of interest to validate the observation error model by considering all kinds of observations without any quality control procedure.For the clear-sky AHI IR radiance DA,strict quality control procedures are required since the retention of erroneous AHI IR radiance data after quality control is likely to ruin the analyses.Erroneous data after quality control tend to further affect the design of the bias correction procedures.Clear-sky AHI IR radiance DA uses a range of quality controls,including:

    (1) removing pixels where satellite zenith angels are larger than 65°;

    (2) removing pixels where channels are over mixed surface types;

    (3) removing pixels where the standard deviation of the neighboring 9 observation points for channel 13 (10.4 μm)is larger than 2 K;

    (4) removing pixels with a gross check,namely |OMB|≥8 K;

    (5) eliminating data if bias-corrected innovation is larger than 3 times the observation error under clear-sky conditions.The clear sky observation errors will be described in section 2.3;

    (6) eliminating cloudy pixels with the cloud detection schemes that will be introduced as below.

    Clear-sky AHI IR radiance data assimilation relies on a cloud detection technique to separate the cloudy AHI IR radiances from the clear-sky AHI IR radiance before assimilation,allowing a better use of available data.The cloud-detection algorithm (named ZZ-cloud detection scheme hereafter) developed by Zhuge and Zou (2016) are applied to search different cloud cases to achieve a composite cloud mask (Wu et al.,2020).In this study,the ZZ-cloud detection scheme is enhanced in the framework of WRFDA system initially.With the 10 tests in the ZZ-cloud detection algorithm,the metrics provided by different tests indicate the signals for discriminating the existence of the cloud.

    2.3.AHI Radiance Observation Error Model

    Among the key issues of the all-sky AHI IR radiance data assimilation technology,the correct estimation of the observation errors has an important impact on the quality of the analysis because the observation error directly affects the determination of the weights for the observations in the analysis procedure.The observation errors should be assigned as situation-dependent,considering the representativeness errors that are smaller in clear-sky conditions and larger in the presence of clouds and precipitation (Geer and Bauer,2011).The observation errors should also be channel dependent since the detection characteristics and the extent of their sensitivities to various cloud and precipitation conditions are different.Based on the observation error model in Harnisch et al.(2016),the observation error related statistics are obtained by calculating the OMB difference from 0000 UTC 4 July 2018 to 1800 UTC 7 July 2018 every 6 hours by using the 3 h forecast as the background.The brightness temperature differences of the clear-sky and the all-sky AHI IR radiance simulated with the CRTM are first calculated from the background.The differences are then binned according to the brightness temperature simulated from the background under the all-sky condition.The clear-sky AHI IR radiance simulation from CRTM needs to remove all hydrometeors while keeping temperature and moisture profiles unchanged.The corresponding brightness temperature value is determined as the brightness temperature threshold (BT) of each channel,when the brightness temperature difference is less than a certain value (0.1 K).Once the BTof the three water vapor channels are determined,the averaged cloud amount

    C

    can be estimated based on both the background and the observation,indicated as:

    Here Obs is AHI IR radiance observation and Bis the simulated all-sky AHI IR radiance from the background.

    C

    and

    C

    are determined from the background and the observation respectively.The parallel relationships between the observation errors and the cloud amount are estimated based on a large amount of OMB samples (observation minus the simulated all-sky AHI IR radiance from the background).The four parameters — the minimum observation errors (SD),the maximum observation errors (SD),the minimum cloud amount

    C

    ,and the maximum cloud amount

    C

    — are determined from the statistics.With the four determined parameters,it is convenient to estimate the cloud amount dependent observation errors for any given

    C

    as follows,

    2.4.Modeling in WRFDA-3DVar

    The climatological background error statistics are computed offline with the commonly used National Meteorological Center (NMC) method (Parrish and Derber,1992).A number of 24 h and 12 h forecasts are launched at different times over roughly one month in July of the previous year in 2017 (Chen et al.,2018).Pairs of forecasts valid at the same time with different forecast lead times are extracted by the NMC method to represent the difference through averaging them over the time and domain,for example with 24 h and 12 h forecasts,respectively.Both 24 h and 12 h forecasts contain errors,while the truth is never known.Further,the

    B

    matrix is decomposed as

    B

    =

    UU

    with the preconditioning transform applied as

    Uv

    =

    x?x

    (Barker et al.,2004).Here

    x

    and

    x

    is the analysis and background respectively,while

    v

    represents the control variable.The Control Variable option 7 (CV7) is applied in this study to include the zonal wind component,the meridional wind component,the full temperature,the full surface pressure,and the pseudo-relative humidity as the control variables (Chen et al.,2016;Li et al.,2016a;Sun et al.,2016;Shen et al.,2019).In addition to the five commonly used control variables,another five hydrometeor-related control variables (snow:

    q

    ,rainwater:

    q

    ,cloud liquid:

    q

    ,ice:

    q

    ,and graupel:

    q

    mixing ratio) are also considered using GEN_BE v2.0 software developed by Descombes et al.(2014).These hydrometeor related control variables represent the microphysical variables for both ice and liquid particles in the cost function so as to effectively assimilate the cloudy AHI IR radiance observations.

    3.Experiment Setup

    3.1.AHI observation

    On 17 October 2014,the Meteorological Satellite Center of Japan Meteorological Agency launched the Himawari-8 satellite successfully into geosynchronous orbit.The AHI instrument onboard Himawari-8 has an observing capability significantly improved compared with its predecessor Himawari-7 series satellites,and is similar to the Advanced Baseline Imager (ABI) onboard U.S.GOESR satellites in terms of the hardware configurations (Schmit et al.,2005,2017).Full-disk imagery is provided every 10 minutes by AHI,covering Japan and China.There are 16 channels including IR channels (nadir resolution 2 km) and near-IR and visible channels (nadir resolution 0.5 or 1 km).Among them,channel 7 to channel 16 are IR channels with central wavelengths ranging from 3.9 μm to 13.3 μm.The specifications of the channel information appear in Bessho et al.(2016).In this study,the three water vapor channels of AHI IR radiance are assimilated,which are sensitive to middle-to-upper tropospheric humidity and less so to the lower troposphere moisture.The three water vapor channels include channel 8 (6.2 μm),channel 9 (6.9 μm),and channel 10 (7.3 μm).For clear-sky IR radiance assimilation,these three channels are applied in Wang et al.(2018b),For all-sky AHI IR radiance assimilation,only one humidity channel is assimilated in Minamide and Zhang (2018) and Honda et al.,(2018b) to exclude the potential correlations among humidity channels.However,the idea of assimilating two humidity channels are stated in Okamoto (2017)with 6.2 μm and 6.9 μm and in Sawada et al.,(2019) with 6.2 μm and 7.3 μm,respectively,for mesoscale systems.In Okamoto (2019),radiances from three water vapor channels and one window channel were assimilated with promising results.In addition,a 9 km resolution thinning mesh is applied to reduce the impact of the potential correlations in the radiance observations,following Wang et al.(2018b).It should be emphasized that 9 km may not be large enough to avoid the inter-channels correlations,especially for the AHI all sky data assimilation.The study of the impact of the thinning mesh on the results from the clear sky and all sky AHI data assimilation is another interesting topic,which is beyond the scope of this study and is planned for future studies.

    3.2.Model Configuration and Experimental Design

    In this study,a severe storm event over north China that occurred during 6?10 July 2018 was selected.The plots for the meteorological factors and weather situations are firstly examined before discussing the model configuration and the experimental design.There is an upper-level low pressure system over northwestern China at 500 hPa at 0000 UTC 6 July 2018 as the main factor for the persistent precipitation in north China.The upper-level low pressure system propagated to north China and became stronger on 8 July and 9 July 2018.Water vapor transport by the southwesterly flow near the subtropical high is strong at 850 hPa for the entire period,leading to a corresponding torrential rain(Fig.1).

    The numerical simulations are based on ARW(Advanced Research WRF) of WRFV4.0.3.The simulation area center over north China is determined as 37°N,105°E in Fig.2a with 649 × 500 grid points (9 km) for domain 1.For domain 2,the grid points are 550 × 424 (3 km).The domains are configured with vertical structure of 57 unequally spaced sigma (non-dimensional pressure) levels with 10 hPa as the cloud top.The initial conditions and boundary conditions for the numerical simulations are from ECMWF analyses with 0.25° × 0.25° resolutions.The following parameterizations are used:New Thompson microphysics scheme (Thompson et al.,2004),YSU planetary boundary layer scheme (Hong et al.,2006),NOAH land surface scheme (Tewari et al.,2004),New Tiedtke scheme cumulus scheme for the larger domain (D01) (Tiedtke,1989),and the Rapid Radiative Transfer Model (RRTM) longwave and shortwave radiation schemes (Mlawer et al.,1997).

    Fig.1.Wind field (vector;m s?1) and geopotential height (shaded;gpm) of the ECMWF analyses (ECA) at 500 hPa valid at (a) 0000 UTC 6 July,(b) 0000 UTC 8 July 2018,Wind field (vector;m s?1) and relative humidity (shaded;%) of the ECA at 850 hPa valid at (c) 0000 UTC 8 July 2018,(d) 0000 UTC 9 July 2018.

    Fig.2.(a) The WRF model domains of D01 and D02,(b) The distribution of the conventional observations valid at 0000 UTC 8 July 2018.The numbers of each observation are marked on the right.

    The cold-start and subsequent cycling data assimilation experiments are completed for the 4-day from 0000 UTC 6 July to 0000 UTC 10 July 2018.A 6 h spin-up run is initialized as the cold start from 1800 UTC to 0000 UTC for each day.Then,the data assimilation and forecast are conducted every 3 hours for the rest of the day.A 24 h deterministic forecast is initialized from the final analysis of each day.

    Two DA experiments are designed to investigate the impact of AHI IR radiances on the analyses and forecasts for the severe storm in north China.The first experiment was the benchmark experiment named as“conv+ahiclr”that assimilated the conventional observations in both the larger domain (D01) and the smaller domain (D02) and the AHI IR radiances under clear-sky conditions in D02.Conventional observations include aircraft reports,ship reports,sounding reports,satellite cloud wind data,GPS data,and ground stations provided by the China Meteorological Administration (CMA).A snapshot of the distribution of the conventional data is provided in Fig.2b.

    The second experiment named as“conv+ahicloud”is similar to the first experiment except for assimilating all-sky AHI IR radiance observations in D02.In D02,conventional observation and all-sky AHI IR radiance data are assimilated in two steps.The analysis with the assimilation of conventional observations from D02 is provided as the first guess for the all-sky AHI IR radiance data assimilation procedure.

    4.Results

    4.1.Radiance Data Assimilation Preprocessing

    4.1.1.

    Quality control for clear-sky AHI IR radiance assimilation

    Figure 3 provides the number of observations from channel 8 before and after the quality control procedure.It is seen that most of the AHI IR radiance data are located over land.With the cloud detection and other quality control procedures,about 1/4 of the total AHI IR radiance observations are kept in the DA system.The counts of the assimilated AHI IR radiance observations are more than 4000 for most of the data assimilation cycles.

    4.1.2.

    Bias correction for the clear-sky radiance assimilation

    Fig.3.The total observation counts before the quality control (black),over land before the quality control(green),over sea before the quality control (blue) and the total counts after the quality control (red) for channel 10.

    The bias correction coefficients of each predictor are estimated using OMB from 0000 UTC on 25 June 2018 to 0000 UTC on 9 July 2018 every 3 hours.The background is interpolated from the ECMWF analyses with 0.25° × 0.25°resolution.The statistically obtained bias correction coefficients are then used as the static coefficients for the first data assimilation cycle.For the clear-sky AHI IR radiance DA,the coefficients for the bias correction are updated by Variational Bias Correction scheme (VARBC) during the minimization,while those for the all-sky AHI IR radiance DA are kept constant.The bias is estimated with the obtained coefficients associated with the corresponding predictors for each DA cycle.Figure 4 shows time series of the 3 hourly OMB statistics with and without bias corrections for channel 9 of the AHI clear-sky radiance simulation from 0000 on 6 July 2018 to 0000 on 10 July 2018 by the conv+ahiclr experiment.Clearly,the bias correction scheme effectively reduces the standard deviation and the mean bias for almost all the analysis times,indicating the validity and rationality of the variational bias correction scheme.For the all-sky data assimilation,it is not straightforward to validate the bias correction scheme by the results of the all-sky bias correction,since the mean and standard deviation of the OMB in all-sky data assimilation include mixed information besides the systematic bias.

    4.1.3.

    Observation error of the cloudy radiance

    As was described in section 2.3,there is no significant difference for the radiative effect between the clear-sky and all-sky simulated radiance when the brightness temperature is greater than BT.The brightness temperature thresholds BTof channels 8,9,and 10 are 233.86 K,242.51 K,and 253.93 K respectively.As is well known,channel 10 is more affected by clouds,since the weighting function shows a maximum at a relatively low height.Thus,the brightness temperature threshold BTof channel 10 is higher than the other two channels.

    Figure 5 shows the OMB distributions normalized by the constant standard deviation of all the samples (blue curve).The statistical results are based on the OMB samples from 0300 UTC on 4 July 2018 to 2100 UTC on 7 July 2018 every 6 hours.The background field for calculating of OMB is the 3 h forecasts initialized from 0000 UTC on 4 July 2018 to 1800 UTC on 7 July 2018 every 6 hours.The distributions of the OMB normalized by the binned standard deviation with different ranges of cloud amount instead of using the constant standard deviation are also presented(green curve).The statistics are obtained from the same OMB samples in Fig.5 for the three water vapor channels.The red curves represent the Gaussian distribution.Note that the frequency of the negative OMB is significantly larger than that of the positive OMB for the blue line with the regular normalization.Following the dynamical normalization achieved by taking the varying cloud amount into account,the frequency of OMB in the negative region is quite comparable to that of the positive OMB,matching the Gaussian curves better.Approaching the Gaussian distribution for the OMB distribution is expected to lead more accurate analyses.

    4.2.Cloudy radiance simulation

    Figures 6a,6b,and 6c show the distribution of the observed AHI radiance along with the simulated radiance from the background and the analysis after the all-sky AHI IR radiance DA respectively,valid at 0000 UTC on 8 July 2018 for channel 10.The difference between the simulated radiance from the analysis and the background is also provided in Fig.6d.As previously mentioned,channel 10 of AHI radiance is rather sensitive to clouds and precipitation.It is noted that clouds mainly exist in the southeast portion of the domain.The brightness temperature simulated with the background is significantly higher than the observed brightness temperatures in most areas,indicating that the hydrometeors in the background state are generally less than the actual situation.Following the all-sky AHI IR radiance assimilation,the amount of cloud is effectively simulated by moistening the background for most of the domain(Figs 6c,6d).The results indicate that the assimilation of AHI water vapor radiance improves the analysis of temperature,humidity,and hydrometeor information,eventually leading to a better simulation of the AHI IR radiance.

    Fig.4.The mean bias and standard deviation of channel 9 for OMB with no bias correction (OMB-nb) and with bias correction (OMB-wb) from 0000 UTC 6 July 2018 to 0000 UTC 10 July 2018.

    Fig.5.The OMB distributions normalized by the standard deviation of all the samples (blue curves) or normalized by the binned standard deviation with different ranges of cloud amount (green curves),and the gaussian curves (red curves) for (a)channel 8 (wavelength:6.2 μm),(b) channel 9 (wavelength:6.9 μm) and (c) channel 10 (wavelength:7.3 μm).

    Figure 7 a shows the cloud height from the AHI level-2 cloud product provided by Japan Aerospace Exploration Agency (JAXA).Figures 7b and 7c present the diagnosed cloud height from the analyses in conv+ahiclr and conv+ahicloud respectively.The Unified Post Processor(UPP) software (Garcia-Reynoso and Mora-Ramirez,2017)is employed to diagnose the cloud height taking the water vapor,the cloud water,the cloud ice,and the precipitation snow into account from the NWP model.It is found that the cloud height is overestimated near (39°N,110°E) in conv+ahiclr.Similarly,the cloud amount is underestimated in the east of the model domain near the area of (40°N,116°E) from conv+ahiclr.Generally,the cloud amount and the cloud height in the analysis with the all-sky AHI IR radiance data assimilation match the AHI level-2 cloud products better.It indicates that the clouds are probably effectively analyzed by fully utilizing the hydrometeor information in both the radiative transfer model and the cost function for the all-sky AHI IR DA.

    4.3.Analysis

    Figure 8 shows the analysis increments of the water vapor (

    q

    ),cloud water (

    q

    ),cloud rain (

    q

    ),and precipitation snow (

    q

    ) water path from the conv+ahiclr and conv+ahicloud.Noticeable water vapor increments exist in Beijing and Hebei provinces.Additional outstanding increments of the water vapor are observed from conv+ahicloud northwest of Beijing close to the region of Hebei and Inner Mongolia with the all-sky DA of AHI IR radiance in Fig.8b.As expected,no visible hydrometeor increment is found in conv+ahiclr,since no hydrometeor is considered in the observation operator in the clear-sky AHI IR radiance data assimilation used in WRFDA.The

    q

    increment exists in the conv+ahicloud experiment,whereas the increment of

    q

    from conv+ahicloud is relatively smaller since precipitation is not observable by IR channels in most cases.The increment of

    q

    from conv+ahicloud is distinct,indicating the existence of snow-phase clouds from the analysis.This large increment of

    q

    explains the noticeable adjustment of the cloud height in the conv+ahicloud experiment in Fig.8,since

    q

    is one of the key inputs of UPP.

    Fig.6.(a) The observed brightness temperature,the simulated brightness temperature (b) from the background,(c)from the analysis of the all-sky radiance data assimilation,and (d) the difference between the simulated radiance from the analysis and the background (units:K) valid at 0000 UTC on 8 July 2018 for channel 10.

    The differences of water vapor flux (WVF) at 850 hPa between analyses of conv+ahicloud and conv+ahiclr are illustrated in Fig.9 along with the winds from the ECMWF analysis.Figure 9a shows the difference (conv+ahicloud minus conv+ahiclr) in the WVF valid at 0000 UTC on 8 July 2018.Compared with the conv+ahiclr,conv+ahicloud generally shows an increase of the WVF to the north and west of Hebei and Beijing,leading to overall intense model-simulated precipitation north of Beijing,especially at 0000 UTC on 8 July 2018.Water vapor can be transported from southwest to the rainfall area over Beijing area,since WVF represents the direction and magnitude of the water vapor transport to maintain a persistent heavy rainfall.It is consistent with Fig.9 showing existence of rainfall from south to north in Beijing.

    4.4.Verification against radar reflectivity and brightness temperature

    The absolute bias and RMSE of the 12 h forecasts against the radar reflectivity (RF) are shown in Fig.10 initialized from 1200 UTC on 6 July 2018 to 0000 UTC on 10 July 2018 every 12 hours.The radar RF observations are from 8 coastal radars in Handan of Hebei,Zhangbei of Hebei,Chengde of Hebei,Beijing,Tianjing,Shijiazhuang of Hebei,Qinhuangdao of Hebei,and Cangzhou of Hebei with total observations counts over 30000 (Fig.10a).Note that the results of applying the all-sky AHI IR radiance data assimilation are superior to those of using the clear-sky AHI IR radiance data assimilation when verified against the radar reflectivity in Fig.10 for both the bias and the RMSE.The RMSE is reduced for all of the forecasts while the improvements for the bias exist for most times.The reduced bias in Fig.10b indicates an improvement forecast by moisturizing the background after the all-sky AHI IR radiance data assimilation.

    Fig.7.Cloud height (units:gpm) for (a) AHI-level2,(b) conv+ahiclr,(c) conv+ahicloud 0000 UTC on 8 July 2018.

    Figure 11 provides the RMSE against the observed brightness temperature for the window channel 13.The RMSEs are calculated with different forecast lead times initialized from 0000 UTC on 6 July,0000 UTC on 7 July,0000 UTC on 8 July,and 0000 UTC on 9 July 2018 respectively.Consistent reduced RMSEs are found from the conv_ahicloud as expected.The reduction of RMSE persists for different forecast lead times although the largest reduction is observed at the analysis time.

    4.5.Rainfall forecasts

    Fig.8.(a,b) The water vapor,(c,d) the cloud water,and (e,f) the cloud rain,(g,h) the precipitation snow water path at 0000 UTC 9 July 2018 for conv+ahiclr (left panel),and conv+ahicloud (right panel).The locations of Beijing (BJ) and Heibei (HB) province are marked in Fig.8c.

    Fig.8.(Continued).

    Fig.9.The water vapor flux (shaded;g cm?1 hPa?1 s?1) difference between analyses of conv+ahicloud and conv+ahiclr at 850 hPa at (a) 0000 UTC 7 July 2018,(b) 0000 UTC 8 July 2018.The vectors show the direction and magnitude of the wind from ECMWF analysis.The locations of Beijing (BJ) and Heibei (HB)province are marked in Fig.9a.

    Fig.10.(a) The counts of radar observations,(b) bias,and (c) the RMSE against the radar reflectivity for the 12 h forecasts valid from 1200 UTC 6 July 2018 to 0000 UTC 10 July 2018.

    Figure12 presents the 24 h accumulated precipitation from 1200 UTC 8 July 2018 to 1200 UTC 9 July 2018 for the two DA experiments.The results initialized from other data assimilation cycles are not shown since similar results are found in terms of the accumulated precipitation.Generally,the precipitation pattern is more compact from conv+ahicloud than that from conv+ahiclear.As can be seen here,there is a rain band extending from the south to the north of Beijing marked as A in the observations.It is apparent that there is a rain band through Beijing from the conv+ahicloud (Fig.12c),at least to a certain extent,while for the clear-sky AHI radiance data assimilation experiment,there is no match for such a rain band.

    In order to quantitatively evaluate the performance of precipitation prediction,the Fractional Skill Scores (FSS,Roberts and Lean,2008) scores for the 1h accumulated rainfall calculated with a threshold of 5 mm are calculated as a function of forecast leading time in Fig.13a averaged from the forecasts initialized from 0000 UTC 7 July 2018 to 1200 UTC 9 July 2018 every 12 hours.The influence distance of the neighborhood used in this study is 9 km.The verification domain is a sub-domain of D02 shown in Fig.12 for the calculating of the FSS.It is found that the FSS decreasees with the forecast leading time.The FSS from conv+ahicloud are consistently higher than those from conv+ahiclr.In addition,the FSSs for the 24 h accumulated precipitation are calculated for the 5 mm,10 mm,15 mm,and 25 mm thresholds,respectively,averaging over the same period mentioned above.It can be seen from Fig.13b that for most of the precipitation thresholds,the scores from the DA experiment with all-sky AHI IR radiance are superior to those from the clear-sky AHI IR radiance DA experiment.The conv+ahicloud experiment yields higher scores than conv+ahiclr experiment does for the thresholds of 10 mm,15 mm,and 25 mm.

    Fig.11.The RMSE against the observed brightness temperature with different forecast leading time initialized from (a) 0000 UTC on 6 July 2018,(b) 0000 UTC on 7 July 2018,(c) 0000 UTC on 8 July 2018 (d) (a) 0000 UTC on 9 July 2018 for channel 13.

    Fig.12.24 h accumulated precipitation (units:mm) for (a) precipitation data from CMA,(b) conv+ahiclr experiment,and (c)conv+ahicloud experiment from 1200 UTC 8 July 2018 to 1200 UTC 9 July 2018.

    Fig.13.FSS (a) with a threshold of 5 mm h?1 as a function of forecast lead time (b) of the 24 h accumulated precipitation for conv+ahiclr and conv+ahicloud experiments averaged from 0000 UTC 7 July 2018 to 1200 UTC 9 July 2018 every 12 hour.

    5.Conclusion

    The effects of assimilating the infrared radiances from Himawari-8 data are studied with WRFDA for a severe storm over northern China.Both clear-sky and all-sky AHI IR radiance data assimilation (DA) experiments are designed.AHI radiance preprocessing is conducted including the bias corrections,cloud detection for clear-sky AHI IR radiance DA and observation error modeling for all-sky AHI IR radiance DA.AHI radiance data are fully applied without any quality control for all-sky AHI radiance data assimilation.

    With the cloud dependent observation error model,the analyzed all-sky AHI IR radiances fit the observations more closely,further improving the cloud heights.Compared to the clear-sky AHI IR radiance assimilation experiment,the all-sky AHI IR radiance assimilation is able to adjust all types of hydrometeor variables,especially for cloud water and precipitation snow.It is found that assimilating all-sky AHI IR radiance in addition to just conventional data is able to improve the hydrometeors when verified against the radar reflectivity.The assimilation of AHI data under the all-sky condition has an overall positive impact on both precipitation locations and precipitation intensity compared to the experiments with only conventional and clear-sky AHI IR radiance data.

    In future work,the statistical model of the background error covariance can be further enhanced to take multivariate correlations into account.In addition,the ensemble-variational DA method (Shen et al.,2015;Xu et al.,2016) is expected to provide better analyses through the use of the flow-dependent background error covariances.Four-dimensional ensemble-variational DA methods (Buehner et al.,2013;Liu and Xiao,2013;Lorenc et al.,2014;Kleist and Ide,2015;Shen et al.,2019) are planned to assimilate AHI IR radiance observations every 10 minutes for severe storm events.Other surface channels are also needed for investigation in order to perform careful treatment of the surface emissivity.How to effectively combine the radar observations and all-sky AHI IR radiance observations simultaneously is another interesting aspect in our future research.

    Acknowledgements

    .This research was primarily supported by the Chinese National Natural Science Foundation of China(G41805016,G41805070),the Chinese National Key R&D Program of China (2018YFC1506404,2018YFC1506603),the research project of Heavy Rain and Drought-Flood Disasters in Plateau and Basin Key Laboratory of Sichuan Province in China(SZKT201901,SZKT20 1904),the research project of the Institute of Atmospheric Environment,China Meteorological Administration,Shenyang in China (2020SYIAE02,2020SYIAE07).

    一区福利在线观看| 国产高清三级在线| 伊人久久精品亚洲午夜| 国产三级中文精品| 一二三四社区在线视频社区8| 国产三级黄色录像| 一区福利在线观看| 色哟哟哟哟哟哟| 国产在视频线在精品| 日本黄色片子视频| 99在线视频只有这里精品首页| 少妇人妻一区二区三区视频| 国产一级毛片七仙女欲春2| 神马国产精品三级电影在线观看| АⅤ资源中文在线天堂| 亚洲成av人片免费观看| 九色国产91popny在线| 国内精品一区二区在线观看| 亚洲精品一区av在线观看| 91麻豆精品激情在线观看国产| 日韩欧美免费精品| 国产精品亚洲av一区麻豆| 国产伦精品一区二区三区四那| 免费无遮挡裸体视频| 最好的美女福利视频网| 叶爱在线成人免费视频播放| 成年免费大片在线观看| 在线播放无遮挡| 亚洲 欧美 日韩 在线 免费| 色吧在线观看| 欧美xxxx黑人xx丫x性爽| 欧美成人免费av一区二区三区| 免费在线观看日本一区| 一区二区三区高清视频在线| 久久精品影院6| 成人无遮挡网站| 欧美极品一区二区三区四区| 国产成人a区在线观看| 免费搜索国产男女视频| 久久国产精品人妻蜜桃| 亚洲精品影视一区二区三区av| 国产主播在线观看一区二区| 国产精品综合久久久久久久免费| 伊人久久精品亚洲午夜| 91九色精品人成在线观看| 人人妻人人看人人澡| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 一个人看的www免费观看视频| 亚洲国产日韩欧美精品在线观看 | 国产高清视频在线观看网站| 久久这里只有精品中国| 69av精品久久久久久| 亚洲一区高清亚洲精品| 国内精品一区二区在线观看| 无遮挡黄片免费观看| 亚洲无线在线观看| 欧美色欧美亚洲另类二区| 国产av麻豆久久久久久久| 不卡一级毛片| 床上黄色一级片| 禁无遮挡网站| 精品国内亚洲2022精品成人| 大型黄色视频在线免费观看| 90打野战视频偷拍视频| 成年免费大片在线观看| 岛国在线观看网站| 国产精品久久久久久精品电影| 天堂√8在线中文| 亚洲一区二区三区不卡视频| 亚洲七黄色美女视频| 国产蜜桃级精品一区二区三区| 国产精品av视频在线免费观看| 女同久久另类99精品国产91| 国产精品 欧美亚洲| 一本一本综合久久| 天天躁日日操中文字幕| 精品人妻一区二区三区麻豆 | 久久亚洲精品不卡| 十八禁人妻一区二区| 久久伊人香网站| 国产黄片美女视频| 黄片小视频在线播放| 久久天躁狠狠躁夜夜2o2o| 亚洲精品456在线播放app | 色尼玛亚洲综合影院| h日本视频在线播放| 久久亚洲精品不卡| 十八禁人妻一区二区| 国产综合懂色| 悠悠久久av| 老汉色av国产亚洲站长工具| 操出白浆在线播放| 国产熟女xx| 婷婷精品国产亚洲av在线| 久久久久久久午夜电影| 欧美色欧美亚洲另类二区| 波野结衣二区三区在线 | 最近最新免费中文字幕在线| 色尼玛亚洲综合影院| 欧美一级a爱片免费观看看| 亚洲国产高清在线一区二区三| 欧美最新免费一区二区三区 | ponron亚洲| 19禁男女啪啪无遮挡网站| 麻豆成人午夜福利视频| 中文字幕av在线有码专区| 国产精品影院久久| 丝袜美腿在线中文| 夜夜夜夜夜久久久久| 中出人妻视频一区二区| 国产精品亚洲一级av第二区| 欧美最新免费一区二区三区 | 午夜精品一区二区三区免费看| 婷婷精品国产亚洲av| 中文字幕人成人乱码亚洲影| 制服丝袜大香蕉在线| 久99久视频精品免费| 中出人妻视频一区二区| 精品久久久久久久人妻蜜臀av| 亚洲中文字幕日韩| 亚洲精品成人久久久久久| 一二三四社区在线视频社区8| 久久精品91无色码中文字幕| 婷婷丁香在线五月| 真人做人爱边吃奶动态| 在线天堂最新版资源| 国产高清视频在线观看网站| 中文字幕精品亚洲无线码一区| 91在线观看av| 国产欧美日韩一区二区三| 1024手机看黄色片| 在线十欧美十亚洲十日本专区| 一区福利在线观看| 无人区码免费观看不卡| 国产乱人伦免费视频| 五月玫瑰六月丁香| 久久久国产成人精品二区| 成年版毛片免费区| 丰满乱子伦码专区| 国产精品亚洲av一区麻豆| 国产黄片美女视频| 久久精品影院6| 99久久精品国产亚洲精品| 亚洲欧美日韩东京热| 男女做爰动态图高潮gif福利片| 免费搜索国产男女视频| 欧美丝袜亚洲另类 | 俄罗斯特黄特色一大片| 国产精品一及| 国产高清有码在线观看视频| 一进一出抽搐动态| 国产成人av激情在线播放| 老熟妇乱子伦视频在线观看| 亚洲最大成人中文| 精品国产亚洲在线| 三级国产精品欧美在线观看| 在线视频色国产色| 亚洲av日韩精品久久久久久密| 久久久色成人| 精品人妻一区二区三区麻豆 | 久久久成人免费电影| 国产精华一区二区三区| 欧美绝顶高潮抽搐喷水| 国产精品av视频在线免费观看| 1024手机看黄色片| 成人国产综合亚洲| 亚洲精品美女久久久久99蜜臀| 中文字幕精品亚洲无线码一区| 欧美激情久久久久久爽电影| 免费电影在线观看免费观看| 国产成人欧美在线观看| 在线观看免费视频日本深夜| av在线天堂中文字幕| 丰满人妻熟妇乱又伦精品不卡| 免费在线观看影片大全网站| 亚洲欧美日韩卡通动漫| 午夜视频国产福利| av欧美777| 麻豆国产av国片精品| 亚洲国产日韩欧美精品在线观看 | 日日夜夜操网爽| 97人妻精品一区二区三区麻豆| 18美女黄网站色大片免费观看| 成人一区二区视频在线观看| 一区福利在线观看| 桃色一区二区三区在线观看| 动漫黄色视频在线观看| 51午夜福利影视在线观看| 最近最新中文字幕大全电影3| 老汉色∧v一级毛片| 97超视频在线观看视频| 天堂av国产一区二区熟女人妻| 国产乱人伦免费视频| 一进一出好大好爽视频| 久久精品国产清高在天天线| 久久香蕉国产精品| 久久久国产成人免费| 国产淫片久久久久久久久 | 国产三级黄色录像| 一级作爱视频免费观看| 淫妇啪啪啪对白视频| av福利片在线观看| 搡老岳熟女国产| 一夜夜www| 欧美3d第一页| 国产伦在线观看视频一区| 国产精品久久久久久亚洲av鲁大| 国产精品三级大全| 无人区码免费观看不卡| 中文字幕高清在线视频| 亚洲激情在线av| 中文字幕精品亚洲无线码一区| 黄色日韩在线| 在线观看美女被高潮喷水网站 | 成人一区二区视频在线观看| 亚洲美女视频黄频| 国产野战对白在线观看| 两个人视频免费观看高清| 午夜精品一区二区三区免费看| 精品久久久久久,| 高清毛片免费观看视频网站| 成人永久免费在线观看视频| 男女床上黄色一级片免费看| 欧美中文综合在线视频| 2021天堂中文幕一二区在线观| 中亚洲国语对白在线视频| 日日摸夜夜添夜夜添小说| av国产免费在线观看| 1024手机看黄色片| 精品福利观看| 淫秽高清视频在线观看| 天堂网av新在线| 亚洲精品成人久久久久久| 精品一区二区三区视频在线 | 男插女下体视频免费在线播放| 搡老熟女国产l中国老女人| 亚洲,欧美精品.| 国产97色在线日韩免费| 亚洲成人久久性| 91久久精品国产一区二区成人 | 一进一出好大好爽视频| 精品日产1卡2卡| 99久久久亚洲精品蜜臀av| 成人午夜高清在线视频| 99久久精品一区二区三区| 热99在线观看视频| 亚洲欧美激情综合另类| 免费在线观看影片大全网站| 嫁个100分男人电影在线观看| 欧美乱色亚洲激情| 淫秽高清视频在线观看| 免费人成在线观看视频色| 国产精品av视频在线免费观看| 国内精品久久久久久久电影| 欧美一区二区精品小视频在线| 亚洲国产精品成人综合色| 色综合婷婷激情| 黄片小视频在线播放| 国产v大片淫在线免费观看| 国产高潮美女av| 欧美在线一区亚洲| 一个人看视频在线观看www免费 | 婷婷精品国产亚洲av| 国产精品自产拍在线观看55亚洲| 看免费av毛片| 亚洲av美国av| 色尼玛亚洲综合影院| 日本与韩国留学比较| 一本一本综合久久| 亚洲成人免费电影在线观看| 中亚洲国语对白在线视频| 国产精品亚洲美女久久久| 长腿黑丝高跟| 真实男女啪啪啪动态图| 九九久久精品国产亚洲av麻豆| 国内精品久久久久久久电影| 亚洲 国产 在线| 国产精华一区二区三区| 在线观看免费视频日本深夜| 亚洲精品影视一区二区三区av| 好看av亚洲va欧美ⅴa在| 精品国产美女av久久久久小说| 久久九九热精品免费| 一卡2卡三卡四卡精品乱码亚洲| 国产高清有码在线观看视频| 99热6这里只有精品| 99热这里只有是精品50| 国产免费一级a男人的天堂| 琪琪午夜伦伦电影理论片6080| 一级黄片播放器| 熟妇人妻久久中文字幕3abv| 久久这里只有精品中国| 亚洲av免费高清在线观看| 亚洲中文日韩欧美视频| 人人妻人人澡欧美一区二区| 最近最新中文字幕大全免费视频| 国产久久久一区二区三区| 欧美日韩黄片免| 九九热线精品视视频播放| 国产精品久久久久久久久免 | 精品人妻一区二区三区麻豆 | 国产99白浆流出| 午夜福利免费观看在线| 免费看日本二区| 久久久久久人人人人人| 毛片女人毛片| 亚洲av二区三区四区| 深爱激情五月婷婷| 99久久99久久久精品蜜桃| 国产黄a三级三级三级人| 久久人妻av系列| 少妇丰满av| 国产97色在线日韩免费| 午夜影院日韩av| 三级毛片av免费| 亚洲成人久久性| 日本一二三区视频观看| 国产精品98久久久久久宅男小说| 欧美黑人巨大hd| 日本熟妇午夜| 日日干狠狠操夜夜爽| 99久久99久久久精品蜜桃| 国产一区二区在线av高清观看| 中文在线观看免费www的网站| 色噜噜av男人的天堂激情| 美女高潮的动态| 国产精品久久久久久亚洲av鲁大| 在线免费观看不下载黄p国产 | 成年版毛片免费区| 99久久无色码亚洲精品果冻| 午夜视频国产福利| av国产免费在线观看| 97人妻精品一区二区三区麻豆| 九九久久精品国产亚洲av麻豆| 亚洲五月天丁香| 中文在线观看免费www的网站| 色哟哟哟哟哟哟| 精品久久久久久久毛片微露脸| 97超级碰碰碰精品色视频在线观看| 国产免费男女视频| 久久99热这里只有精品18| 成年人黄色毛片网站| 搡老熟女国产l中国老女人| 国产蜜桃级精品一区二区三区| 午夜免费激情av| 欧美高清成人免费视频www| 欧美av亚洲av综合av国产av| 久久6这里有精品| 精品久久久久久久久久久久久| 每晚都被弄得嗷嗷叫到高潮| 国产精品女同一区二区软件 | 午夜免费成人在线视频| 观看美女的网站| 国产aⅴ精品一区二区三区波| 麻豆久久精品国产亚洲av| 成年女人永久免费观看视频| 国产精品嫩草影院av在线观看 | 亚洲欧美日韩卡通动漫| 激情在线观看视频在线高清| 久久亚洲真实| 亚洲国产精品合色在线| 男女那种视频在线观看| 18禁黄网站禁片免费观看直播| 99久久精品国产亚洲精品| 神马国产精品三级电影在线观看| 夜夜夜夜夜久久久久| 国产男靠女视频免费网站| 在线免费观看不下载黄p国产 | 久久久久国产精品人妻aⅴ院| 在线观看美女被高潮喷水网站 | 黄色成人免费大全| 午夜老司机福利剧场| av福利片在线观看| 男人的好看免费观看在线视频| 国产视频一区二区在线看| 亚洲精品久久国产高清桃花| 成年女人看的毛片在线观看| 国产欧美日韩一区二区精品| 欧美+亚洲+日韩+国产| 久99久视频精品免费| 日本免费一区二区三区高清不卡| 女警被强在线播放| 在线天堂最新版资源| 午夜亚洲福利在线播放| 精品午夜福利视频在线观看一区| 三级毛片av免费| 免费观看精品视频网站| 窝窝影院91人妻| 女人高潮潮喷娇喘18禁视频| 欧美激情在线99| 欧美最黄视频在线播放免费| 99在线人妻在线中文字幕| av专区在线播放| 成年女人永久免费观看视频| 欧美日本视频| 91在线精品国自产拍蜜月 | 深爱激情五月婷婷| 久久久久国产精品人妻aⅴ院| 欧美zozozo另类| 欧美不卡视频在线免费观看| 日韩欧美精品免费久久 | 欧美黄色片欧美黄色片| 免费人成在线观看视频色| 国产亚洲欧美在线一区二区| 中文资源天堂在线| 91在线精品国自产拍蜜月 | 一级作爱视频免费观看| 欧美又色又爽又黄视频| 亚洲一区高清亚洲精品| 午夜a级毛片| 99久久成人亚洲精品观看| 国产欧美日韩一区二区三| 真实男女啪啪啪动态图| 亚洲18禁久久av| 日韩欧美在线二视频| 亚洲国产精品999在线| 可以在线观看的亚洲视频| 男人舔奶头视频| 欧美黄色淫秽网站| av欧美777| 中文在线观看免费www的网站| 男女床上黄色一级片免费看| 中文字幕av在线有码专区| 免费大片18禁| 国产精品一区二区三区四区免费观看 | 1024手机看黄色片| 男人舔奶头视频| 国产成人欧美在线观看| 99久久精品热视频| 免费观看精品视频网站| 极品教师在线免费播放| 亚洲中文日韩欧美视频| 岛国视频午夜一区免费看| 蜜桃亚洲精品一区二区三区| 熟女电影av网| 亚洲精品国产精品久久久不卡| 观看免费一级毛片| 国产高潮美女av| 亚洲国产精品成人综合色| 内地一区二区视频在线| 午夜福利高清视频| 午夜精品在线福利| 国产高清视频在线播放一区| 久久久久性生活片| 啦啦啦免费观看视频1| 欧美国产日韩亚洲一区| 国产国拍精品亚洲av在线观看 | 国产精品久久久人人做人人爽| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 亚洲精品国产精品久久久不卡| netflix在线观看网站| 亚洲av电影在线进入| 免费人成在线观看视频色| 在线观看66精品国产| 亚洲国产精品999在线| av女优亚洲男人天堂| av天堂中文字幕网| 女人十人毛片免费观看3o分钟| 欧美绝顶高潮抽搐喷水| 国产一区二区亚洲精品在线观看| 制服人妻中文乱码| 国产午夜精品久久久久久一区二区三区 | 国模一区二区三区四区视频| 国产高潮美女av| 色视频www国产| 国产精品av视频在线免费观看| 99热这里只有是精品50| 乱人视频在线观看| 国产精品国产高清国产av| 桃红色精品国产亚洲av| 久久欧美精品欧美久久欧美| 熟女少妇亚洲综合色aaa.| 日韩欧美三级三区| 91在线观看av| 中亚洲国语对白在线视频| 我的老师免费观看完整版| av片东京热男人的天堂| 色综合亚洲欧美另类图片| 欧美最黄视频在线播放免费| 国产伦人伦偷精品视频| 午夜两性在线视频| bbb黄色大片| 午夜免费成人在线视频| aaaaa片日本免费| 中文字幕av成人在线电影| 午夜激情福利司机影院| 久久精品夜夜夜夜夜久久蜜豆| 黑人欧美特级aaaaaa片| a级毛片a级免费在线| 欧美日韩瑟瑟在线播放| 久久午夜亚洲精品久久| 日韩成人在线观看一区二区三区| 国产综合懂色| 亚洲性夜色夜夜综合| 成人国产综合亚洲| 亚洲国产色片| 欧美在线黄色| 亚洲在线自拍视频| 少妇裸体淫交视频免费看高清| 久久国产精品影院| 欧美日韩一级在线毛片| av国产免费在线观看| 两个人看的免费小视频| 国产一区二区三区在线臀色熟女| 波野结衣二区三区在线 | 国产老妇女一区| 国产麻豆成人av免费视频| 色噜噜av男人的天堂激情| 午夜精品一区二区三区免费看| 欧美一区二区亚洲| 麻豆成人av在线观看| 51国产日韩欧美| 男女下面进入的视频免费午夜| 怎么达到女性高潮| 性色av乱码一区二区三区2| 亚洲人成网站在线播放欧美日韩| 国内精品久久久久久久电影| 丰满的人妻完整版| 99久久99久久久精品蜜桃| 亚洲精品影视一区二区三区av| 精品电影一区二区在线| 亚洲av不卡在线观看| 亚洲av免费在线观看| 亚洲一区二区三区不卡视频| 日本一二三区视频观看| 亚洲自拍偷在线| 两个人看的免费小视频| 制服丝袜大香蕉在线| 亚洲电影在线观看av| 欧美丝袜亚洲另类 | 一个人免费在线观看的高清视频| 国产精品98久久久久久宅男小说| 51国产日韩欧美| 久久这里只有精品中国| 久久久久久久久中文| 亚洲五月天丁香| 男人舔奶头视频| 黄色丝袜av网址大全| 国产一区二区三区在线臀色熟女| 国产老妇女一区| 日本撒尿小便嘘嘘汇集6| 婷婷六月久久综合丁香| 色综合婷婷激情| 免费看十八禁软件| 精品人妻1区二区| 制服人妻中文乱码| 免费大片18禁| 日本黄大片高清| 91麻豆av在线| 亚洲内射少妇av| 亚洲av五月六月丁香网| 国产成人aa在线观看| 久久婷婷人人爽人人干人人爱| 99久久精品一区二区三区| 男女视频在线观看网站免费| 亚洲精品亚洲一区二区| 国产在视频线在精品| 极品教师在线免费播放| 久久精品91蜜桃| 老司机午夜福利在线观看视频| 国产精品一区二区三区四区久久| www日本在线高清视频| 欧美色视频一区免费| 男人舔女人下体高潮全视频| 小蜜桃在线观看免费完整版高清| 成人无遮挡网站| 国产野战对白在线观看| 午夜精品一区二区三区免费看| 最近最新中文字幕大全免费视频| 午夜福利高清视频| 高清在线国产一区| 舔av片在线| 动漫黄色视频在线观看| 精品日产1卡2卡| 欧美一区二区国产精品久久精品| 久久国产精品影院| 国产97色在线日韩免费| 亚洲精品美女久久久久99蜜臀| 少妇高潮的动态图| 成人av在线播放网站| 九九久久精品国产亚洲av麻豆| 亚洲国产日韩欧美精品在线观看 | 国产av在哪里看| 最新美女视频免费是黄的| 国产黄色小视频在线观看| 国产高潮美女av| 日本在线视频免费播放| 欧美一区二区亚洲| 免费在线观看影片大全网站| 国产久久久一区二区三区| 国产精品一区二区三区四区久久| 亚洲午夜理论影院| 免费人成视频x8x8入口观看| 精品电影一区二区在线| 亚洲天堂国产精品一区在线| а√天堂www在线а√下载| 国产精品综合久久久久久久免费| 欧美国产日韩亚洲一区| 免费大片18禁| 老熟妇仑乱视频hdxx| av黄色大香蕉| 精品免费久久久久久久清纯| 老熟妇仑乱视频hdxx| 午夜久久久久精精品| 91在线精品国自产拍蜜月 | 精品99又大又爽又粗少妇毛片 | 亚洲av免费在线观看| 国产精品一区二区三区四区免费观看 | 手机成人av网站| 国产日本99.免费观看| 亚洲av成人av| 网址你懂的国产日韩在线| 老司机午夜十八禁免费视频| 欧美黄色片欧美黄色片| 亚洲av五月六月丁香网|