SAMPLE OPTIMIZATION OF ENSEMBLE FORECAST TO SIMULATE TROPICAL STORMS (MERBOK, MAWAR, AND GUCHOL) USING THE OBSERVED TRACK

Ensemble forecasting is a method in which the initial state of a deterministic forecast is considered perturbed, and a number of ensemble members are integrated ^[1]. In recent years, ensemble forecasting has gained popularity in numerical weather prediction (NWP) because its application provides more significant economic and social benefits than a single "best guess"forecast does. However, a few studies have found that ensemble builds based on initial perturbations generally have insufficient dispersion ^[7]. Furthermore, some studies suggest that prediction uncertainty is best simulated using ensembles containing both modeling uncertainties and the influence of initial conditions ^[8]. However, there is usually no reliable reference point to represent initial uncertainties because uncertainties simulated by different ensemble systems are different ^[9].

With limited members, an ensemble may not produce a perfect Gaussian probability distribution; however, each ensemble member can be regarded as a sample from an approximate Gaussian probability distribution, including those with close to zero probability when a sample is outside the distribution range ^[9]. Ideally, samples should reflect the true state of the atmosphere (observation); however, in practice, samples do not represent the true state. Especially, some samples are outliers and small probability events. At times, these samples largely deviate from verified observations. These impossible or outlier samples are meaningful for statistical analysis and may downgrade the accuracy of ensemble forecasts. One method to solve this problem is to transform the probability density function (PDF) for the truth to a PDF for observations to be verified ^[1]. In this study, observation data were used to limit the probability distribution of samples. This method is different from data assimilation. In general, equal weights are calculated using a simple arithmetic average, whereas unequal weights are determined using more complex techniques. This principle can be applied to obtain optimal ensembles, and it has been applied successfully to forecast both warm-season meso-scale convective systems^[10] and cold-season topographical forcing ^[11].

Track forecasts and maximum surface wind speed forecasts can be used to quantitatively assess risks from tropical storms (TSs), thereby allowing for earlier and more appropriate decisions regarding coastal evacuations ^[12]. Typically, TSs originate far from land; therefore, data assimilation is not very effective. During the past 20 years, track forecast errors from day 1 to 5 have reduced by more than 50%. Nevertheless, improvements in NWP are still required to decrease TS tracks. Sample optimization is also known as observation-based ensemble subsetting, and it is also the objective of previous works^[13-16]. Qi et al. ^[14] first proposed a method for obtaining ensemble mean track forecasts by choosing members weighted as a function of their short-term track forecasted error. Dong and Zhang ^[15] applied this method to a single ensemble after modifying the subset selection. The China Meteorological Administration discovered that for a given forecast ensemble, some samples produce small track forecast errors ("good" samples), but some produce poor forecast tracks ("bad" samples) that deviate substantially from the truth. Therefore, good members should be objectively identified by some standards. In this study, samples close to expected values (observation) were identified as"good", and those with small probabilities were regarded as"bad". The bad samples were replaced by good ones to reflect the true state as far as possible. To determine the effect of sample optimization on ensemble forecasting, some studies, such as Dong and Zhang^[15] used the best-track position of hurricanes that occurred during 2012-2013. However, this best track data is not available in time during runtime. Therefore, TS track data from the Central Meteorological Station (http://typhoon.nmc.cn/web.html) were used in this study as the criterion for restricting the probability distribution.

Houtekamer ^[17] showed that a spread has the most predictive value when it is"extreme", being very large or small compared to the corresponding mean value. Whitaker ^[18] pointed out that the more a spread departs from the corresponding climatological mean value, the more useful it is as a predictor of skill. Grimit and Mass ^[19] used a statistical model based on an ensemble prediction system to obtain a perfect forecast assumption. Generally, a larger (smaller) ensemble dispersion implies more (less) uncertainty in forecast ^[20]. However, an excessive number of bad samples being replaced by good ones may lead to an extremely small ensemble spread, which will produce inaccurate ensemble mean values. In contrast, insufficient replacement of bad samples may prohibit the adjustment and upgrade of the ensemble mean. Therefore, the number of samples to be replaced should be carefully determined so that the ensemble mean is upgraded and the ensemble spread is not excessively downgraded.

Relevant work has been done by Li ^[21] in this field. However, the subject of study was typhoon (a stronger tropical cyclone) in that work. To balance the relationship between ensemble spread and track errors, Li replaced approximately ten to twenty-five samples and it showed good performance ^[21]. Therefore, in this study, twenty bad samples were replaced by good ones for convenience, and the subject of the research was a tropical storm (a weaker tropical cyclone). The effectiveness of sample optimization was tested.

This paper is organized as follows. The scheme, datasets, and analysis of the cycles are presented in Section 2. The analysis of the results is described in Section 3, and concluding remarks are given in Section 4.

According to Li ^[21], two main steps were followed for sample optimization for ensemble forecasts, i. e., sample selection and supplementation.

In the first step, samples were selected, and their qualities were evaluated. Because track is a very important indicator of TSs, the observed track data obtained from the China Meteorological Administration was selected as the standard and the samples were compared to this standard. Track was used to restrict the probability density distribution of the samples, with the principle of absolute track forecast error as the basis; absolute track forecast error is defined as the great-circle distance between the forecast position and observed-track position of a TS at a given time. According to this definition, a sample is considered good when it is close to the observed position. The inherent premise of this technique is that samples with smaller track error result in better performance for longer forecast lead times on average. It is worth mentioning that only track was chosen as the standard among many other criteria for sample selection, such as radar echo and intensity. Samples selected based on other standards will be investigated in future studies. In this study, we focused only on whether good and bad samples can influence ensemble forecasts; we used the following procedures to do so. First, the track error between each sample and observation was calculated and classified. Through comparative analysis, samples with a lower track error were identified as good samples, and samples with a higher track error were identified as bad samples. Good samples were copied, and the bad ones were removed. Neither type of sample is called a moderate sample.

In the second step, after elimination of the bad samples, good samples were supplemented with their copies to maintain the total number for the ensemble forecast. Further, perturbations were added to the supplemented samples using EnKF (EnSRF), that is, the copied samples acquired assimilation increment.

Figure 1 is a schematic diagram of the sample optimization process for ensemble prediction in one cycle. A 1-h lag was observed between the time when the sample optimization spins up and the time at which TS track forecast data were first made available at the observed position. Although the ensemble comprised a total of N samples in this study (N=60), for convenience the schematic diagram shows the results of only 8 samples (N=8), and the forecast track is a simple arithmetic average of the positions of selected samples or M samples. Conventional ensemble forecasts are called SOno, and ensemble forecasts with sample optimization are called SO. Ensemble members 4 and 5 (i. e., M= 2 and N= 8) are good samples, based on the observed position 1 h after the spin-up time. The exact number of selected good samples is based on past performance and this number is constant for each experiment (introduced in Section 2 b).The method used in this study is different from that used by Dong and Zhang ^[15], in which ensemble subsetting only contains selected good samples and both M and N vary. The method used in this study has a fixed N and only M is variable. M bad samples are replaced by M good samples, and there are N-2M medium samples.

Figure 1.  Schematic illustration of the 1-h difference between the time when ensemble members arrive and the time when the ensemble mean arrives.

In this study, three TSs were considered: Merbok, Mawar, and Guchol. These were the second, sixteenth, seventeenth TSs in 2017, respectively, and all did not reach typhoon level.

Data for Merbok were compiled by the Central Meteorological Observatory at 0000 UTC June 10, 2017. Merbok then intensified into a TS at 0600 UTC June 11, 2017 and landed on the coast of the Dapeng Peninsula in Shenzhen at around 1500 UTC June 12, 2017. At that time, it had a maximum wind speed of 25 m s^-1.

Data for Mawar was compiled on September 1, 2017. Around 21: 30 on September 3, 2017, Mawar landed in Lufeng City, Guangdong Province.

From September 6 to 7, Guchol crossed the Bashi Channel and entered the Taiwan Strait. There, it combined with the cold air front and turned into an extratropical cyclone.

To examine the impact of sample optimization on ensemble forecasts, seven experiments were designed for each TS. All experiments were performed using the Advanced Research Weather and Research Forecasting (ARW-WRF) model, version 3.4.1. In this study, the ensemble was initialized using ECM / ECMWF 0.125° × 0.125° analysis data for each TS.

Figure 2 shows the tracks of the three TCs. The horizontal grid spacing is 3 km. Merbok, Mawar, and Guchol have 298 × 667, 304 × 262, and 409 × 445 grid points, respectively, all with 35 vertical levels. The Grell Devenyi cumulus scheme^[22], WRF single-moment (WSM) 6-class microphysics scheme with graupel^[23], and Yonsei State University (YSU) scheme^[24] for planetary boundary layer (PBL) processes were used in this research.

Figure 2.  The tracks of three TCs. (a) Merbok, (b) Mawar, and (c) Guchol.

In particular, the sample optimization experiments comprised of three steps:

(1) Spin-up (6h forecast without data assimilation). The initial value (created using the WRF-3DVAR method) and the major physical parameters were the same as those in Zhu et al. ^[25]. The WRF-based EnKF system was first developed for regional-scale data assimilation by Meng and Zhang ^[26-27]. The control variables are stream function, velocity potential, and unbalanced pressure. The initial ensemble members were generated by adding perturbations that were randomly sampled from the default"cv3"background error covariance option in the WRF 3DVAR package ^[28]. Similar perturbations were used to represent the boundary condition uncertainties of the ensemble. The covariance relaxation method was the same as that of Zhu et al. ^[25].

The perturbed variables included horizontal wind components (u, v), potential temperature, mixing ratio for water vapor, and a standard deviation of 2 m s^-1 for wind, 1 K for temperature, and 0.5 g kg^-1 for mixing ratio ^[26-27]. Similar perturbations were used to represent the boundary condition uncertainties of the ensemble. The covariance relaxation method, with a relaxation coefficient of 0.8, proposed by Zhang ^[29] was used to infl ate the background error covariance. The prognostic variables of perturbation potential temperature (T), vertical velocity (W), horizontal wind components (u, and v), mixing ratio for water vapor (QVAPOR), cloud water (QCLOUD), rainwater (QRAIN), perturbation geopotential (PH), perturbation dry air mass in column (MU), surface pressure (PSFC), and perturbation pressure (P) were updated according to Zhu et al. ^[25]. The horizontal length of the covariance localization was set to 30 km and the vertical length was set to 6 layers.

(2) Data assimilation was performed using the WRF-EnKF method for 6h cycles after the spin-up. The configuration of this process was the same as that of Zhu et al. ^[25], with the only difference being the inclusion of three steps in each cycle: selection, assimilation, and integration. In other words, the replacement of good and bad samples (proposed in Section 2 a: Sample optimization for ensemble forecasts) was conducted before the application of EnKF, similar to Zhu et al. ^[25].

In EnKF, α is a constant. HP^bH^T and R are scalars representing the background and observational-error variance at the observation location ^[30].

In this work, some differences remained between the supplemental samples and good ones. When assimilated, α is not equal to 0 for good samples but is set to 0 for the supplemental ones.

According to Table 1, SOno was conventional ensemble forecast experiment, and SO20 (track) was conducted for sample optimization. In other words, SO20 (track) is the case in which twenty good samples were replicated, twenty bad samples were deleted, and twenty moderate samples were retained. Furthermore, the sample number was maintained at 60 throughout the sample optimization process.

(3) The initial value (ensemble average) was generated via the final analysis in step 2, and then forecast to the end. Fig. 3 is a flow chart of the process.

Figure 3.  Flowchart for sample optimization experiments.

Table 2 shows the start and end time of ensemble forecast and sample optimization. For example, 2017061012 means 1200 UTC on June 10, 2017. Merbok began at 1200 UTC on June 10, 2017 and ended at 0000 UTC on June 13, 2017. Mawar started at 0000 UTC on September 1, 2017 and ended at 0000 UTC on September 4, 2017. Guchol started at 1200 UTC on September 4, 2017 and ended at 1200 UTC on September 7, 2017. Their sample optimization began at 1800 UTC on June 10, 2017, 0600 UTC on September 1, 2017, and 1800 UTC September 4, 2017, and ended at 0000 UTC on June 11, 2017, 1200 UTC on September 1, 2017, and 0000 UTC on September 5, 2017, respectively, in all experiments.

The influence of sample optimization on ensemble spread was also investigated. Fig. 4 shows the evolution of ensemble spread and the 6-h evolution of mean track errors in three experiments.

Figure 4.  (a) 6-h evolution of ensemble spread and (b) 6-h evolution of mean track errors. Merbok (1900 UTC on June 10, 2017 to 0000 UTC on June 11, 2017), Mawar (0700 UTC to 1200 UTC on September 1, 2017) and Guchol (1900 UTC on September 4, 2017 to 0000 UTC September 5, 2017).

According to the figure, sample optimization can certainly make ensemble spread reduce. Next, it can be easily seen that the ensemble spread reduced very quickly during the first 3-h cycles but decreased relatively mildly in the next 3-h cycles. Finally, the ensemble spread was very close in the last cycle for all three experiments.

It can be seen from Fig. 4b that sample optimization can reduce the track error in each cycle. The variation trend of track errors seen in Fig. 4b is very similar to that of ensemble spread seen in Fig. 4a. With further sample optimization, the magnitudes of the track errors in the last three cycles reduce to values smaller than those in the first three cycles.

Figure 5 displays the relationship between ensemble spread and track error for each cycle. It can be clearly seen in the figure that the ensemble spread and track error reduce at the same time. With passage of time, the declines in the ensemble spread and track error reduce. Both of them have a value of approximately 20 km in the last two cycles. There is little need for further sample optimization, and 6-h cycles may be enough.

Figure 5.  The relationship between ensemble spread (km) and track error (km) in 6-h cycles. (a) to (f) correspond to the first to the sixty cycle. For the icons of the same color, the lower left one represents those before optimization, and the upper right one represents those after optimization in the same figure.

Figure 6 shows the change in temperature anomaly with time and height for 6-h cycles. The intensities of the warm cores were effectively weakened by sample optimization. A comparison of Fig. 6a and 6b shows that there is an obvious positive temperature anomaly under 3 km both in the case of SOno and that of SO20(track); however, there is still a stronger positive temperature anomaly at the height of 12 km for SOno. Further, the stronger positive temperature anomaly at the height of 12 km disappeared after sample optimization. Analogously, in Fig. 6c, large scale positive temperature anomalies are seen in the middle and upper layers for SOno. However, the positive temperature anomaly above 1 km completely disappeared in the case of SO20 (track). In this respect, the change for Guchol is similar to those for Merbok and Mawar, and detailed discussion of this point is not required here.

Figure 6.  The change in temperature anomaly (℃) with time and height (m). (a) (SOno) and (b) (SO20 (track)) represent Merbok (1800 UTC on June 10, 2017 to 0000 UTC on June 11, 2017). (c) (SOno) and (d) (SO20 (track)) represent Mawar (0600 UTC to 1200 UTC on September 1, 2017). (e) (SOno) and (f) (SO20 (track)) represent Guchol (1800 UTC on September 4, 2107 to 0000 UTC September 5, 2017). x-axis is time (hour), and y-axis is height (m).

Figure 7 shows the change in temperature anomaly (℃) with height (hPa) by reanalysis data. It should be noted that there are only reanalysis data of 0000 and 1200. In other words, for a particular TS, the time seen in Fig. 7 corresponds to the last time seen in Fig. 6.

Figure 7.  The change in temperature anomaly (℃) with height (hPa) determined by analysis. (a) Merbok (0000 UTC June 11, 2017). (b) Mawar (1200 UTC September 1, 2017). (c) Guchol (0000 UTC September 5, 2017).

The relationships between air pressure and altitude are roughly as follows: 1000 hPa for approximately 0 km, 700 hPa for approximately 3 km, 500 hPa for approximately 5.5 km, 300 hPa for approximately 9 km, 200 hPa for approximately 12 km, and 100 hPa for approximately 16 km.

It can be seen from a comparison of Figs. 6a, 6b, and 7a that the simulated temperature anomaly was larger than the observed temperature anomaly for the reanalysis data below 2 km; however, the temperature anomaly simulated in SOno was much larger. Between 6 km to 9 km, the temperature anomaly simulated in SO20 (track) was larger than that of analysis data. Further, the temperature anomaly of the analysis data was larger than that simulated in SOno. For 12 km to 15 km, the temperature anomaly of the analysis data was approximately in the range 2~4℃; that in SO20 (track) was close to this range, but that in SOno was smaller than this range.

On comparison of Fig. 6a and 6d with Fig. 7b, it can be seen that the temperature anomaly simulated in SOno above 1 km slightly changed and then remained in the 1.5~2℃ range. Moreover, the temperature anomaly simulated in SOno was about 5~6℃ below 1 km. These values are largely deviated from the analysis data. The temperature anomaly simulated in SO20 (track) was similar to the analysis data below 1 km and that around 9 km, but larger than the analysis data around 3 km and 12 km.

From Figs. 6e, 6f, and 7c, it is obvious that the temperature anomaly simulated in SO20 (track) was much larger than that simulated in SOno, and close to the analysis data, above 9 km.

Therefore, by sample optimization, the vertical distribution of TS can be brought closer to reality.

In summary, sample optimization is helpful for adjusting the track, intensity, and structure of TSs.

A comparative analysis of the prediction results was performed. Fig. 8 shows the observed and simulated TS tracks. It can be seen that the tracks simulated in SO20 (track) are closer to the observation than those simulated in SOno, especially at the starting time (in these figures).

Figure 8.  Tracks of (a) Merbok (0000 UTC Jun 11, 2017 to 0000 UTC Jun 13, 2017), (b) Mawar (1200 UTC Sep 1, 2017 to 0000 UTC Sep 4, 2017), and (c) Guchol (0000 UTC Sep 5, 2017 to 0900 UTC Sep 7, 2017) as per WRF data.

Table 3 shows the average track and intensity errors between the observations and simulations. All tracks errors in the simulations in SO20 (track) were smaller than those for SOno. The track error for Guchol decreased by approximately 42% after sample optimization.

The SO20 (track) track errors for typhoons Rammasun, Nida, and Megi were observed to have decreased by approximately 12.7%, 4.3%, 16.4% respectively in the study by Li (2018). Furthermore, the track errors for TSs Merbok, Mawar, and Guchol reduced by approximately 29.2%, 7.1%, 42.1%, respectively.

To some extent, the track errors for weaker tropical cyclones reduced more significantly than those for stronger tropical cyclones in the case of SO20 (track).

Also, in the majority of cases, the intensity error reduced in SO20 (track).

Figures 9 and 10 show the development of TS intensity (minimum sea level pressure (SLP) and maximum wind speed). In the two figures, it can be seen that the intensity adjustment is not perfect. At first, the minimum SLP by simulated was lower than that by observation in for the three TSs. In the end, the minimum SLP was higher than the SO20 (track) observation for these cases. Furthermore, the evolution of the minimum SLP by SO20 (track) was similar to that of the observation. In Fig. 10a and 10b, it can be seen that the simulated maximum wind speed was much larger than the observation, but through sample optimization, the simulated maximum wind speed was improved successfully. The exception to this can be seen in Fig. 10c where the simulated maximum wind speed in SO20 (track) was far from the observation. Therefore, SO20 (track) can improve the accuracy of track prediction but cannot definitely improve intensity prediction. Hence, there is room for further development of this sample optimization scheme.

Figure 9.  The minimum SLP (hPa) errors for (a) Merbok (0000 UTC Jun 11, 2017 to 0000 UTC Jun 13, 2017), (b) Mawar (1200 UTC September 1, 2017 to 0000 UTC September 4, 2017), and (c) Guchol (0000 UTC September 5, 2017 to 0900 UTC September 7, 2017) as per WRF data.

Figure 10.  The maximum wind speed (m s^-1) errors for (a) Merbok (0000 UTC Jun 11, 2017 to 0000 UTC Jun 13, 2017), (b) Mawar (1200 UTC September 1, 2017 to 0000 UTC September 4, 2017), and (c) Guchol (0000 UTC September 5, 2017 to 0900 UTC September 7, 2017) according to WRF data.

In this study, we performed a series of simulations to achieve samples that are as close to the true state as possible. We took the observed track data from the China Meteorological Administration as standard during the process of restricting the probability distribution of the samples. We then investigated the influence of sample optimization on ensemble prediction. Finally, we conducted a comparative analysis of the results to determine the effect of sample optimization on the vertical structure of TSs.

In previous works, typhoons (stronger tropical cyclones) were the objects of research, and good results had been obtained. Therefore, in order to prove the applicability of sample optimization, we considered three TSs in our simulations.

In 6-h cycles, the track, intensity, and vertical structure were all found to have adjusted better and were closer to the observation. It was found that the vertical structures were more reasonable and closer to the reanalysis data in the case of SO20(track) than in the case of SOno. Therefore, it may be concluded that the various thermal structures that result from sample optimization lead to different TS evolutions.

The prediction results show that sample optimization outperforms conventional ensemble forecasts (SOno) with respect to TS track simulations. Thus, sample optimization was found to be beneficial for updating tracks of TSs.

Sample optimization can be used to simulate both stronger tropical cyclones (such as typhoons) and weaker ones (just like tropical storms).

However, despite the better results shown in this study, there are potential drawbacks to this new optimization method. For instance, the criterion of sample optimization is based on the distance between the forecast position of a TS and the corresponding observed track position. Also, the intensity of TSs was not improved significantly. Thus, intensity or other indicators should be included in the future work.

In short, double standard (such as track plus intensity) may be applied in future.

Acknowledgments: The authors would like to thank all the members of the Guangzhou Institute of Tropical and Marine Meteorology for their constructive remarks.