Comparisons of Initial Condition Perturbation Methods for Regional Ensemble Forecasts of Wind Speed in Gansu of China

1. Introduction

According to the Global Wind Energy Council (GWEC) [1], the total global capacity of wind power reached a milestone of 837 GW. Wind energy will likely continue to grow strongly and play a leading role in achieving a low-carbon or net-zero future. However, the large scale penetration of wind power also brings many challenges due to fluctuations and intermittency of wind power generation [2,3]. The wind forecast plays a critical role in overcoming these challenges [4,5,6]. One of the essential benefits of accurate wind forecast is to reduce grid stress and reserve requirements [7,8].

The wind forecast for energy generation and power system operations mainly focuses on the immediate short term of seconds to minutes, the short term of hours up to two days, and the medium term of 2 to 7 days as the power system operations are carried out within these time windows [9]. Nowadays, there are mainly three classes of wind forecasting methods. The first is the statistical method based on historical data, and the second is physics-based numerical weather prediction (NWP) models. The last one is the hybrid approach combining different approaches, such as combinations of statistical and physics-based approaches. The importance of different methods varies with the forecast lead time [10]. Hanifi et al. [11] systematically reviewed these approaches and concluded that NWP models significantly benefits for forecasts beyond 6 hours. However, NWP simulations are uncertain due to uncertainties in the initial conditions, limited understanding of the atmosphere's physical processes, and the atmospheric flow's chaotic nature [12,13]. The ensemble prediction system (EPS) is a promising way to estimate the forecast uncertainties [14,15].

The EPS is often made by perturbing initial conditions, and the perturbation methods applied by different global operational ensemble forecasting centers are different. At the National Centers for Environmental Prediction [NCEP, previously the National Meteorological Center (NMC)), Toth and Kalnay [16] introduced the breeding of growth modes (BGM) method based on the argument that fast-growing perturbations develop naturally in a data assimilation cycle [17]. The European Centre for Medium-range Weather Forecast (ECMWF) developed and implemented the Singular Vector (SV) method to identify the directions of the fastest growth [18,19]. The Canadian Meteorological Centre (CMC) developed an ensemble data assimilation method to produce different initial conditions for ensemble forecasts [20,21]. As a result of limited computing resources, the global EPS is usually run at coarser resolutions than deterministic forecasts. The NCEP global ensemble forecast system (GEFS) uses a 34 km horizontal resolution [22], while the CMC Global Ensemble Prediction System (GEPS) runs at a horizontal resolution of 39 km [23]. The ECMWF EPS uses the highest horizontal resolution at 18 km and 91 vertical levels, containing one control member and 50 perturbed members [24].

The region of interest for this study, Gansu province in China, is rich in wind resources [25] and has the world’s most giant onshore wind farms [26]. Unfortunately, Gansu is also the province with the most severe wind curtailment in China, discarding 10.4 TW h of potential wind power in 2016 [27]. According to Lew et al. [28], a 10% improve ment in wind forecasts could lead to a 4% reduction in curtailment and operation costs. Therefore, improving wind forecasts for wind farms in Gansu is very important. In addition, Gansu is located in northwest China and has a complex topography, requiring higher spatial resolutions to resolve topographic impact [8,29].

Given that the horizontal resolution of the global EPS is too coarse, it is necessary to build a regional EPS for an accurate wind forecast. The construction of the initial condition perturbations and lateral boundary condition (LBC) perturbations are crucial for a skillful regional EPS. The most common approach is the dynamical downscaling of a global EPS to the regional domain [30]. Because of its simplicity and low computational costs, this method is implemented by many NWP centers for regional operational EPS [31,32,33]. However, the dynamical downscaling method fails to represent the small-scale uncertainties resolved by the regional model [34]. Thus, researchers use regional versions of traditional perturbation methods such as BGM, SV, and ensemble transform Kalman filter (ETKF) [35,36] and that produce more information on small-scale uncertainties. Also, Caron [37] found that the mismatches between the initial condition perturbations and the LBC perturbations cause spurious perturbations. Therefore, a blending method was proposed to combine the regional model-based small-scale initial condition perturbations with large-scale perturbations from a global EPS [38,39]. Wang et al. [39] described the blending method implemented in the regional EPS, i.e. Aire Limitée Adaption dynamique Développement International-Limited Area Ensemble Forecasting (ALADIN-LAEF), and demonstrated that the blending method outperforms the dynamical downscaling and breeding method. Zhang et al. [40] also showed that the breeding method improved the ensemble spread and forecast skills of the Global/Regional Assimilation and Prediction Enhanced System (GRAPES) Regional EPS (GRAPES-REPS).

In this study, we use the Weather Research and Forecasting (WRF) model for dynamical downscaling of ECMWF EPS to generate large-scale perturbations and the BGM method to generate small-scale initial condition perturbations due to its clear meaning and low computational cost. As the BGM method requires using forecast error for calculating scaling factor and observations are not always available, we proposed an alternative to calculating the scaling factor. Additionally, we apply the blending method to combine perturbations of different scales and compare the wind forecast performance of dynamical downscaling, BGM, and blending in Gansu.

The paper is structured as follows: Section 2 describes the WRF model setup and regional EPS using dynamical downscaling, BGM, and blending methods. Section 3 introduces data and metrics for evaluation. Section 4 presents the evaluation results of day-ahead and ultra-short wind forecasts for one month. Finally, section 5 concludes the study with suggestions for future work from the results.

2. WRF Model Setup and Perturbation Methods Description

2.1. WRF Model Setup

The WRF model version 3.9.1 is used. As shown in Figure 1, the WRF model is configured with one single domain at a horizontal grid resolution of 8 km centered at 38° N in latitude and 101° E in longitude. The model used a terrain-following vertical coordinate with 55 vertical levels and a model top at 50 hPa. The ECMWF EPS is run four times per day, and only the ensemble forecast data initialized at 12 UTC is used to run the WRF model for 54 hours of forecasts. The physics parameterizations are selected following the same methodology presented in the [41]. The selected schemes include WRF Single-Moment 6-Class for microphysics [42], Bougeault-Lacarrère (BouLac) planetary boundary layer (PBL) scheme [43], Pleim-Xiu land surface model [44,45], Kain-Fritsch scheme for cumulus [46], and New Goddard for both longwave and shortwave radiation [47].

2.2. Description of the Initial Condition Perturbation Methods

As listed in Table 1, this study compared three initial condition perturbation methods: dynamical downscaling, BGM, and blending. The WRF model is used for dynamically downscaling all the 51 ensemble members of ECMWF EPS for initial condition and LBC perturbations.

2.2.1. BGM Method

The BGM method is an initial condition perturbation method widely applied for ensemble prediction. In this work , the classic BGM method is implemented following [48]. Firstly, an initial random perturbation determined by Equation (1) is added at a given time $t$:

$$P\left(z\right)=\omega RE\left(z\right)$$

(1)

where $z$ is the state variable of the NWP model ($z$ could be $T,$ $U$, or $V$ representing the thermal and dynamic fields, respectively), and $\omega$ is an adjustment coefficient to control the magnitude of the initial random perturbations. $R$ is a uniform random number in the interval [−1, 1]. $E\left(z\right)$ is the root mean square error (RMSE) for variables at each layer.

Then, the model with both the unperturbed and perturbed initial condition is run for a short period (i.e., 12 hours for all the experiments in this paper). The difference between the control prediction with an unperturbed initial condition and the perturbation prediction with an unperturbed initial condition is calculated. As the difference increases with time, the difference is scaled down to have the same norm (e.g., amplitude) as the initial perturbation. This scaled difference is added to the analysis of the following time $t+1$. The process is repeated forward in time to generate the final perturbations until the growth rate of perturbations is saturated. The detailed calculations of perturbations are as follows:

$$p_{t+1}=c_t\left(f_t^a-f_t^s\right)$$

(2)

$$c_t=E\left(p_t\right)/E\left(p_{t+1}^{\text{'}}\right)$$

(3)

where $p_{t+1}$ is the perturbation at time $t+1$ of the next cycle; $f_{t+1}^a$ is the perturbed prediction, and $f_{t+1}^s$ is the control prediction. $c_t$ is the scaling factor; $E\left(p_t\right)$ and $E\left(p_{t+1}^{\text{'}}\right)$ are RMSE at the beginning and end of the cycle period.

T Traditionally, observation data is needed to calculate the scaling factor. However, as observation data is not always available, we proposed to use the min-max scaling technique, a widely used normalization technique in machine learning, to scale the perturbations between $-a$ and $a$ without observations:

$$p_{t+1}=\frac{p_{t+1}^{\text{'}}-\min \left(p_{t+1}^{\text{'}}\right)}{\max \left(p_{t+1}^{\text{'}}\right)-\min \left(p_{t+1}^{\text{'}}\right)}2a-a$$

(4)

where $a=\omega E\left(z\right)$ so that perturbations at the end of the cycle are of the same amplitude as initial perturbations.

2.2.1. Blending

The blending method combines the small-scale uncertainties resolved by WRF BGM and large-scale features from the global ensemble. Also, including perturbations from the global ensemble in the initial condition perturbation ensure that the initial condition is consistent with LBC provided by the global ensemble.

Here, the blending is implemented following the steps: (1) downscaling the global ensemble initial condition to the regional domain using the WRF model to get large-scale uncertainties; (2) sum the downscaled initial condition with perturbations generated by the WRF BGM described in Section 2.2.1.

3. Data and Metrics for Evaluation

3.1. Observation Data

The study period is from October 1st, 2020, to October 31st, 2020, for evaluation. We use hourly wind speed observations by anemometers installed at wind turbine hub height from 28 wind farms in Gansu, China. The WRF model output is interpolated to the turbine hub height to compare with the observation data. In addition, the wind speed forecast of the nearest grid point from the WRF domain is extracted to compare with wind turbine observations.

3.2. Evaluation Methods

According to the National Energy Bureau (NEB) regulation, day-ahead short-term and four-hours ahead wind forecasts are required [49]. Using the 12 UTC ECMWF data, the 28 to 51 hours forecast horizon represents the day-ahead forecast, while the 10 to 13 hours forecast horizon corresponds to the ultra-short term forecast.

To compare the performance of ensemble forecasts generated by the three initial condition perturbation methods, we applied several verification methods: RMSE and mean bias error (MBE) for the ensemble mean, ensemble spread [i.e., standard deviation (std)], continuous ranked probability score (CRPS) [50,51], and rank histograms [52]. These measures are calculated at the 1-hour frequency for wind speed forecasts following the equations below:

$$RMSE=\sqrt{\frac{1}{N}{{\displaystyle \sum }}^{}{\left(W{S}_{pred}-W{S}_{obs}\right)}^2}$$

(5)

$$MBE=\frac{1}{N}{{\displaystyle \sum }}^{}\left(W{S}_{pred}-W{S}_{obs}\right)$$

(6)

$$std=\sqrt{\frac{1}{N}{{\displaystyle \sum }}^{}{\left(x_{pred}-\overline{x_{pred}}\right)}^2}$$

(7)

where $W{S}_{pred}$ is the predicted wind speed ($WS$), $W{S}_{obs}$ is the observed wind speed, and $N$ is the number of pairs of forecast and observation. $x_{pred}$ is the ensemble member forecast, and $\overline{x_{pred}}$ refers to the ensemble mean. RMSE and MBE of the ensemble mean describe the deterministic skill of the ensemble forecast. The ensemble spread is the std of the ensemble members with respect to the ensemble mean. When observations are unavailable, the ensemble spread can be used as a predictor of skills of the ensemble mean [53]. Small spread indicates low uncertainties, otherwise indicates high uncertainties.

The CRPS is defined as:

$$CRPS={\displaystyle {{\displaystyle \int }}_{-\infty }^{\infty }}{\left[F\left(y\right)-F_o\left(y\right)\right]}^{2}dy$$

(8)

where $F_o\left(y\right)$ is an indicator function of which the value is 0 if the forecast variable $y$ is less than the observation and one otherwise [54].

$$F_o\left(y\right)=\{\begin{array}{c}0, \mathrm{y}<\mathrm{o}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d} \mathrm{v}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e}\\ 1, \mathrm{y}\ge \mathrm{o}\mathrm{b}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{d} \mathrm{v}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e}\end{array}$$

(9)

The CRPS is a penalty score, and a smaller value of CRPS indicates a better ensemble forecast. The properscoring package in Python is used for calculating the CRPS in this work [55].

Additionally, the rank histogram is a valuable tool for examining ensemble reliability. It is generated by tallying the rank of verification (usually observation) relative to the forecast values of each ensemble member sorted in ascending order [52]. For a perfect EPS, the rank histogram should be flat. A U-shaped rank histogram generally indicates a lack of variability in the ensemble forecast, and an asymmetric shape (either J or L-shaped) indicates bias.

4. Results and Discussion

Figure 2 illustrates the one-month averaged CRPS and ensemble spread of wind speed forecast at turbine height for downscaling (solid purple lines), breeding (dashed green lines), and blending (dashed red lines) of ECMWF-EPS as a function of the forecast horizon from 10 to 54 hours. Overall, the blending ensemble shows the best performance with the smallest values of CRPS and the largest spread, especially within the forecast lead time between 10 and 25 hours. The BGM ensemble performs slightly better than the downscaling ensemble, with smaller values of RMSE and a larger spread, in the early forecast lead time. As all three ensembles use the same LBCS from ECMWF-EPS forecasts, it can be concluded that the BGM method is superior to the downscaling method, and the blending method takes advantage of both perturbation methods. After 25 hours, the difference in CRPS and ensemble spread among downscaling, BGM, and blending are almost negligible. These results suggest that the physics and boundary conditions dominate over initial condition perturbations in the long-term forecast in the study domain. As demonstrated in Figure 3, the RMSE and MBE of the ensemble mean of the downscaling, BGM, and blending ensemble for wind speed are very similar. However, during the earlier forecast lead time, the blending ensemble shows slightly smaller values than the downscaling and BGM ensemble.

Figure 4 compares the rank histogram for the wind speed forecast of the forecast lead time from 10 to 54 hours among the downscaling (blue), BGM (green), and blending (red) ensemble. The U-shaped rank histogram of all three ensembles means that the three ensemble forecasts are under dispersive. However, the blending ensemble is flatter than the downscaling and BGM ensemble, indicating that the frequency that observations lay inside the ensemble is highest for the blending ensemble. The BGM ensemble is also comparatively flatter than the downscaling.

Table 2 summarizes the averaged RMSE and MBE of the ensemble mean of wind speed forecasts for the downscaling, BGM, and blending ensemble over the forecast horizons of 10 to 13 hours and 28 to 51 hours over one month. As seen from Table 2, the RMSE of the blending ensemble is smaller than the other two ensembles, with more difference over the earlier forecast lead time from 10 to 13 hours than the forecast lead time of 28 to 51 hours, consistent with Figure 3.

The analysis above demonstrates the overall improvement of the blending ensemble over the downscaling and BGM ensemble. Also, the effect of BGM and blending is evident mainly within an earlier forecast lead time of up to 25 hours. However, the effect could be extended to a longer forecast lead time if the domain of the study increases, as it will take longer for LBC to dominate over perturbations in initial conditions (this is beyond the scope of this study, thus not shown). Therefore, although the improvement in RMSE is not significant, it is still worth applying the blending method to improve the ensemble reliability of regional EPS [40].

5. Conclusions

In this study, we compared three methods to generate initial condition perturbations for building ensemble wind speed forecasts for the Gansu province of China, including dynamical downscaling, BGM and blending. The dynamical downscaling can produce large-scale perturbations by directly downscaling the ECMWF EPS using the WRF model. As we lack the observation data that the BGM method needs for generating the scaling factor periodically, we propose to use the min-max scaling method to scale the small-scale perturbations. The WRF BGM method can resolve small-scale perturbations, while the blending method can combine large-scale and small perturbations. Over the one-month testing period from October 1st to October 31st, 2020, the blending ensemble shows the best forecast skill, especially at the early forecast lead time. The CRPS, ensemble spread, and rank histogram of wind speed forecast demonstrate that the blending ensemble has a better probabilistic forecast. Furthermore, the difference among the three perturbation methods decreases as the forecast lead time becomes negligible after 25 hours. The RMSE and MBE comparisons of wind speed forecast also indicate that the blending ensemble is superior to downscaling and BGM ensemble.