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

    Evaluation of a Regional Ensemble Data Assimilation System for Typhoon Prediction※

    2022-12-07 10:27:40LiliLEIYangjinxiGEZheMinTANYiZHANGKekuanCHUXinQIUandQifengQIAN
    Advances in Atmospheric Sciences 2022年11期

    Lili LEI, Yangjinxi GE, Zhe-Min TAN*, Yi ZHANG, Kekuan CHU, Xin QIU, and Qifeng QIAN

    1Key Laboratory of Mesoscale Severe Weather/Ministry of Education, School of Atmospheric Sciences,Nanjing University, Nanjing 210063, China

    2National Meteorological Center, China Meteorological Administration, Beijing 100081, China

    ABSTRACT An ensemble Kalman filter (EnKF) combined with the Advanced Research Weather Research and Forecasting model(WRF) is cycled and evaluated for western North Pacific (WNP) typhoons of year 2016. Conventional in situ data, radiance observations, and tropical cyclone (TC) minimum sea level pressure (SLP) are assimilated every 6 h using an 80-member ensemble. For all TC categories, the 6-h ensemble priors from the WRF/EnKF system have an appropriate amount of variance for TC tracks but have insufficient variance for TC intensity. The 6-h ensemble priors from the WRF/EnKF system tend to overestimate the intensity for weak storms but underestimate the intensity for strong storms. The 5-d deterministic forecasts launched from the ensemble mean analyses of WRF/EnKF are compared to the NCEP and ECMWF operational control forecasts. Results show that the WRF/EnKF forecasts generally have larger track errors than the NCEP and ECMWF forecasts for all TC categories because the regional simulation cannot represent the large-scale environment better than the global simulation. The WRF/EnKF forecasts produce smaller intensity errors and biases than the NCEP and ECMWF forecasts for typhoons, but the opposite is true for tropical storms and severe tropical storms. The 5-d ensemble forecasts from the WRF/EnKF system for seven typhoon cases show appropriate variance for TC track and intensity with short forecast lead times but have insufficient spread with long forecast lead times. The WRF/EnKF system provides better ensemble forecasts and higher predictability for TC intensity than the NCEP and ECMWF ensemble forecasts.

    Key words:ensemble Kalman filter, typhoon prediction, ensemble forecast

    1.Introduction

    Over the last few decades, there has been a steady decrease in tropical cyclone (TC) track forecast errors; but there have been minimal changes in TC intensity forecast errors over the same period (e.g., Rogers et al., 2006; Rappaport et al., 2009). TC track and intensity forecasts from numerical weather prediction models are limited by initial condition errors that are associated with the large-scale environments and TC structures, as well as model errors introduced by grid resolution and physical parameterization schemes.While TC motion is mostly controlled by the large-scale environment, TC intensity depends on the large-scale kinematic and thermodynamic environment, the inner-core dynamics,and the lower boundary condition including the surface sea temperature, ocean heat content, and land surface (e.g.,Wang and Wu, 2004).

    One feature of TCs is the large gradients in the mass and wind fields, which are often difficult to solve due to relatively coarse grid resolutions compared to the scale of the TCs. Therefore, several different techniques have been developed to relocate the storm to the observed position and construct initial conditions with more realistic TC structures.One kind of technique splits the forecast field into basic and disturbance fields and simply relocates the TC circulation to the observed position, then obtains the updated forecast filed with a correct TC position (e.g., Liu et al., 2000; Hsiao et al., 2010). Another kind of technique generates a more realistic TC initial structure by adding a TC-like vortex into the analysis field through data assimilation schemes (e.g., Kurihara et al., 1995; Zou and Xiao, 2000; Wang et al., 2008).To avoid the interference between a bogus vortex and actual observations, there are techniques that first implant a bogus vortex and then assimilate the profile data into the bogus vortex (e.g., Chou and Wu, 2008). The third kind of technique is dynamical initialization that incorporates the three-dimensional initial TC structure through the numerical forecast model (e.g., Kurihara et al., 1993; Cha and Wang, 2013; Hendricks et al., 2013; Liu and Tan, 2016). These techniques are somewhat limited by their assumptions, so they may not be optimal at all times.

    Compared to the previously discussed TC initialization methods, advanced data assimilation methods, which are not specialized for TC initialization, have been shown to make great impacts on TC forecasts. Among the widely applied data assimilation methods, ensemble-based assimilation approaches, such as the ensemble Kalman filter (EnKF;Burgers et al., 1998), have shown great promise for atmospheric data assimilation in both global (e.g., Whitaker et al.,2008; Buehner et al., 2010a, b; Houtekamer et al., 2014)and regional applications (e.g., Dowell et al., 2004; Tong and Xue, 2005; Meng and Zhang 2008; Aksoy et al., 2009).As a Monte Carlo approximation to the Kalman filter(Kalman, 1960), the EnKF uses flow-dependent error statistics estimated from short-term ensemble forecasts to determine the weight given to observations relative to model forecasts and spread observation information to model state variables. The usage of flow-dependent error statistics should provide more effective assimilation of observations near TCs,and the set of ensemble analyses can naturally provide initial conditions for TC ensemble forecasting.

    Previous work has shown that EnKF assimilation is able to provide dynamically consistent TC state estimation and improve TC track and intensity forecasts. Hamill et al.(2011) initialized the National Centers for Environmental Prediction (NCEP) Global Forecasting System (GFS) with EnKF analyses and obtained improved TC track forecasts compared to the operational ensemble data assimilation systems at the time (e.g., ensemble transform technique (Wei et al., 2008) and ensemble transform Kalman filter (Hunt et al.,2007)). Torn and Hakim (2009) cycled an EnKF over the lifetime of Hurricane Katrina (2005) and obtained a 50% reduction in TC track and intensity errors compared to the NCEP GFS and National Hurricane Center (NHC) official forecasts. Moreover, Zhang et al. (2009, 2011b) showed that assimilating Doppler radar data from both land-based and reconnaissance aircraft platforms with an EnKF led to improved TC intensity forecasts. Assimilating microwave radiances with an EnKF was also found to be beneficial for TC predictions (Schwartz et al., 2012). Besides these individual case studies with relatively short periods over which observations are assimilated, Torn (2010) cycled an EnKF over the life cycle of multiple TCs ranging from marginal to intense TCs in the Atlantic basin and found that cycling with an EnKF system was particularly effective for weak TCs. Cavallo et al. (2013) evaluated a cycling EnKF for the 2009 North Atlantic hurricane season and obtained systematically reduced TC track and intensity errors by assimilating observations, except for strong TCs. Xue et al. (2013) systematically compared the EnKF and three-dimensional variational data assimilation (3DVAR; Kleist et al., 2009) for the 2010 North Atlantic hurricane season and found significantly improved TC intensity forecasts initialized from the EnKF compared to those initialized from the 3DVAR.

    Although previous studies have been encouraging, systematic work on regional cycling EnKF systems for the western North Pacific (WNP) basin is limited. Besides the internal dynamical and physical processes that have important impacts on TC track and intensity, TCs over the WNP are strongly influenced by complex environmental conditions,including vertical wind shear, synoptic-scale weather systems, Asia monsoons, easterly waves, Rossby wave energy dispersion, and so on (e.g., Fudeyasu and Yoshida, 2018;Ma et al., 2019). The interactions across scales impose challenges on the predictability, data assimilation, and forecasts for typhoons. Thus, this study describes the results of a cycling mesoscale EnKF system combined with the Advanced Research Weather Research and Forecasting model (WRF) applied for most of the 2016 WNP typhoon season. The WRF/EnKF system produces an 80-member ensemble analysis every 6 h. For each storm during the experimental period, a 5-d deterministic forecast is launched from the ensemble mean analysis every 6 h; and a 5-d ensemble forecast is obtained from the 80-member ensemble analyses for seven typhoons whose intensities are underestimated. The application of the combined data assimilation and regional forecasting system over a nearly complete typhoon season provides a large sample for assessing its potential benefit to TC forecasts.

    This paper is organized as follows. Section 2 describes the modeling and data assimilation system. An overview of the 21 storms studied for the 2016 WNP typhoon season is provided in section 3. Section 4 discusses verifications for the cycling data assimilation system, while section 5 provides verifications for the 5-d deterministic forecasts. The performances of ensemble forecasts are evaluated in section 6. A summary and conclusion are provided in section 7.

    2.Experimental design

    Ensemble analyses are generated every 6 h by cycling a WRF/EnKF system from 0000 UTC 1 July to 1200 UTC 21 October. During this period, 21 TCs are simulated, whose categories and starting and ending dates are given by Table 1 and tracks are shown in Fig. 1. Within the duration of each TC, a 5-d forecast is launched every 6 h from the ensemble mean analysis. Due to limited computational resources, 5-d ensemble forecasts are launched for the typhoons whose rapid intensification processes are not well captured by the WRF/EnKF system and the ensemble mean analysis of the minimum sea level pressure (SLP) is higher than the observed value. Table 2 lists the typhoons and initial times for the 5-d ensemble forecasts.

    Table 1. The name, category, and duration for each TC in the 2016 WNP typhoon season assimilation.

    Table 2. Forecast initialization times for the typhoons whose ensemble mean analysis of the minimal sea level pressure is higher than the observation.

    The simulation of the 2016 WNP typhoon season uses WRF V3.4 (Skamarock et al., 2008). The model’s domain 1,shown in Fig. 1, covers most of China and the WNP basin.This domain could include the track of all storms in the WNP basin. 18-km horizontal grid spacing is used for domain 1, with 520 × 660 grid points. When the Joint Typhoon Warning Center (JTWC) announces an advisory TC position, a vortex following domain 2 of 6-km horizontal grid spacing and 180 × 180 grid points is set for 6-hourly cycling. When a 5-day forecast is produced for both the ensemble mean and the ensemble members, an additional high-resolution domain is set into the vortex following domain 2. This high-resolution domain has 2-km horizontal grid spacing and 300 × 300 grid points. There are 57 vertical levels, with the model top at 10 hPa. The implementation of WRF has the following components: the WRF 6-class microphysics scheme (Hong and Lim, 2006), the Yonsei University(YSU) planetary boundary layer (PBL) scheme (Hong et al.,2006), the Noah land surface model (Ek et al., 2003), the Rapid Radiative Transfer Model (RRTM) longwave scheme(Mlawer et al., 1997), and the RRTM shortwave scheme(Iacono et al., 2008). The cumulus parameterization uses the Tiedtke cumulus scheme (Tiedtke, 1989; Zhang et al.,2011a) and is only used for the 18-km domain 1.

    Fig. 1. Tropical cyclone tracks for each of the 21 TCs studied here. See Table 1 for a detailed list of storms. The colors are used to differentiate the TC tracks.

    The ensemble initial conditions (ICs) at the starting date 0000 UTC 1 July are generated from the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) analysis of 0.25° resolution using the fixed-covariance perturbation technique of Torn et al. (2006). This technique produces random perturbations that sample the NCEP background error covariance by use of the WRFDA-3DVAR (Barker et al., 2012), and the initial ensembles are generated by adding these random perturbations to the GFS analysis. Ensemble lateral boundary conditions (LBCs) are generated in a similar manner to the ensemble ICs, except that the NCEP GFS forecasts valid at the appropriate times are used. LBCs for the 5-d forecast launched from the ensemble mean analysis are produced using the NCEP GFS forecasts launched from the same analysis date but at the appropriate times.

    Observations from the NCEP global data assimilation system (GDAS) are assimilated every 6 h, including conventional in situ observations, cloud-motion vectors, and remotely sensed satellite radiances from the Advanced Microwave Sounding Unit-A (AMSUA), High resolution IR Sounder (HIRS), Atmospheric IR sounder (AIRS), and Microwave Humidity Sounder (MHS). For the assimilated radiance observations, a thinning mesh of 60 km is used, at which the radiance observation errors are assumed to be uncorrelated (Lin et al., 2017). In addition, the JTWC advisory minimum SLP at the observed position (latitude and longitude of the lowest sea level pressure) are assimilated. The observation error variances are the same as those used in the NCEP GDAS.

    An ensemble square-root filter (EnSRF; Whitaker et al.,2008) that is very similar to the NOAA operational EnKF for the NCEP GFS in 2016 (https://dtcenter.org/sites/default/files/community-code/enkf/docs/users-guide/EnKF_UserGuide_v1.3.pdf) is used to assimilate the observations.Note that, recently, a localized ensemble transform Kalman filter with model space localization (Hunt et al., 2007; Lei et al., 2018) has been implemented for the NCEP GFS. Ensemble size is 80. The observation priorHxb, whereHis the observation forward operator andxbis the model ensemble background or prior, is computed by the “observer” portion of the Gridpoint Statistical Interpolation system (GSI; Wu et al., 2002; Kleist et al., 2009). The “observer” runs separately for the ensemble mean and the 80 ensemble members,and the obtained observation ensemble priors are used by the EnKF. Bias correction of radiance observations is computed using the EnKF based on iterated analysis error covariance (Miyoshi et al., 2010).

    The EnKF system requires additional steps designed to overcome sampling errors that result from using a limited ensemble size and also account for model error. Sample covariances derived from the ensembles are localized by the Gaspari and Cohn (1999) localization function that is an approximately Gaussian fifth-order piecewise continuous polynomial function. Observation space localization(Houtekamer and Mitchell, 1998) is applied, by which observation impact is tapered to 0 at 1000 km in the horizontal and 1.5 ln (hPa) in the vertical. The vertical locations of the non-local radiance observations are assigned to the vertical level at which the weighting function maximizes. At each assimilation time, the ensemble deviations from the ensemble mean are inflated posterior to assimilation using the relaxation posterior ensemble spread to prior ensemble spread (relaxation-to-prior spread; Whitaker and Hamill, 2012) with a relaxation coefficient of 1.15. To maintain an appropriate ensemble spread, a relaxation coefficient larger than 1.0 is necessary, which forces posterior ensemble spread to be larger than prior ensemble spread (Schwartz and Liu, 2014;Schwartz, 2016). During a 10-day test period, a group of sensitivity experiments with varying localization and inflation parameters are conducted. The localization and inflation parameters are chosen based on the sensitivity experiments that provide the smallest 6-h prior errors comparing to the conventional observations and JTWC advisory TC data.

    To overcome the underestimation of TC initial intensity given an 18-km horizontal resolution (Torn, 2010) and avoid the technical challenges associated with moving nested domains for each ensemble member (like each ensemble member having a different nest location and interpolation from the coarse grid), the non-cycled nesting procedure of Cavallo et al. (2013) is adopted. Given our domain configuration, the EnKF data assimilation is done on the 18-km resolution domain 1, but not on the 6-km vortex-following domain 2. However, the 18-km resolution domain 1 benefits from the increased resolution for storm intensity, since the 6-h forecast on the 18-km grid is averaged from the 6-km grid where the 6-km grid exists. Cavallo et al. (2013) noted a 22% error reduction in the 6-h prior of TC minimum SLP and a moderate error reduction in the 6-h prior of TC maximum wind speed when using non-cycled nests within the data assimilation compared to no nests within the data assimilation.

    3.Overview of cases

    A short review of the 21 TCs in the 2016 WNP typhoon season during the assimilation period is provided in this section (Fig. 1). Most of the TCs formed in the main development region, while a few formed at relatively western longitudes. More than half of the TCs moved westward, and some of them turned poleward at some point in their lives,while several TCs directly moved poleward. Besides these TC tracks, Typhoon Lionrock first moved southwest, then turned northeast, and at last turned towards the northwest due to a high pressure system located east of Japan.

    Fig. 2. (a) 6-h forecast ensemble mean RMS error (dark gray bar) and ensemble spread (light gray bar) of TC positions as a function of TC intensity. The number of verification times is given along the top. (b) 6-h forecast ensemble mean bias of TC positions for tropical storm (TS, blue line), severe tropical storm (STS, green line), typhoon (TY, red line), and all storms(ALL, black line). The range rings denote 10-km intervals.Error bars denote the 5% and 95% percentiles determined from bootstrap resampling.

    Fig. 3. (a) 6-h forecast ensemble mean RMS error (dark gray bar) and ensemble spread(light gray bar) of TC minimum SLP as a function of TC intensity. The number of verification times is given along the top. (b) 6-h forecast ensemble mean bias of TC minimum SLP as a function of TC intensity. (c) As in (a), but for the TC maximum wind speed. (d) As in (b), but for the TC maximum wind speed. Error bars denote the 5% and 95% percentiles determined from bootstrap resampling.

    Fig. 4. RMS error of TC positions for 5-d forecast as a function of forecast hour for (a) TS, (b) STS, (c) TY,and (d) ALL. The blue solid line denotes the WRF/EnKF forecast launched from ensemble mean analysis,and the red and black solid lines are the NCEP GFS forecast and ECMWF forecast, respectively. The number of verification times for WRF/EnKF and NCEP is given by the first row along the top, and the number of verification times for ECMWF is given by the second row along the top. Error bars denote the 5% and 95%percentiles determined from bootstrap resampling.

    The lifetime durations and categories of the 21 TCs are given in Table 1. The lifetime durations of the storms ranged from 24 h (Rai) to 11 days (Lionrock). There were short-lived storms, such as Rai, Nida, Lupit, and Aere, and long-lived storms, such as Sarika, Chaba, Malakas, and Lionrock. According to the Regional Specialized Meteorological Centre (RSMC) Tokyo’s tropical cyclone intensity scale(Typhoon Committee, 2015), there were five tropical storms(TS) (Lupit, Conson, Dianmu, Kompasu, and Rai), whose maximum wind speeds were less than 24.2 m s-1; and there were five severe tropical storms (STS) (Mirinae, Nida,Omais, Chanthu, and Aere), whose maximum wind speeds were less than 32.4 m s-1. The remaining storms reached typhoon (TY) scale, with maximum wind speeds larger than 32.4 m s-1. Based on the Kaplan-DeMaria criteria (Kaplan and DeMaria, 2003), all the typhoons were characterized by at least one instance of rapid intensification (maximum wind speed increasing more than 15.4 m s-1per 24 h). Both Malakas and Songda met the rapid intensification criteria twice during their lifetimes. A wide variety of TC intensities are captured in this study.

    4.Cycling verification

    In this section, the cycling WRF/EnKF system is evaluated. Ensemble priors, i.e., 6-h ensemble forecasts, are verified against the JTWC advisory TC track position and intensity estimates that are the same as the assimilated quantities.The simulated TC position is given by the location of the minimum SLP. Figure 2 shows the TC track RMS errors, ensemble spread, and biases (observation minus forecast) for tropical storms (TS), severe tropical storms (STS), typhoons(TY), and all TCs (ALL), respectively. Figure 3 represents the RMS errors and biases of TC minimum SLP and maximum wind speed for different TC categories. The statistical significances of RMS errors and biases are determined using a bootstrap resampling with replacement method(Efron and Tibshirani, 1993). Each error distribution is resampled 1000 times with a bootstrap process. Error bars denote 5% and 95% percentiles; thus, only significance levels greater than 90% are considered statistically significant.

    For all TC categories, the ensemble spread of TC position is slightly smaller than the RMS error. As shown by Anderson(1996) and Hamill (2001), verification of an ensemble forecast can be affected by the used imperfect observations. The ensemble spread here does not take into account the observation uncertainties; thus, the generally consistent TC position error with ensemble spread indicates a well calibrated ensemble being provided by the WRF/EnKF system (Houtekamer et al., 2005). The TC position error and spread of ALL are dominated by those of TY, because TY has many more samples than TS and STS. The position error and spread are inversely proportional to the TC intensity. This may be a result of poorly defined TC centers for weak TCs. Note that the 6-h track errors are larger than the operational forecast errors because there are ensemble members with large track errors, especially when the vortex is weak (figures are not shown). The zonal and meridional position biases of TY are also much smaller than those of TS and STS. For all instances, there are larger meridional position biases than zonal position biases, which is consistent with the persistently positivev-wind biases (figures are not shown).

    Unlike for TC position, the ensemble spread of minimum SLP is insufficient to present the ensemble mean error, especially for TY. Similarly, the magnitude of TC maximum wind speed suffers from a large mismatch between the ensemble mean error and ensemble spread. Thus, the ensemble spread is not calibrated well for TC minimum SLP and maximum wind speed. TY has a much larger minimum SLP error than TS and STS, but the ensemble spread does not increase much from TS/STS to TY. Compared to TS and STS, the larger minimum SLP error of TY is also accompanied by a larger bias. TS and STS have similar errors of approximately 7 hPa, and similar biases of approximately +/- 0.75 hPa; but TY has an increased error of 22.21 hPa and an increased bias of -13.12 hPa. Although TS has a similar minimum SLP error and absolute bias to STS, TS has a larger maximum wind speed error and bias than STS. The error (bias) of maximum wind speed is 12.59 m s-1(-9.06 m s-1) for TS, 7.37 m s-1(-3.24 m s-1) for STS, and 14.61 m s-1(9.5 m s-1) for TY. Therefore, given a 6-km horizontal grid spacing, the WRF model is still unable to resolve the large gradients of TY’s wind and mass fields.The WRF/EnKF system overestimates the maximum wind speed for weak storms like TS and STS, but underestimates the minimum SLP and maximum wind speed for strong storms like TY.

    5.Forecast verification

    In this section, WRF/EnKF TC track and intensity deterministic forecasts launched from the ensemble mean analyses every 6 h are evaluated against the JTWC advisory TC track position and intensity estimates. For comparison, errors of the NCEP Global Ensemble Forecast System (GEFS) forecast(Zhou et al., 2017) and European Centre for Medium-Range Weather Forecasts (ECMWF) Ensemble Prediction System(EPS) global forecast (details are available online at https://www.ecmwf.int/sites/default/files/elibrary/2012/14557-ecmwf-ensemble-prediction-system.pdf) are also computed at the equivalent time. The NCEP GEFS has one control member and 20 ensemble members, and the 21 total ensemble members are integrated with the GFS model at approximately 34-km horizontal grid spacing for the first 8-d forecasts.The ECMWF EPS contains one control member and 50 ensemble members, and the 51 total ensemble members have a grid spacing of about 32 km through the 10-d forecast. Since WRF/EnKF uses a single forecast from the ensemble mean analysis, deterministic forecasts from the control members of the NCEP GEFS and ECMWF EPS are used for comparison in this section. Please note that WRF/EnKF has much higher horizontal grid resolutions (18 km, 6 km,and 2 km) than NCEP and ECMWF for the vortexes.

    Instead of launching a deterministic forecast from the ensemble mean analysis, Torn (2010) and Cavallo et al.(2013) used a single member of the WRF/EnKF analysis ensemble for the deterministic forecast in order to avoid the overly smoothed depiction of the TC mass and wind fields which result from ensemble averaging. The single-member initial condition is chosen as the closest member to the ensemble mean analysis based on normalized latitude, longitude,and minimum SLP differences. To show the differences of a single forecast from the ensemble mean analysis or from a single ensemble member, a 5-d deterministic forecast is run from the ensemble mean analysis and a randomly chosen ensemble member every 6 h for 5 TS, 5 STS, and 8 TY.Results show that the deterministic forecast from the ensemble mean analysis and that from a single ensemble member have similar errors of minimum SLP and maximum wind speed through the 5-d forecast period (figures are not shown). The forecast from the ensemble mean analysis has track RMS errors similar to those from a single ensemble member within 48 h, and the former has smaller track RMS errors than the latter beyond 48 h, although the differences are not statistically significant at a 90% confidence interval.Therefore, the deterministic forecast from the WRF/EnKF ensemble mean analysis is used here and compared to the control forecast of NCEP GEFS and ECMWF EPS.

    Figure 4 shows the track RMS errors as a function of forecast hours for TS, STS, TY, and ALL, respectively. The statistical significance of RMS error is denoted by the error bars showing 5% and 95% percentiles. The number of verification times is denoted on the top of the figure. The sample sizes of WRF/EnKF and NCEP are larger than that of ECMWF because WRF/EnKF and NCEP GEFS forecasts are computed every 6 h while the ECMWF EPS forecast is computed every 12 h. For TS, STS, TY, and ALL, WRF/EnKF deterministic forecasts generally have significantly larger track RMS errors than NCEP and ECMWF forecasts at short forecast lead times. The large initial track errors of WRF/EnKF are mainly resulted from the ensemble mean initial conditions, which have degraded initial position due to the ensemble averaging and have a subsequently fast error growth rate for the first 6 h due to imbalanced mass and wind fields of the ensemble mean. WRF/EnKF has insignificantly larger track errors or significantly smaller track errors than NCEP and ECMWF at long forecast lead times. Moreover, the meridional and zonal track biases of WRF/EnKF,NCEP, and ECMWF generally increase with forecast lead times, while the meridional track biases are larger than the according zonal track biases, especially at longer forecast lead times for TS and STS (Fig. 5). The meridional and zonal position error differences at longer forecast lead times are consistent with those from 6-h priors, which is possibly due to the persistently positivev-wind biases. Consistent with track RMS errors, WRF/EnKF deterministic forecasts in general have smaller track biases than NCEP and ECMWF at long forecast lead times. Thus, the regional simulation cannot better represent the large-scale environment compared to the global simulation, especially at short forecast lead times.

    Fig. 5. Biases of TC positions from 5-d forecasts for (a) TS, (b) STS, and (c) TY. The blue dot, green plus,and red square denote WRF/EnKF forecasts launched from the ensemble mean analyses, NCEP GEFS control forecasts, and ECMWF EPS control forecasts, respectively. The range rings denote 200-km intervals.

    Fig. 6. Same as Fig. 4, except for RMS error of TC minimum SLP.

    Fig. 7. Same as Fig. 4, except for bias of TC minimum SLP.

    Figures 6 and 7 display RMS errors and biases of minimum SLP as a function of forecast hour for different categories, respectively. For TS, WRF/EnKF deterministic forecasts produce larger RMS errors of minimum SLP than NCEP and ECMWF until about 60 h, and the error differences between WRF/EnKF and NCEP/ECMWF are not statistically significant at long forecast lead times. WRF/EnKF and NCEP (ECMWF) forecasts have negative (positive) biases of minimum SLP at short forecast lead times, while WRF/EnKF and ECMWF (NCEP) forecasts have negative (positive) biases for forecast lead times longer than 48 h. For STS, NCEP and ECMWF forecasts obtain similar RMS errors of minimum SLP, which are smaller than those of WRF/EnKF. Consistently, NCEP and ECMWF forecasts have similar biases of minimum SLP, which have smaller magnitudes than the positive biases of WRF/EnKF. For TY,WRF/EnKF deterministic forecasts generally have significantly smaller RMS errors of minimum SLP than NCEP and ECMWF until forecast lead times of 42 h and 72 h,respectively. Consistent with the RMS errors, NCEP and ECMWF forecasts have similar negative biases of minimum SLP, which have much larger magnitudes than the biases of WRF/EnKF. The RMS errors and biases of minimum SLP for ALL are similar to those of TY since the samples from TY are larger than those from TS and STS and thus the features of TY dominate.

    The RMS errors and biases of maximum wind speed at 10 m as a function of forecast hour for different categories are shown in Figs. 8 and 9, respectively. For TS, WRF/EnKF deterministic forecasts have similar RMS errors of maximum wind speed to NCEP and ECMWF within 36 h but larger RMS errors than NCEP and ECMWF beyond 48 h.Similar biases of maximum wind speed are obtained for WRF/EnKF, NCEP, and ECMWF forecasts. For STS, WRF/EnKF deterministic forecasts have larger RMS errors of maximum wind speed than NCEP and ECMWF. Meanwhile,WRF/EnKF deterministic forecasts have negative biases of maximum wind speed, and the magnitudes are larger than those from NCEP and ECMWF. Different from the weak storms, WRF/EnKF deterministic forecasts produce significantly smaller RMS errors of maximum wind speed than NCEP and ECMWF for TY. The error differences among WRF/EnKF, NCEP, and ECMWF beyond 90 h are generally not statistically significant. Consistently, WRF/EnKF deterministic forecasts produce positive biases of maximum wind speed, and the magnitudes are much less than those of NCEP and ECMWF. Similar results to the comparisons for TY are obtained for ALL.

    Therefore, although WRF/EnKF deterministic forecasts generally have larger RMS errors of minimum SLP and maximum wind speed than NCEP and ECMWF for weak storms like TS and STS, WRF/EnKF deterministic forecasts produce smaller RMS errors of intensity than NCEP and ECMWF for strong storms like TY. For weak storms, WRF/EnKF deterministic forecasts often have positive biases of minimum SLP and negative biases of maximum wind speed, which indicates overestimation of the intensity. For strong storms,WRF/EnKF deterministic forecasts often have negative biases of minimum SLP and positive biases of maximum wind speed, which indicates underestimation of the intensity. The underestimation of TY intensity is much better mitigated for WRF/EnKF compared to NCEP and ECMWF.There are environmental factors that may have influences on the vortex intensity, including temperature, specific humidity, vertical wind shear, etc. Figure 10 shows the profiles of differences for the mean specific humidity between WRF/EnKF and ECMWF forecasts. For each forecast, the mean specific humidity is averaged over an annulus that is centered around each vortex location with outer and inner circle radii of 5° and 2°, respectively. For weak storms, WRF/EnKF produces larger values of specific humidity below 500 hPa than ECMWF, which contributes to the intensity overestimation of weak vortexes for the WRF/EnKF. For strong storms,WRF/EnKF also has moister conditions than ECMWF,which explains the much better mitigated underestimation of TY intensity of WRF/EnKF compared to ECMWF.

    6.Ensemble forecast performance

    Performances of the ensemble forecasts for the typhoons listed in Table 2 are evaluated in this section. The mean absolute errors of the ensemble-mean position and intensity forecasts are compared with the mean ensemble spread.The consistency between these two quantities indicates that the ensemble contains an appropriate amount of variance,while the inconsistency between these two quantities indicates a lack of growing modes or reflecting model biases.For comparison, the mean absolute errors and ensemble spreads of NCEP GEFS and ECMWF EPS forecasts are also computed, as shown in Fig. 11.

    Fig. 8. Same as Fig. 4, except for RMS error of TC maximum wind speed.

    Fig. 10. Profiles of the differences of the mean specific humidity between WRF/EnKF and ECMWF 48-h forecasts for (a) STS and (b) TY. For each forecast, the mean specific humidity is the averaged specific humidity over an outer circle centered around each vortex with a 5° radius minus that over an inner circle with a 2° radius.

    Fig. 11. The mean absolute error (solid) and ensemble spread(dashed) from ensemble forecasts of typhoons listed in Table 2 as a function of forecast hour for (a) track, (b) minimum SLP,and (c) maximum wind speed. The blue lines denote the WRF/EnKF forecast, and the red and black lines displays the NCEP GEFS and ECMWF EPS forecasts, respectively.

    Fig. 12. 5-d ensemble forecast (a) tracks, (b) minimum SLP,and (c) maximum wind speed for typhoon Meranti from 0000 UTC 10 September. The thin blue, red, and green lines show the forecast values of the WRF/EnKF, NCEP, and ECMWF ensemble forecasts, respectively; and the thick lines denote the according ensemble mean. The black line denotes the observed value.

    Fig. 13. Same as Fig. 12, except for typhoon Sarika from 1200 UTC 12 October.

    For TC track, the mean absolute errors and ensemble spreads of WRF/EnKF, NCEP, and ECMWF ensemble forecasts are comparable within 30 h. Beyond that, the mean absolute errors of WRF/EnKF and NCEP, especially NCEP,increase at a higher rate than the ensemble spread; but the ensemble track forecasts of ECMWF have an even larger ensemble spread than the error. Thus, the ensemble track forecasts of WRF/EnKF and NCEP are underdispersive, while the ensemble track forecasts of ECMWF are overdispersive.For minimum SLP and maximum wind speed, WRF/EnKF,NCEP, and ECMWF ensemble forecasts contain appropriate variance when the forecast hour is less than 18 h. Beyond that, all three ensemble forecasts have spread deficiency.Compared to the ensemble forecasts of NCEP and ECMWF,WRF/EnKF ensemble forecasts have smaller intensity errors but contain larger ensemble spread. Therefore, WRF/EnKF ensemble forecasts lack variance for TC track (Puri et al., 2001; Magnusson et al., 2008; Torn, 2010), but WRF/EnKF provides better intensity ensemble forecasts than NCEP and ECMWF. Numerous factors have impacts on the ensemble forecasts, like the ensemble initial conditions,numerical model and parameterization schemes for different physical processes, model error representations, etc. A systematic investigation for the ensemble forecasts of different ensemble systems will be reported in a future study.

    To illustrate the ensemble performance in detail, ensemble forecasts for typhoon Meranti from 0000 UTC 10 September and typhoon Sarika from 1200 UTC 12 October by the WRF/EnKF, NCEP, and ECMWF are shown in Figs. 12 and 13, respectively. The track ensemble forecasts for typhoon Meranti by WRF/EnKF, NCEP, and ECMWF all have the ensemble plumes centered on the observed value;thus, typhoon Meranti appears predictable. The ensemble forecasts of the minimum SLP and maximum wind speed by the NCEP and ECMWF show some degree of confidence, but they are far from the observed intensity. The intensity ensemble forecasts by WRF/EnKF are much closer to the observed values than the ensemble forecasts by NCEP and ECMWF; WRF/EnKF predicts a central pressure that is approximately 40 hPa lower than NCEP and ECMWF and a wind speed that is about 20 m s-1stronger than NCEP and ECMWF. Typhoon Sarika also appears predictable from the track ensemble forecasts, although WRF/EnKF shows less predictability than NCEP and ECMWF. Similar to typhoon Meranti, the intensity ensemble forecasts for typhoon Sarika by NCEP and ECMWF are far from the observed value, and all ensemble members fail to capture the re-intensification process. However, the intensity ensemble forecasts by WRF/EnKF have the ensemble plume around the observed intensity, and most ensemble members predict the re-intensification process. At forecast lead times from 12 h to 120 h,WRF/EnKF ensemble members for typhoons Meranti and Sarika have moister conditions than ECMWF at low levels(figures are not shown), which could explain the stronger vortexes predicted by WRF/EnKF. Therefore, consistent with previous error statistics based on the mean absolute error and ensemble spread, WRF/EnKF has better intensity predictability than NCEP and ECMWF.

    7.Conclusions

    This study describes the performance of a cycling WRF/EnKF system during most of the 2016 WNP typhoon season. Conventional in situ data, radiance observations, and TC minimum SLP are assimilated every 6 h using an 80-member ensemble. For the 21 storms during the experimental period, a 5-d deterministic forecast is launched from the ensemble mean analysis every 6 h within the duration of each storm; and a 5-d ensemble forecast is produced from the ensemble analyses for 7 typhoons whose intensities are underestimated. The forecast errors are compared to the ECMWF and NCEP operational models.

    For all TC categories, the 6-h ensemble prior estimates of TC position from the WRF/EnKF system contain an appropriate amount of variance, while the TC intensity estimates are variance-deficient for all intensities, especially for TY.The TC position RMS error and spread are inversely proportional to the TC intensity, which may be a result of poorly defined TC centers for weak TCs. All TC instances are characterized by larger meridional position bias than zonal position bias. Category TY has much larger minimum SLP errors and biases than categories TS and STS, which indicates that a 6-km horizontal grid spacing is still unable to resolve the large gradients of TC wind and mass fields. Maximum wind speed errors and biases indicate that the WRF/EnKF system tends to overestimate maximum wind speed for TS and STS,but underestimate maximum wind speed for TY.

    Compared to the NCEP and ECMWF operational control forecasts, the WRF/EnKF deterministic forecasts from the ensemble mean analyses often has larger TC track errors for all categories because the regional simulation cannot better represent the large-scale environment compared to the global simulation. A blending method that merges the analyses of global and regional models can be beneficial for TC track forecasting (Hsiao et al., 2015). The meridional track biases of WRF/EnKF, NCEP, and ECMWF are generally larger than the according zonal track biases for TS and STS,while WRF/EnKF and ECMWF produce much smaller zonal and meridional track biases than NCEP for TY. The WRF/EnKF deterministic forecasts exhibit smaller TC intensity errors for TY than the NCEP and ECMWF control forecasts, which is due to the higher grid resolution of the WRF/EnKF system; but the WRF/EnKF forecasts have larger TC intensity errors for TS and STS. The WRF/EnKF deterministic forecasts often have positive biases of minimum SLP and positive biases of maximum wind speed for weak storms, which means overestimation of the intensity. The NCEP and ECMWF control forecasts have negative biases of minimum SLP and positive biases of maximum wind speed for strong storms, which means underestimation of the intensity, but the WRF/EnKF deterministic forecasts produce smaller intensity biases than the NCEP and ECMWF control forecasts. Profiles of specific humidity differences between WRF/EnKF and ECMWF show that WRF/EnKF produces moister conditions than ECMWF for both weak and strong storms, which explains how WRF/EnKF overestimates intensity for weak storms and why WRF/EnKF has better mitigated underestimation of intensity for strong storms.

    The ensemble forecasts from the WRF/EnKF system contain appropriate variance for TC track and intensity with short forecast lead times. With long forecast lead times, the ensemble forecasts of WRF/EnKF and NCEP are underdispersive for TC track, while the ensemble forecasts of ECMWF are overdispersive for TC track. The WRF/EnKF ensemble forecasts have smaller intensity errors but larger ensemble spread than the NCEP and ECMWF ensemble forecasts;thus, the WRF/EnKF system provides better TC intensity ensemble forecasts than NCEP and ECMWF, in terms of the comparison between amount of ensemble spread and mean absolute error. Moreover, the ensemble forecasts of WRF/EnKF can better capture the detailed intensity evolution than those of NCEP and ECMWF; thus, WRF/EnKF shows better intensity predictability than NCEP and ECMWF.

    The large initial track errors of WRF/EnKF are possibly a result of fast error growth due to imbalances caused by data assimilation, which could be mitigated by appropriate initialization methods and will be reported in a separate study.To improve the WRF/EnKF ensemble forecasts with enlarged ensemble spread, advanced data assimilation that updates ensemble perturbations with hybrid background error covariance (Lei et al., 2021) and additive inflation that can represent model uncertainties (Whitaker and Hamill,2012) need be further studied. Moreover, the cycling ensembles and 5-d ensemble forecasts here provide a unique dataset for studying TC structure, dynamics, genesis, and predictability. Ensemble sensitivity analysis (e.g., Torn and Hakim, 2008; Lei and Hacker, 2015) using this output to understand the dynamical processes that limit the predictability of TC track, intensity, and structure will be presented in a future study. Data assimilation algorithms that can capture multiscale features of TCs and improve TC predictability will also be investigated.

    Acknowledgements. We thank the editor and three anonymous reviewers for their insightful and constructive comments and suggestions. This work is jointly sponsored by the National Key R&D Program of China through Grant No. 2017YFC1501603, and the National Natural Science Foundation of China through Grant Nos. 41675052 and 41775057.

    Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

    精品久久久噜噜| 亚洲av成人精品一二三区| 国产熟女欧美一区二区| 国产成人午夜福利电影在线观看| 91成人精品电影| 免费看日本二区| 国产美女午夜福利| 一区二区三区免费毛片| 国产成人免费无遮挡视频| 亚洲情色 制服丝袜| 国产精品欧美亚洲77777| 大片电影免费在线观看免费| 最新中文字幕久久久久| 欧美日韩精品成人综合77777| 亚洲av国产av综合av卡| 日韩欧美 国产精品| 欧美成人午夜免费资源| 成年av动漫网址| 日韩亚洲欧美综合| 最后的刺客免费高清国语| 亚洲伊人久久精品综合| 国产精品女同一区二区软件| 国产中年淑女户外野战色| 中文在线观看免费www的网站| 欧美激情国产日韩精品一区| 蜜臀久久99精品久久宅男| 久久99蜜桃精品久久| 国产男人的电影天堂91| 插阴视频在线观看视频| 99热这里只有精品一区| 久久精品国产亚洲av涩爱| 日本黄大片高清| 欧美三级亚洲精品| 精品国产乱码久久久久久小说| av网站免费在线观看视频| 人人妻人人澡人人爽人人夜夜| av卡一久久| 国产成人精品婷婷| 97超视频在线观看视频| 欧美精品一区二区免费开放| 一区在线观看完整版| 日韩av在线免费看完整版不卡| 在线观看一区二区三区激情| 香蕉精品网在线| 国产黄色视频一区二区在线观看| 免费av中文字幕在线| 国产精品国产三级国产av玫瑰| 国内揄拍国产精品人妻在线| 高清av免费在线| 国产日韩一区二区三区精品不卡 | 综合色丁香网| 美女国产视频在线观看| 一级片'在线观看视频| av专区在线播放| 尾随美女入室| 国产在线免费精品| 精品亚洲乱码少妇综合久久| 高清毛片免费看| 亚洲美女搞黄在线观看| 男人狂女人下面高潮的视频| 少妇熟女欧美另类| 久久精品久久久久久久性| 久久韩国三级中文字幕| av一本久久久久| 男女边摸边吃奶| a级毛片在线看网站| 欧美日韩精品成人综合77777| 亚洲伊人久久精品综合| 少妇的逼好多水| 自拍欧美九色日韩亚洲蝌蚪91 | 日本av免费视频播放| 少妇高潮的动态图| 如日韩欧美国产精品一区二区三区 | 乱码一卡2卡4卡精品| 久久久久久久久久成人| 人妻系列 视频| 国产乱来视频区| 欧美一级a爱片免费观看看| 国产午夜精品一二区理论片| 国产一区二区三区综合在线观看 | 哪个播放器可以免费观看大片| 国产成人freesex在线| 色94色欧美一区二区| 国产精品.久久久| 97超碰精品成人国产| 免费大片黄手机在线观看| 欧美三级亚洲精品| 亚洲精品日韩在线中文字幕| 又大又黄又爽视频免费| 国产成人免费观看mmmm| 夫妻午夜视频| 久久人人爽av亚洲精品天堂| 精品少妇内射三级| 婷婷色综合www| 热99国产精品久久久久久7| 狂野欧美白嫩少妇大欣赏| 香蕉精品网在线| 久久狼人影院| 国产精品不卡视频一区二区| 一个人看视频在线观看www免费| 日本欧美视频一区| 国产成人午夜福利电影在线观看| 亚洲国产毛片av蜜桃av| 国产伦在线观看视频一区| 欧美变态另类bdsm刘玥| 国产在线视频一区二区| 制服丝袜香蕉在线| 中文欧美无线码| 国产黄色免费在线视频| 又大又黄又爽视频免费| 婷婷色综合大香蕉| 国产成人精品婷婷| av在线播放精品| 99久国产av精品国产电影| 亚洲图色成人| 午夜视频国产福利| 国产成人精品福利久久| 中文字幕人妻丝袜制服| 在线观看av片永久免费下载| 日日摸夜夜添夜夜添av毛片| 久久久久久伊人网av| 草草在线视频免费看| 久久99精品国语久久久| av免费在线看不卡| 日韩强制内射视频| 午夜老司机福利剧场| 亚洲av免费高清在线观看| 婷婷色av中文字幕| 国产免费又黄又爽又色| 91精品一卡2卡3卡4卡| h日本视频在线播放| 免费观看av网站的网址| 成人亚洲欧美一区二区av| 日韩av不卡免费在线播放| 久久久久精品久久久久真实原创| 美女主播在线视频| 一本—道久久a久久精品蜜桃钙片| 91久久精品国产一区二区三区| 六月丁香七月| 久久午夜福利片| 建设人人有责人人尽责人人享有的| 麻豆乱淫一区二区| 中文字幕久久专区| 熟女电影av网| 日本wwww免费看| av国产精品久久久久影院| 精品一区二区三区视频在线| 久久久久久久亚洲中文字幕| 国产伦理片在线播放av一区| 久久久午夜欧美精品| 欧美另类一区| 欧美日韩av久久| 亚洲第一av免费看| 人体艺术视频欧美日本| av免费观看日本| 观看美女的网站| 亚洲人与动物交配视频| 国产日韩欧美亚洲二区| 久久亚洲国产成人精品v| 在线观看av片永久免费下载| 精品久久久久久久久av| 大香蕉97超碰在线| 十八禁高潮呻吟视频 | 91aial.com中文字幕在线观看| 亚洲国产av新网站| 欧美日韩视频精品一区| xxx大片免费视频| 一区二区三区四区激情视频| 亚洲国产精品成人久久小说| 国产色爽女视频免费观看| 亚洲久久久国产精品| 亚洲欧美一区二区三区黑人 | 日韩精品免费视频一区二区三区 | 麻豆成人av视频| 久久女婷五月综合色啪小说| av在线观看视频网站免费| 国产一区二区在线观看av| 91久久精品国产一区二区三区| videossex国产| 欧美xxxx性猛交bbbb| 高清在线视频一区二区三区| 少妇被粗大的猛进出69影院 | 国产成人精品一,二区| 内射极品少妇av片p| 插逼视频在线观看| 中文天堂在线官网| av免费观看日本| 一区二区三区精品91| 99精国产麻豆久久婷婷| 成人毛片a级毛片在线播放| 中文欧美无线码| 午夜av观看不卡| 天天躁夜夜躁狠狠久久av| 视频区图区小说| 18禁裸乳无遮挡动漫免费视频| 免费观看性生交大片5| videossex国产| 精品少妇久久久久久888优播| 91精品国产国语对白视频| 免费观看的影片在线观看| 免费在线观看成人毛片| 亚洲av在线观看美女高潮| 午夜视频国产福利| 久久婷婷青草| 亚洲精华国产精华液的使用体验| kizo精华| 九九在线视频观看精品| 2022亚洲国产成人精品| 久久久久视频综合| 亚洲av电影在线观看一区二区三区| 午夜福利影视在线免费观看| 大香蕉97超碰在线| 亚洲欧洲日产国产| √禁漫天堂资源中文www| 久久青草综合色| 老司机影院成人| 另类精品久久| 国产日韩欧美在线精品| 午夜福利,免费看| 黄色一级大片看看| 亚洲精品视频女| 有码 亚洲区| 大片电影免费在线观看免费| 青青草视频在线视频观看| 亚洲国产精品国产精品| av专区在线播放| 人人妻人人澡人人爽人人夜夜| 天天躁夜夜躁狠狠久久av| 69精品国产乱码久久久| 精品一区在线观看国产| 国产成人免费无遮挡视频| 97精品久久久久久久久久精品| 久久精品国产a三级三级三级| 赤兔流量卡办理| av天堂中文字幕网| 国产成人精品无人区| 美女主播在线视频| 免费观看性生交大片5| 女的被弄到高潮叫床怎么办| 国产精品.久久久| 五月伊人婷婷丁香| 在线观看国产h片| 麻豆成人午夜福利视频| 欧美一级a爱片免费观看看| av.在线天堂| 亚洲国产精品专区欧美| 黄色毛片三级朝国网站 | 人妻一区二区av| 男人舔奶头视频| 亚洲内射少妇av| 久久精品夜色国产| 两个人免费观看高清视频 | 美女内射精品一级片tv| 久久久久久人妻| 色网站视频免费| 18禁裸乳无遮挡动漫免费视频| 午夜福利影视在线免费观看| 亚洲av在线观看美女高潮| 国产av码专区亚洲av| 婷婷色麻豆天堂久久| 久久97久久精品| 欧美成人午夜免费资源| kizo精华| 这个男人来自地球电影免费观看 | 久久亚洲国产成人精品v| 国产黄色免费在线视频| 久久ye,这里只有精品| 乱系列少妇在线播放| 国产淫语在线视频| 国产永久视频网站| 爱豆传媒免费全集在线观看| 好男人视频免费观看在线| 高清在线视频一区二区三区| 五月开心婷婷网| 久久国产乱子免费精品| 午夜免费鲁丝| 简卡轻食公司| 国产成人精品一,二区| 欧美精品人与动牲交sv欧美| 久久ye,这里只有精品| 妹子高潮喷水视频| 免费看光身美女| 久久精品久久久久久噜噜老黄| 在线观看三级黄色| 大香蕉97超碰在线| 日韩视频在线欧美| 在线天堂最新版资源| av又黄又爽大尺度在线免费看| 春色校园在线视频观看| av在线老鸭窝| 91久久精品电影网| 9色porny在线观看| 亚洲国产成人一精品久久久| 女性生殖器流出的白浆| 亚洲综合色惰| 熟女人妻精品中文字幕| 国产av一区二区精品久久| 国产永久视频网站| 有码 亚洲区| 韩国av在线不卡| 日韩精品免费视频一区二区三区 | 一级毛片 在线播放| 五月天丁香电影| 男女无遮挡免费网站观看| 成人毛片60女人毛片免费| 夫妻午夜视频| 亚洲欧美日韩另类电影网站| 国产伦理片在线播放av一区| 观看av在线不卡| 18禁在线无遮挡免费观看视频| 国产精品久久久久久精品古装| 性高湖久久久久久久久免费观看| 色婷婷av一区二区三区视频| 简卡轻食公司| 亚洲欧美成人精品一区二区| 国内揄拍国产精品人妻在线| 亚洲婷婷狠狠爱综合网| 中文字幕久久专区| 街头女战士在线观看网站| 精品久久国产蜜桃| 精品久久久久久久久av| 麻豆乱淫一区二区| 精品人妻熟女毛片av久久网站| 免费不卡的大黄色大毛片视频在线观看| 水蜜桃什么品种好| 亚洲久久久国产精品| 日韩大片免费观看网站| 午夜av观看不卡| 国产黄片视频在线免费观看| 新久久久久国产一级毛片| 国产精品一区二区性色av| 国产日韩一区二区三区精品不卡 | 免费观看无遮挡的男女| 亚洲国产毛片av蜜桃av| 欧美亚洲 丝袜 人妻 在线| 插阴视频在线观看视频| 久久久精品免费免费高清| 日韩中字成人| 亚洲国产色片| 九草在线视频观看| 自拍欧美九色日韩亚洲蝌蚪91 | 国产精品国产三级专区第一集| 国产日韩一区二区三区精品不卡 | 国产老妇伦熟女老妇高清| 尾随美女入室| 国产伦精品一区二区三区视频9| 亚洲精品中文字幕在线视频 | 亚洲怡红院男人天堂| 一区二区三区乱码不卡18| 蜜桃在线观看..| 亚洲美女视频黄频| 曰老女人黄片| 国产成人freesex在线| 一本色道久久久久久精品综合| 简卡轻食公司| 五月开心婷婷网| 777米奇影视久久| 成人二区视频| 精品熟女少妇av免费看| 一级二级三级毛片免费看| 97在线人人人人妻| 男女啪啪激烈高潮av片| 国产精品久久久久久av不卡| 日韩电影二区| 一区在线观看完整版| 麻豆成人午夜福利视频| 精品少妇黑人巨大在线播放| 97超视频在线观看视频| 妹子高潮喷水视频| 国产精品久久久久久久电影| 激情五月婷婷亚洲| 亚洲美女搞黄在线观看| 男女边摸边吃奶| 一级毛片我不卡| 中国美白少妇内射xxxbb| www.色视频.com| 最近的中文字幕免费完整| 免费av不卡在线播放| 亚洲国产精品一区三区| 在线观看一区二区三区激情| av在线观看视频网站免费| 日韩免费高清中文字幕av| 亚洲在久久综合| 97超视频在线观看视频| 国产乱来视频区| 中文字幕人妻丝袜制服| 久久影院123| 欧美 日韩 精品 国产| 男的添女的下面高潮视频| 免费观看性生交大片5| 久久婷婷青草| av福利片在线| 熟女人妻精品中文字幕| 亚洲av二区三区四区| 99九九线精品视频在线观看视频| 色婷婷久久久亚洲欧美| 国产淫语在线视频| 又大又黄又爽视频免费| 国产精品99久久99久久久不卡 | 高清在线视频一区二区三区| 国产黄色免费在线视频| 国产亚洲精品久久久com| 国产在线视频一区二区| 国产视频首页在线观看| 免费黄频网站在线观看国产| 日韩中字成人| 国产中年淑女户外野战色| 国产成人精品无人区| 国精品久久久久久国模美| av视频免费观看在线观看| 国产精品秋霞免费鲁丝片| 国产 精品1| www.av在线官网国产| 黑人猛操日本美女一级片| 国产乱人偷精品视频| 黄片无遮挡物在线观看| 熟女人妻精品中文字幕| 欧美变态另类bdsm刘玥| 亚洲av.av天堂| 极品教师在线视频| 美女大奶头黄色视频| 性色avwww在线观看| 九草在线视频观看| 丰满迷人的少妇在线观看| 欧美97在线视频| 我的女老师完整版在线观看| 成人二区视频| 亚洲激情五月婷婷啪啪| 一本大道久久a久久精品| 国产在线免费精品| 黑人猛操日本美女一级片| 欧美精品高潮呻吟av久久| 亚洲欧美成人综合另类久久久| 麻豆乱淫一区二区| 夫妻性生交免费视频一级片| 观看美女的网站| 日韩精品免费视频一区二区三区 | 国产黄片视频在线免费观看| 国产在线男女| 国产精品免费大片| 老女人水多毛片| 色婷婷av一区二区三区视频| 亚洲av不卡在线观看| 国产熟女午夜一区二区三区 | 精品亚洲成a人片在线观看| 久久99蜜桃精品久久| 国产精品99久久99久久久不卡 | 免费看光身美女| 久久人人爽av亚洲精品天堂| 中国美白少妇内射xxxbb| 国产精品欧美亚洲77777| 在线亚洲精品国产二区图片欧美 | 人妻夜夜爽99麻豆av| 国语对白做爰xxxⅹ性视频网站| 国产欧美亚洲国产| 亚洲av电影在线观看一区二区三区| 中文字幕久久专区| 亚洲激情五月婷婷啪啪| 久久国产亚洲av麻豆专区| 秋霞伦理黄片| 精品午夜福利在线看| 国产日韩欧美亚洲二区| 在线观看免费视频网站a站| 精品一区二区免费观看| 国产综合精华液| 亚洲欧美成人精品一区二区| 亚洲欧美中文字幕日韩二区| 久久久久久久亚洲中文字幕| 日日撸夜夜添| 丰满饥渴人妻一区二区三| 一级,二级,三级黄色视频| 王馨瑶露胸无遮挡在线观看| 精品酒店卫生间| 18禁在线播放成人免费| 国产又色又爽无遮挡免| 国产男女内射视频| 中文字幕久久专区| 日韩中文字幕视频在线看片| videos熟女内射| 少妇人妻 视频| 久久久久久久国产电影| 久久国产乱子免费精品| 一级av片app| 一区在线观看完整版| 欧美成人精品欧美一级黄| 日本与韩国留学比较| 欧美3d第一页| 日本av免费视频播放| 妹子高潮喷水视频| 中文欧美无线码| 丰满饥渴人妻一区二区三| 国产黄色免费在线视频| 一本一本综合久久| 国产免费福利视频在线观看| 一区二区三区乱码不卡18| 最近中文字幕2019免费版| 国产高清有码在线观看视频| 亚洲在久久综合| 中文字幕av电影在线播放| 国产探花极品一区二区| 亚洲精品久久午夜乱码| 91精品一卡2卡3卡4卡| 欧美性感艳星| 少妇被粗大的猛进出69影院 | 精品久久久噜噜| 2018国产大陆天天弄谢| 一级毛片aaaaaa免费看小| 亚洲欧洲日产国产| 91精品一卡2卡3卡4卡| 久久ye,这里只有精品| 18+在线观看网站| 男女边吃奶边做爰视频| 久久久欧美国产精品| 男女国产视频网站| 国产成人午夜福利电影在线观看| 成人特级av手机在线观看| 国产又色又爽无遮挡免| 一级二级三级毛片免费看| 亚洲性久久影院| 国精品久久久久久国模美| 亚洲中文av在线| 丰满少妇做爰视频| 一级毛片黄色毛片免费观看视频| 国产亚洲精品久久久com| 成年美女黄网站色视频大全免费 | 啦啦啦啦在线视频资源| 日本午夜av视频| 99热这里只有精品一区| 少妇人妻 视频| 亚洲av在线观看美女高潮| 99热网站在线观看| 在线观看一区二区三区激情| 18禁裸乳无遮挡动漫免费视频| 99热这里只有是精品50| 国产男人的电影天堂91| 亚洲丝袜综合中文字幕| 精品久久久久久久久av| 日韩电影二区| 国产亚洲午夜精品一区二区久久| 两个人免费观看高清视频 | 亚洲精品一二三| 免费大片18禁| 亚洲中文av在线| 欧美日韩国产mv在线观看视频| 高清午夜精品一区二区三区| 日韩亚洲欧美综合| 伊人亚洲综合成人网| 欧美精品亚洲一区二区| 欧美亚洲 丝袜 人妻 在线| 婷婷色av中文字幕| 国内揄拍国产精品人妻在线| 久久久久久久久久久丰满| 亚洲精品乱久久久久久| 中文字幕久久专区| 丰满乱子伦码专区| 久久久久久久久久成人| 成人18禁高潮啪啪吃奶动态图 | 在线观看人妻少妇| 免费不卡的大黄色大毛片视频在线观看| 夜夜爽夜夜爽视频| 国产精品三级大全| 亚洲人成网站在线观看播放| 在线观看国产h片| 在线精品无人区一区二区三| av在线播放精品| 国产av精品麻豆| 亚洲av免费高清在线观看| 久久免费观看电影| 97精品久久久久久久久久精品| 如日韩欧美国产精品一区二区三区 | 亚洲av男天堂| av在线app专区| 在线观看三级黄色| 午夜视频国产福利| 一级毛片久久久久久久久女| 看十八女毛片水多多多| 9色porny在线观看| 看十八女毛片水多多多| 午夜福利影视在线免费观看| av网站免费在线观看视频| 国产日韩欧美视频二区| 午夜av观看不卡| 久久久久久人妻| 一级毛片 在线播放| 免费人成在线观看视频色| 国产高清国产精品国产三级| 国产欧美亚洲国产| 一本久久精品| 久久99一区二区三区| 丝袜脚勾引网站| 日本av免费视频播放| 免费黄网站久久成人精品| 国产成人91sexporn| 在线观看一区二区三区激情| 国产精品三级大全| 国产在线视频一区二区| 欧美日韩国产mv在线观看视频| 人体艺术视频欧美日本| 亚洲综合精品二区| 成人毛片a级毛片在线播放| 青春草国产在线视频| 在线观看免费日韩欧美大片 | 超碰97精品在线观看| 国内揄拍国产精品人妻在线| 在线观看一区二区三区激情| 91久久精品国产一区二区三区| 国产一区二区三区综合在线观看 | 在线观看av片永久免费下载| 91精品国产九色| 日日撸夜夜添| 丝袜喷水一区| 97在线视频观看| 久久鲁丝午夜福利片| videossex国产| 国产黄片美女视频|