A Review on Development, Challenges, and Future Perspectives of Ensemble Forecast

1.   Introduction

Numerical weather prediction (NWP) is one of the greatest scientific advancements of the 20th century, often described as a “quiet revolution” (Bauer et al., 2015; Shen et al., 2020). Among its key achievements, four-dimensional variational data assimilation (4DVar) and ensemble forecasting are regarded as two major milestones. The atmospheric system is inherently nonlinear, characterized by chaos that is sensitive to minor initial perturbations. Even with ideal numerical models, uncertainties in the initial states—arising from errors in observations and assimilation systems—lead to significant forecast errors as the forecast lead time increases (Lorenz, 1963; Chou, 1986; Mu et al., 2003, 2011). Additionally, numerical models employ discretization schemes, sub-grid parameterizations, and other approximations to represent the dynamical and physical processes of the atmosphere. Consequently, single deterministic NWP has inherent forecast uncertainty (Li et al., 2000; Arakawa, 2004; Chen et al., 2004). Quantifying this uncertainty has thus become a critical focus for NWP. The ensemble prediction system (EPS), which may transform single deterministic numerical forecasts into probability density function (PDF) distributions and estimate the uncertainty of deterministic numerical forecasts, was first developed through the Monte Carlo method (Epstein, 1969; Leith, 1974).

Theoretically, the evolution of the PDF over time can be described by the Liouville equation. However, solving this equation is highly challenging and practically infeasible, even for a nonlinear system with only a few degrees of freedom. To address this issue, the ensemble forecasting method was developed (Molteni et al., 1996). This method utilizes specific mathematical techniques to generate a set of initial conditions, each potentially representing the true state of the atmosphere. These initial conditions are then integrated by a numerical model to produce an ensemble of forecasts, from which the PDF distribution of future weather variables can be derived. When the initial condition ensemble accurately captures the distribution of initial analysis errors and the model has adequate precision, the ensemble forecast can reliably approximate the atmospheric state’s PDF.

The core challenge in ensemble forecasting lies in appropriately representing error sources and accurately estimating uncertainties in NWP. The well-documented sources of NWP errors originate from both initial conditions and the model’s inherent limitations. Correspondingly, ensemble forecasting methods are derived from the first (initial condition) and second (model) predictability theories. Various perturbation methods are designed to represent initial condition errors and model uncertainties, leading to the development of initial condition and model perturbation methods (Chen et al., 2002; Du et al., 2018). In terms of initial condition perturbation, methods such as the singular vector (SV; Buizza, 1997), the breeding of growing mode (BGM; Toth and Kalnay, 1993), the perturbed observations (Houtekamer and Derome, 1995), the ensemble transform Kalman filter (ETKF; Wang and Bishop, 2003; Wei et al., 2006), the ensemble transform with rescaling (ETR; Wei et al., 2008), the ensemble square root filter (EnSRF; Whitaker and Hamill, 2002; Zhou et al., 2022), and hybrid perturbations combining SV perturbations with ensemble data assimilation (EDA; Isaksen et al., 2010) have been developed. In terms of model perturbation, proposed methods include the multi-physics scheme (Chen et al., 2003), the stochastically perturbed parameterization tendencies (SPPT; Buizza et al., 1999; Tan et al., 2013), the stochastic kinetic energy backscatter scheme (SKEB; Shutts, 2005; Peng et al., 2019), and the stochastically perturbed parameterization (SPP; Jankov et al., 2019; Xu et al., 2019). These studies have significantly contributed to the continuous improvement of operational EPSs in terms of overall performance and reliability.

Since the 1990s, significant progress has been made in operational ensemble forecasting, driven by advancements in ensemble forecasting theory and technology. In 1992, the ECMWF and the NCEP established the first global medium-range EPS, marking a major milestone in operational ensemble forecasting (Toth and Kalnay, 1993; Molteni et al., 1996; Buizza, 1997). Following remarkable improvements in model technology and computer power, ECMWF took a pioneering step in June 2023 by upgrading its operational NWP system. The traditional EPS combined a high-resolution deterministic model with a lower-resolution ensemble forecast model. The new system, however, employs high-resolution models for both deterministic and ensemble forecasts, offering a more accurate representation of forecast uncertainty at higher resolutions. In 2006, the World Meteorological Organization (WMO) launched the Observing System Research and Predictability Experiment (THORPEX) Interactive Grand Global Ensemble (TIGGE) as an important component of the THORPEX initiative (Jiao, 2010). In recent years, ensemble forecasting has expanded to fields such as data assimilation, physical process parameterization, flood forecasting, severe weather prediction, and typhoon track forecasting. The ensemble forecasting approach has become an integral part of the entire operational NWP process (Buizza et al., 2018).

Chinese scientists have also been actively engaged in ensemble forecasting research and operations since the early stages (Li et al., 1997; Chen et al., 2002; Li and Chen, 2002). In terms of initial condition perturbation, methods such as conditional nonlinear optimal perturbation (CNOP; Mu et al., 2003; Duan and Huo, 2016), the heterogeneous physical mode method (Chen et al., 2005), nonlinear local Lyapunov vectors (NLLV; Feng et al., 2014), multi-scale hybrid initial perturbation (Zhang et al., 2015), and multi-scale SV (Ye et al., 2020) have been proposed. Regarding model perturbation, studies have focused on methods such as the multi-physics scheme (Chen et al., 2003), SPPT (Tan et al., 2013), nonlinear forcing SV (NFSV; Duan et al., 2013), SPP (Xu et al., 2019), and a model tendency perturbation method combining systematic and random errors (Han et al., 2023). China’s operational EPS has evolved from initially relying on imported spectral models to the self-developed GRAPES ensemble prediction system (GRAPES-EPS; Chen and Li, 2020). The initial perturbation methods have advanced from the ETKF (Ma et al., 2008; Long et al., 2011; Zhang et al., 2017) to the multi-scale SV initial perturbation method (Liu et al., 2011; Ye et al., 2020). Furthermore, a combination of SPPT and SKEB has been implemented to generate model perturbations. These advancements collectively represent a milestone in the development of China’s EPSs (Chen and Li, 2020; Shen et al., 2020).

This paper first reviews the ensemble weather forecasting methods and the primary techniques employed in the main EPSs designed by China and other countries, emphasizing the development of China Meteorological Administration (CMA) global and regional EPSs, including operational technologies, initial and model perturbation methods, and the applications of ensemble forecast. Besides, a comparison of the EPS settings and verification results between the CMA EPS and other major international NWP centers is presented. Furthermore, the paper addresses the key challenges in ensemble forecast technologies from the aspects of operation, science, integrated with Artificial Intelligence (AI), the integration of weather and climate, and user requirements. The summaries and future perspectives of ensemble forecast are concluded at the end.

2.   Review of ensemble forecast methods

2.1   Sources of initial errors and advances in initialperturbation methods

Initial errors play a critical role in numerical simulations, primarily arising from inaccuracies in observation, as well as errors introduced during the processes of data assimilation. The initial condition perturbation (ICP) method, which characterizes these initial errors, has been developed within global medium-range EPSs. The development of ICP methods has gone through several stages, from the initial human-generated perturbations to statistically-based perturbation techniques, such as the Monte Carlo, BGM, SV, ETKF, and CNOP. These methods have been widely applied in the construction of various EPSs, providing effective information about forecasting uncertainty and enhancing prediction skill.

2.1.1   Initial perturbation for estimating assimilationanalysis errors

The primary function of these methods is to estimate the error probability distribution in assimilated analyses. For example, the Monte Carlo-based random perturbation method superimposed random noise as perturbed field on initial condition (Hollingsworth, 1979); the time-lagged average method used previous analysis of interval time as perturbed initial condition (Hoffman and Kalnay, 1983); and the first version of Global Ensemble Prediction System developed by the Canadian Meteorological Centre (CMC) used the Monte Carlo-based random perturbations of observation to generate initial conditions (Houtekamer et al., 1996). With the development of data assimilation systems and EPSs, ensemble Kalman filter (EnKF) initial perturbation method, which randomlyperturbed observation data, is used by Canada’s EPS (Houtekamer and Mitchell, 2005) and the NCEP (Zhou et al., 2016). While these perturbation methods demonstrate advantages in effectively representing uncertainties in initial analysis errors and enabling generation of substantial ensemble members, they exhibit limitations including constrained growth of perturbation energy, insufficient ensemble spread, and inadequate characterization of forecast uncertainties. Recent innovations by Pan et al. (2021) and Wang et al. (2023) introduced an analysis-constrained method that effectively incorporates characteristics of assimilation analysis increments. This method improved initial perturbation quality and the consistency between ensemble spread and forecast error.

2.1.2   Estimation of initial perturbations withthe fastest error growth

These types of initial perturbations are based on research findings in atmospheric predictability. By analyzing the direction and speed of numerical prediction error growth in phase space, perturbations are introduced along the most unstable directions in initial conditions. Representative methods include the SV developed based on nonlinear dynamical finite-time instability theory (Molteni et al., 1996; Yang et al., 2002; Liu et al., 2011) and the BGM employed by NCEP, in which random perturbations are cycled through a nonlinear model and scaled to generate the final fields (Toth and Kalnay, 1993). Both the two methods aim to capture initial errors in atmospheric baroclinic instability regions, representing the forecast uncertainty of large-scale unstable waves. The SV, often developed in conjunction with a 4D variational system, has clear physical significance and can achieve favorable results. It has been successfully applied in medium-range weather ensemble forecasting by ECMWF, Japan, China, and others. Moreover, Zhang et al. (2020) proposed an initial perturbation method combining ensemble sensitivity analysis (ESA) and BGM to represent the characteristics of initial analysis errors. The distribution of perturbations, which is derived from sensitivity modes and rapidly growing perturbations calculated by BGM, adapts to changes in weather conditions and provides accurate simulations of the location and intensity of convective systems. The BGM method is favored by operational forecasting centers due to its effectiveness, simplicity, and low computational cost. Additionally, NCEP developed the ETR (Wei et al., 2008), a method similar to the BGM but with perturbation structures that are nearly orthogonal and may reduce the correlation among ensemble members.

The SV method has achieved significant success in operational ensemble forecasting. However, since it relies on a linear approximation of the nonlinear model, it fails to fully capture the impact of nonlinear physical processes, which limits further improvements in ensemble forecast skill. Mu et al. (2003) proposed the CNOP method, which comprehensively considers the effects of nonlinear processes and identifies the fastest-growing initial perturbations in nonlinear models, thereby addressing the limitations of the SV. Duan and Huo (2016) extended the CNOP to different perturbation phase spaces, developing an orthogonal scheme for ensemble forecast initial perturbations. This method has been successfully applied in the study of typhoon track ensemble forecasting, demonstrating that the initial perturbations generated by orthogonal CNOP significantly improve the typhoon track forecasting skill, especially for abnormal tracks, compared to those generated by the SV and the BGM methods (Zhang H. et al., 2023).

To address the rapid growth of errors caused by convective instability in meso- and small-scale systems, Chen et al. (2005) proposed a new method, known as the Heterogeneous Physical Mode Method, for constructing initial perturbations in ensemble forecasting that targets convective instability and exhibits mesoscale motion characteristics. By identifying forecast deviations from different cumulus convection parameterization schemes, this method identifies regions sensitive to convection, extracts perturbation variables, structures, and magnitudes, and constructs perturbed initial conditions. Unlike traditional initial perturbation methods such as BGM or SV, this approach perturbs the initial values in convectively unstable regions, thereby promoting the rapid growth of unstable perturbations related to convective motion.

2.1.3   Estimating initial perturbations for multi-scalesystem errors

With the continuous improvement of numerical weather prediction resolution, however, the rapid growth of multiscale initial errors continues to significantly impact forecast skill. Limited-area models face dual challenges, i.e., insufficient large-scale initial error information and inadequate characterization of small-scale errors from dynamical downscaling. In response to these issues, Chinese scientists developed the hybrid initial condition perturbation method for enhanced ensemble perturbation generation (Wang et al., 2014; Zhang et al., 2015; Zhuang et al., 2017; Ma et al., 2018). This method uses filtering and spectral analysis techniques to extract small-scale perturbations from regional ensemble forecasts and large-scale perturbations from global ensemble forecasts, and combines the two to generate the ICPs. The hybrid ICPs effectively represent both large- and small-scale uncertainties in the analysis, aligning more closely with the boundary perturbations provided by global EPSs, thereby significantly improving the forecast skill of regional EPSs. Furthermore, the CMA (Ye et al., 2020; Liu et al., 2024) and the Japan Meteorological Agency (JMA; Ono et al., 2021) have conducted research on multiscale SV ICP methods based on regional and global ensemble prediction models. These studies have generated ICPs containing multiscale initial uncertainty, which more effectively reflect the multiscale variations in initial errors and improve both regional and global ensemble prediction performance.

In general, ensemble forecasts based on ICPs formed by analysis errors tend to show slightly higher forecast accuracy but a lower spread, while ensemble forecasts based on initial perturbations derived from dynamical instability growth theory generally exhibit an inverse pattern, i.e., a slightly higher spread but marginally lower accuracy (Pauluis and Schumacher, 2013).

2.2   Model error sources and development of modelperturbation methods

Early ensemble forecasting studies primarily focused on initial perturbations. However, as research advanced, it became evident that addressing initial condition uncertainty alone could lead to drawbacks such as insufficient ensemble spread, resulting in systematically biased forecasts (Palmer et al., 2009). Therefore, it is necessary to consider the uncertainty caused by model errors or defects in ensemble forecasting. The model errors primarily originate from inadequate representations of both physical and dynamical processes in the numerical model (Mu et al., 2011). Currently, the predominant model errors considered by ensemble forecasting systems stem from imperfect representations of subgrid-scale parameterization processes.

2.2.1   Estimation of model uncertainty

The uncertainty of the physical process of the model is described by a combination of different physical parameterization schemes within the framework of a single model (Houtekamer et al., 1996), which is called a multiphysics scheme. This scheme was applied to the global EPS of Environment Canada in the early days and has also been commonly used in regional ensemble forecasting in recent years (Stensrud et al., 2000; Zhi et al., 2013) to improve the probabilistic forecasting skill. Another model perturbation method is the multi-model method (Krishnamurti et al., 1999), which can represent the systematic bias in the forecasts of different models, so that the ensemble members can produce a larger ensemble spread, which greatly improves the accuracy of medium-term forecasts and small- and medium-scale ensemble forecasts in limited areas (Zhang et al., 2017). However, since the ensemble members in multi-model or multi-physics schemes employ different models or physical processes, their physical consistency is compromised. This makes it challenging to satisfy the equal-likelihood requirement for members and increases the complexity of determining member weights for probability calculations (Berner et al., 2015). Consequently, the model perturbation method has been increasingly developed in recent years to estimate the uncertainty of the parameterization scheme of the model physical process.

2.2.2   Estimation of physical parameterization schemeuncertainty

Stochastic physical perturbation is now a standard technique for assessing physical parameterization uncertainty in modeling. The theoretical basis is to introduce a random process or factor to perturb some parameter values or related terms of a model (such as tendency and diffusion) to characterize the uncertainty of the model. This type of perturbation scheme mainly includes methods such as the stochastically perturbed physics parameterization tendency term (SPPT; Buizza et al., 1999; Yuan et al., 2016), stochastic kinetic energy backscattering (SKEB; Shutts, 2005), random parameters (RP), and stochastically perturbed parameterization (SPP; Chen et al., 2003; Bowler et al., 2008; Tan et al., 2013; Jankov et al., 2017).

The SPPT scheme employs spatiotemporally continuous random numbers conforming to a uniform distribution to multiplicatively perturb parameterized subgrid physical tendencies. By applying distinct perturbation values to different tendency terms, this method effectively enhances ensemble spread and improves probabilistic forecast reliability in ensemble forecasting systems (Li et al., 2008). The SPPT is currently the most widely adopted method among global NWP centers (Charron et al., 2010; Sanchez et al., 2016; Peng et al., 2020), yet it also has some drawbacks. For instance, perturbation amplitude may cause discontinuity and energy non-conservation in model top and near-surface fluxes, and the long-term integration will produce systematic bias (Leutbecher et al., 2017). Yuan et al. (2016) introduces SPPT technology to the CMA regional EPS, which significantly improves the precipitation forecast of heavy rainfall in the late medium-range forecast. Qiao et al. (2017) developed a novel parameterization scheme for stochastic perturbation dissipation based on the uncertainty of the model dissipation term. It follows the same procedure as the SPPT method, but uses a recursive filter to generate smooth perturbations. It also uses horizontal and vertical localization to maintain the influence of perturbation in the strong shear region, which effectively improves the systematically weak vortex intensity. Similar results have been demonstrated in perturbation studies of microphysical processes (Qiao et al., 2018).

The SKEB scheme targets excessive momentum dissipation errors in the numerical model’s dynamic framework at truncated scales, compensating for dissipated energy effects on resolvable scales through stochastic perturbations. It has been applied in global medium-range ensemble forecasting (Berner et al., 2009; Zhou et al., 2017) and is being progressively implemented in regional convective ensemble forecasting (Cai et al., 2017; Yang et al., 2024). Peng et al. (2020) found that the combined application of SKEB and SPPT in the GRAPES global ensemble forecast can effectively improve the ensemble spread in the tropics. The SKEB scheme can excite possible instabilities in atmospheric motion (Berner et al., 2015), thereby producing larger ensemble spread and capturing low-probability events. However, due to the higher computational cost of the SKEB compared to the SPPT and the unsatisfactory balance between spread and root mean square error (RMSE), the value of the SKEB scheme has gradually decreased. It is reported that the ECMWF has discontinued the operational implementation of this scheme (Chen et al., 2021).

The RP/SPP method characterizes model uncertainty by performing random perturbations of subjective empirical and semi-empirical parameters within physical parameterization schemes (Bowler et al., 2008), which represents the predicted uncertainty of small-scale system changes. Although this method improves the ensemble forecasting, it has certain irrationalities. For instance, the parameter values will not change suddenly at a certain time (Jankov et al., 2019). In addition, it is found that the SPP method has a more significant effect on the parameters in cumulus convection and planetary boundary layer schemes (Xu et al., 2019). However, the RP/SPP method yields inadequate ensemble spread for medium- and extended-range forecasts (Leutbecher et al., 2017).

Beyond these three prevalent stochastic physical perturbation schemes, researchers have developed various approaches to address uncertainty sources in other model components, such as the stochastic trigger of convection (STC; Li et al., 2015), stochastic boundary-layer humidity (SHUM; Whitaker et al., 2013), vorticity confinement (VC; Steinhoff and Underhill, 1994), meso- and medium-scale terrain perturbation schemes (Li et al., 2017), and model tendency perturbation methods combining model systematic bias and random errors (Han et al., 2023). Zhao and Torn (2022) applied an independent stochastic perturbation parameterization scheme to a single parametric perturbation, and the results showed that turbulent mixed stochastic perturbation can increase the ensemble standard deviation of typhoon intensity, while the impact of stochastic perturbation on microphysics, radiation, and cumulus tendencies is negligible. The SHUMscheme improves the consistency between ensemble spread and forecast error in the tropical regions, and reduces the ensemble mean forecast error. The SHUM is applied to the ensemble data assimilation module of the NCEP global forecast system EnKF/3DVar. The VC method is applied to the NCEP global EPS, which increases the subtropical ensemble spread. The model tendency perturbation method combining model systematic bias and random error reduces the systematic bias of the model, which is applied to the CMA global medium-term EPS.

2.2.3   Nonlinear characteristics of model errorgrowth

Accurately characterizing the nonlinear nature of model error growth remains a critical challenge for cutting-edge high-resolution ensemble forecasting systems. Wang et al. (2020) implemented the CNOP perturbation (CNOP-P) method in the GRAPES model for convection-scale ensemble forecasting. The experiments demonstrated enhanced ensemble spread in tropospheric humidity and temperature predictions, along with improved forecast reliability for near-surface variables and precipitation. Xu et al. (2022a) applied a nonlinear forcing SV method (NFSV), also known as CNOP forcing (CNOP-F). Duan and Zhou (2013) described the combined effects of model errors from different sources by superimposing NFSV perturbations with specific structures on top of SPPT perturbations, which better characterizes the impact of model uncertainty on convection-scale ensemble forecasting (Xu et al., 2022b). In addition, Zhang Y. C. et al. (2023) extended NFSV to the orthogonal subspace of inclined perturbations and developed an orthogonal NFSV ensemble prediction model perturbation method. When applied to typhoon intensity ensemble prediction, this method demonstrated significantly higher forecasting skill than SPPT and SKEB, particularly in providing earlier warning information for predicting typhoon rapid intensification processes.

In summary, the multi-model and multi-physics schemes have been gradually eliminated due to the average error of members, the inequality of meteorological significance, and the complexity of probability calculation, while the stochastic physics perturbation method has received more and more attention. However, whether the representativeness and flow dependence of reasonable stochastic functions and their spatiotemporal correlation scales of random noises truly reflect the uncertainty of physical processes still need to be studied.

2.3   Post-processing and verification assessment ofensemble forecast uncertainty information

The EPS generates a vast array of data products. In contrast to a single deterministic forecast, ensemble forecasting can provide a probability distribution of possible future atmospheric states and characteristics of forecast uncertainty, demonstrating superior advantages in the early warning of extreme weather events (Du and Chen, 2010; Gao et al., 2019). The critical task in the post-processing of ensemble forecasts is to design appropriate post-processing techniques that provide meteorologically meaningful interpretations of the distribution of ensemble forecast products, as well as to extract similarities, discrepancies, and extreme information among the ensemble forecast members (Williams et al., 2014). The primary techniques currently employed in ensemble forecast post-processing include obtaining probabilistic forecast information to describe forecast uncertainty, correcting systematic model errors, and extracting extreme weather forecast information from the tails of the ensemble forecast probability distribution.

2.3.1   Probability forecast information based onensemble prediction

The ensemble mean and ensemble spread are two fundamental products of ensemble forecasting. The ensemble mean represents the primary level of information in ensemble forecast products. It filters out unpredictable elements from the ensemble members, providing an overall forecast trend. However, due to its inherent smoothing effect, the ensemble mean might not completely retain extreme values in weather events. The ensemble spread is a measure of the uncertainty in ensemble forecasting or the amplitude of variation of ensemble members relative to the ensemble mean. It can be measured by the standard deviation of the ensemble members relative to the ensemble mean forecast, or by the average anomaly correlation coefficient (ACC) of each member relative to the overall mean field. The spread, to some extent, can represent the skill of ensemble forecasting. Generally speaking, a smaller spread indicates higher forecast confidence with potentially higher forecast skill; however, a larger spread indicates lower forecast confidence but does not necessarily imply lower forecast skill.

Probabilistic forecast of weather elements is a crucial ensemble forecast product. By calculating the forecast probabilities for different thresholds of weather elements such as precipitation, temperature, and wind, it represents the occurrence likelihood of specific weather conditions. For instance, precipitation can be categorized into thresholds such as greater than 1, 5, 10, and 20 mm day⁻¹, or the forecast probabilities for specific weather elements at certain stations can be provided. Another probabilistic forecast product is the spaghetti plot, which involves selecting a specific characteristic contour line and plotting forecasted contours from all ensemble members on the same chart. Generally, the degree of divergence among contour lines roughly indicates forecast confidence, and more concentrated contours correspond to higher confidence.

Ensemble forecast clustering products utilize clustering analysis to classify ensemble member forecasts into distinct clusters, offering various possible outputs. Common clustering analysis algorithms include hierarchical, non-hierarchical, and tubing clustering approaches. Through clustering classification, the number of members in each cluster can be compared, allowing for the determination of the most probable weather scenarios. Additionally, two commonly used visualization tools, i.e., plume chart and spaghetti chart, apply clustering analysis principles to classify 500-hPa geopotential height (Z500) patterns, providing forecasters with operationally actionable guidance.

2.3.2   Post-processing methods for systematic errorcorrection in ensemble forecasting

Common post-processing methods for systematic bias in ensemble forecasting can be broadly categorized into parametric and non-parametric approaches (Mylne et al., 2022). Parametric post-processing methods generally assume that the target variable follows a specific probability distribution and conduct regression on its parameters. Examples include ensemble model output statistics (EMOS), Bayesian model averaging (BMA), and Bayesian joint probability (BJP). Additionally, correction methods based on different probability distributions have been developed for various physical quantities. For instance, Gaussian distribution is used for temperature and pressure, truncated normal distribution and lognormal distribution for wind, and Gamma (and truncated Gamma) distribution for precipitation. Furthermore, ensemble forecast correction based on precipitation categories has been proposed and applied to improve the forecast skill for different levels of precipitation, particularly extreme precipitation (Ji et al., 2023).

Non-parametric post-processing methods eliminate the need for predefined probability distributions of target variables, thereby providing enhanced flexibility. Common non-parametric methods primarily include quantile mapping (QM), frequency matching (FM), probability matching (PM), member by member (MBM), the optimal percentile method, and isotonic distributional regression (IDR).

Meanwhile, by combining effective forecast information obtained from multiple models for multi-model ensemble forecasting, non-equal-weight ensemble methods based on error analysis, such as super ensemble and Kalman filtering, have been developed (Zhi et al., 2013; Zhu et al., 2021). The concept of a sliding training period (Zhi et al., 2012) has been widely applied in the correction of forecasts for tropical cyclones, temperature, precipitation, wind, and other variables. To reduce the positional and structural biases in the target system forecasts, various post-processing methods based on spatial features have been proposed, such as neighborhood methods, modal projection methods, and object-based correction techniques (Ji et al., 2020).

In recent years, machine learning techniques have been widely adopted due to their high nonlinearity and strong robustness. In recent years, deep learning algorithms such as random forests, support vector machines, convolutional neural networks, long short-term memory neural networks, and U-Net neural networks have also been extensively applied in the post-processing of EPSs, achieving remarkable results (Zhi et al., 2020; Yang et al., 2022; Lyu et al., 2023). At the same time, by integrating machine learning algorithms with traditional statistical methods, hybrid machine learning and statistical approaches have been designed to further enhance the reliability of ensemble forecasting. Detailed descriptions of the aforementioned systematic error post-processing methods can be found in the ensemble forecast guidance manual published by the WMO (Mylne et al., 2022).

2.3.3   Extreme forecast information extractionbased on ensembles

To extract early warning information for extreme weather from ensemble forecasts, the Extreme Forecast Index (EFI) developed through ECMWF ensemble forecast results was established (Lalaurette, 2003). This index characterizes the continuous difference between the cumulative probability distribution of the ensemble forecast results and the model climate cumulative probability distribution. The larger the value, the greater the deviation of the forecast from the model climate. Such deviations indicate a higher probability of extreme forecasts, and consequently, an increased probability of extreme weather occurring in reality. NCEP has also conducted studies on extreme weather forecasting based on EFI. Guan and Zhu (2017) pointed out that the EFI index based on NCEP global ensemble forecasts has certain medium-range forecasting ability for extreme precipitation and extreme low-temperature weather in the United States during the winter of 2013/14. In addition to conventional weather elements such as temperature, precipitation, and 10-m wind, the EFI index can also be applied to diagnostic quantities such as convective parameters (e.g., convective available potential energy) and total column water vapor flux (Tsonevsky et al., 2018). The application range of the EFI index has gradually expanded from large-scale medium-range ensemble forecasts to convective-scale and sub-seasonal to seasonal ensemble forecasts (Dutra et al., 2013; Raynaud et al., 2018). Many domestic scholars have also studied the EFI index (Xia and Chen, 2012; Liu et al., 2018; Peng et al., 2022). Research on extreme weather forecasting methods based on EFI using the CMA’s global ensemble forecasts has shown that EFI has certain recognition ability for extreme low temperatures and extreme heavy precipitation in China, with good medium-range forecasting skill. Another similar method is the Ensemble Anomaly Forecasting by ranking ensemble forecasts relative to observed climatology rather than model climatology (Du et al., 2014).

2.3.4   Verification and evaluation metrics

The accuracy and reliability of ensemble forecasting are crucial in meteorological prediction, driving the continuous development and refinement of verification and evaluation methods for ensemble forecasts. The works of Wilks (2011) and Jolliffe and Stephenson (2012) comprehensively summarize the verification methods for meteorological forecasts, covering deterministic, probabilistic, and qualitative forecasts, along with their applications. For deterministic forecasts based on ensemble means or individual ensemble members, commonly used verification metrics include the RMSE and spatial correlation coefficients. In the realm of probabilistic forecast verification, primary evaluation methods encompass the Brier skill score (BSS), continuous rank probability skill score (CRPSS; Hersbach, 2000), reliability diagrams, relative operating characteristics (ROC) curve, and potential economic value (Zhu et al., 2002).

In quantifying uncertainty, ensemble spread and the Talagrand distribution are key metrics for assessing the uncertainty of forecast systems. The spread–error relationship, commonly referred to as the spread–skill relationship, is often employed to evaluate the predictability of EPSs. However, real-world systems often exhibit under-dispersion in ensemble spread, particularly in quantitative precipitation forecasting (Li, 2001; Li et al., 2009; Su et al., 2014). Uncertainty quantification typically relies on the bootstrapping method (Hamill, 1999), with error bars representing the confidence intervals of the scores. Additionally, Torn and Hakim (2008) introduced the ESA method to assess the sensitivity of forecast outcomes to minor changes in initial conditions or other input parameters.

Su et al. (2014) proposed an area-weighted verification method, applying spatial weighting to probabilistic precipitation forecast scores. In the field of precipitation forecasting, some researchers have integrated these spatial verification methods into ensemble forecast evaluation. For instance, Ji et al. (2020) explored the integration of the Method for Object-based Diagnostic Evaluation (MODE) with ensemble forecasting techniques, finding significant advantages in evaluating and improving the spatial structural characteristics of precipitation forecasts. Chen et al. (2018) investigated the spatial relationship between ensemble spread and forecast error for Meiyu precipitation forecasts in the Yangtze–Huaihe region using the precipitation consistency scale method, revealing the impact of positional errors on forecast skill under different precipitation thresholds. Furthermore, Du and Deng (2020) introduced two new forecast evaluation metrics, i.e., the measure of forecast challenge (MFC) and the predictability horizon diagram index (PHDX). The MFC, which may comprehensively consider forecast error and uncertainty, serves as a novel metric for assessing forecast difficulty. The PHDX further investigates how the temporal evolution of forecasts influences decision-making processes, aiming to better represent the dynamic nature of ensemble forecast information.

It is important to note that the accuracy of scores in ensemble forecast verification may be influenced by reference values. For example, the calculation of BSS by using high-frequency grid-based climatological samples could potentially lead to score underestimation (Hamill and Juras, 2006). Additionally, the uncertainty and quality of verification data can significantly impact evaluation outcomes, particularly in precipitation forecast verification (Yuan et al., 2005). With the ongoing advancement of high-resolution ensemble forecasting and observational data, ensemble forecast verification methods are expected to become more refined and scientifically robust in the future.

3.   Chinese operational ensemble prediction system and its applications

3.1   Development history of global/regional ensembleprediction systems

3.1.1   Global ensemble prediction system ofspectral model

In the 1990s, the National Meteorological Center (NMC) of CMA developed the global EPS (Li et al., 1997). In 1996, a T63L16 EPS had been implemented using the time-lagged method, consisting of 12 ensemble members with a 10-day forecast lead time. In 1999, a T106L19 EPS was developed by the SV method (Li and Chen, 2002), which includes 32 ensemble members with a 10-day forecast lead time. In 2007, based on the BGM method (Tian et al., 2007), a T213 EPS was established, which includes 15 ensemble members with a 10-day forecast lead time. The year 2014 saw the upgrade of CMA’s operational EPS from T213 to T639, which for the first time incorporated SPPT to generate model perturbations (Tan et al., 2013).

Exploratory efforts toward regional EPS commenced in the early 21st century. In 2005, based on the mesoscale model MM5, a heterogeneous physical mode method was proposed to generate initial perturbation (Chen et al., 2005), alongside experimental studies incorporating multi-physics configurations and stochastic physical parameter perturbations (Chen et al., 2003; Feng et al., 2006; Wang et al., 2007). By 2006, the NMC developed a regional mesoscale EPS over China using the Weather Research and Forecasting (WRF) model, integrating the BGM method and multi-physics ensemble method (Deng et al., 2010). This system became operational EPS on 15 November 2010, operating at a horizontal resolution of 15 km with 15 ensemble members and running 4 daily forecast cycles (Ma et al., 2008).

3.1.2   GRAPES EPS

The GRAPES model is a new-generation numerical weather prediction system developed through self-innovation by the CMA (Chen et al., 2008). Since 2005, research and development efforts have been dedicated to advancing ensemble forecast techniques within the GRAPES framework. The Chinese Academy of Meteorological Sciences has successively implemented both the BGM and ETKF methods to investigate GRAPES regional ensemble prediction techniques and systems (Tan and Chen, 2007; Tian and Zhuang, 2008; Ji et al., 2011), ultimately developing its regional EPS (GRAPES-REPS). The model’s version 1.0 became operational in 2014 at NMC, utilizing ETKF for initial condition perturbations and superseding the WRF regional EPS. This system featured a 15-km horizontal resolution with 15 ensemble members. In 2015, the initial condition perturbation scheme was upgraded to a multi-scale blending (MSB) method (Zhang et al., 2015) that combined large-scale initial perturbations derived from the T639 EPS by using the BGM method with meso-scale initial perturbations obtained through the ETKF method within the GRAPES-REPS. This upgrade culminated in the release of version 2.0. In September 2019, the GRAPES global medium-range EPS, employing the SV method for initial condition perturbations, was operationally implemented (Liu et al., 2011) to supersede the T639 EPS. The lateral boundary perturbations for the new GRAPES-REPS were sourced from the GRAPES EPS. The parallel enhancement of initial condition perturbation and model perturbation methods led to the operational system’s evolution to version 3.0 (Chen and Li, 2020).

The development of the GRAPES global EPS began in 2008, with initial research focused on the global ETKF method (Ma et al., 2008). With the advancement of the GRAPES global 4DVar data assimilation technology, the GRAPES SV ICP technique was developed (Li and Liu, 2019). In December 2018, the GRAPES global medium-range EPS, based on the SV ICP, became operational and replaced the T639 EPS. This achievement represents China’s first successful implementation of both global and regional operational EPSs utilizing fully independent and self-controlled technologies (Chen and Li, 2020).

To characterize uncertainties in the GRAPES model, researchers have investigated multiple model perturbation methods through GRAPES regional and global EPSs, building upon prior studies of the NMC’s T213 and T639 global medium-range ensemble models. Employed methods include the multi-physics (MP; Feng et al., 2006), SPPT (Yuan et al., 2016), SKEB (Peng et al., 2019), and SPP (Xu et al., 2019). The SPPT scheme became operational in the GRAPES regional EPS in 2015, followed by implementation of both SPPT and SKEB in the global EPS in 2018 (Chen and Li, 2020).

Since 2018, significant methodological improvements have been achieved in the tropical cyclone SV technique, blended SV-EDA perturbation scheme, and sea surface temperature perturbation scheme (Huo et al., 2020; Qi et al., 2022). Based on the GRAPES-REPS regional forecast model, extensive advances in initial condition and model perturbations include stochastic parameter perturbation (Xu et al., 2019), conditional typhoon vortex relocation (Wu et al., 2020), radar reflectivity algorithm (Chen et al., 2021), and model tendency perturbation combined with systematic bias and random errors (Han et al., 2023).

Since 2020, active research on convective-scale ensemble forecasting has been conducted using the GRAPES-Meso 3-km model. Researchers have performed various experiments involving both model and initial condition perturbations, including multi-scale initial condition perturbation experiments (Ma et al., 2023) and the construction of model perturbations through conditional nonlinear optimal perturbations combined with stochastic physical perturbations (Xu et al., 2022a, b). The established 3-km regional ensemble forecast experimental system has been applied in major events such as the Beijing Winter Olympics and the Hangzhou Asian Games. Under the WMO Research Demonstration Project for the Hangzhou Asian Games, a fully independent 3-km convective-scale ensemble prediction system was successfully developed and implemented in Hangzhou, providing ensemble forecast services for the broader East China region. This system commenced operational service on 21 May 2023, delivering comprehensive probabilistic forecasts and multifaceted decision-support capabilities to forecasters.

3.2   Technology of CMA EPS

3.2.1   Global/regional integrated multi-scale SVinitial perturbation technique

The SV initial perturbations effectively capture both the primary error information and the sensitive features of analysis errors that grow and intensify within the basic flow. The development of the Tangent Linear Model (TLM) and Adjoint Model (ADM) within the CMA global 4DVar system has provided the necessary foundation for the research and advancement of SV perturbation techniques in the CMA global EPS (Liu et al., 2017). The SV initial perturbation calculation scheme for the CMA global EPS is composed of two parts, i.e., calculations of SVs for the extratropical regions and for tropical cyclones, respectively. The specific mathematical treatment scheme is outlined as follows.

3.2.1.1   SV calculation scheme for extratropical regions

In the CMA global EPS, the SV calculation can be mathematically formulated as maximizing the ratio of the evolved perturbation vector norm to the initial perturbation vector norm, as shown in Eq. (1):

where L represents TLM in the CMA global 4DVar assimilation system, L^T indicates the corresponding adjoint model, E is the weight matrix that measures the perturbation magnitude, x denotes the SV, λ is the corresponding singular value, and P stands for a projection operator that sets the SV perturbations outside the target region to zero. From Eq. (1), it is evident that the SV structure is influenced by three key factors.

(1) Definition of the weight matrix E. This matrix determines the relative importance of different components of the perturbation vector.

(2) Characteristics of the TLM (L) and ADM (L^T). These models primarily involve the use of linearized physical process parameterization schemes.

(3) Integration of time length for the TLM (L) forward and the ADM (L^T) backward, referred to as the optimal time interval.

The CMA global SVs employ the dry total energy norm to define the weight matrix E (Liu et al., 2013). Let the prognostic variables of the CMA global TLM and ADM include the horizontal wind components (i.e., u and v), perturbation potential temperature (${\theta }'$), and perturbation dimensionless pressure (${\Pi }'$). The corresponding perturbation quantities can be expressed as ${u}'$, ${v}'$, $({\theta }'{)}'$, and $({\Pi }'{)}'$. The total energy norm E is calculated as follows:

In Eq. (2), the sums of the first two and last two terms represent the perturbation kinetic energy (KE) and the perturbation potential energy (PE) norms, respectively. The third and fourth terms correspond to the contributions of perturbation potential temperature $({\theta }'{)}'$ and perturbation dimensionless pressure $({\Pi }'{)}'$ to the PE norm. Meanwhile, $dV=d\lambda d\varphi d\widehat{z}$, among which $\widehat{z}$ denotes the terrain-following coordinate, and $\lambda$ and $\varphi$ represent longitude and latitude in the spherical coordinate system of the model. In addition, ${C}_{p}$ is the specific heat capacity of dry air at constant pressure; and ${T}_{r}$, ${\theta }_{r}$, ${\Pi }_{r}$, and ${\rho }_{r}$ denote the reference temperature, reference potential temperature, reference dimensionless pressure, and reference density, respectively.

Based on the aforementioned CMA global SV calculation technique, three types of SVs with different scales are defined by varying resolutions and optimal time intervals. Large-scale SV (LSV) has a horizontal resolution of 2.5° and an optimal time interval of 48 h. It is computed by using dry linearized physical processes and is primarily used to capture uncertainty information at the synoptic scale. Meso-scale SV (MSV), with a horizontal resolution of 1.5° and an optimal time interval of 24 h, also employs dry linearized physical processes and is designed to capture uncertainty at the meso-α scale. Small-scale SV (SSV), which has a horizontal resolution of 0.5° and an optimal time interval of 6 h, incorporates moist linearized physical processes in its computation compared to LSV and MSV, and is mainly used to capture uncertainty at the meso-β scale. The specific configurations of the CMA global/regional integrated multi-scale SVs are summarized in Table 1.

Table  1.  The settings of calculation for multi-scale SVs

Download as CSV

	LSV	MSV	SSV
Target region	30°–80°N, 30°–80°S	30°–80°N, 30°–80°S	20°–50°N, 105°–125°E
TLM resolution/optimization time	2.5°/48 h	1.5°/24 h	0.5°/6 h
Energy norm	Dry energy	Dry energy	Dry energy
Linearized physical process	Subgrid topographic drag, vertical diffusion	Subgrid topographic drag, vertical diffusion	Subgrid topographic drag, vertical diffusion, large-scale condensation
Number of SVs	10	10	10

 | Show Table

DownLoad: CSV

3.2.1.2   Tropical cyclone SV calculation scheme

Considering the unique characteristics of tropical cyclones, a specialized tropical cyclone SV calculation scheme has been developed for the CMA global EPS (Huo et al., 2020). Specifically, for tropical cyclone cases, the target region for SV computation is defined as a 10°×10° (latitude × longitude) domain centered on the cyclone position. A maximum of six tropical cyclone SVs can be calculated.

As a result, the target regions for SV computation in the CMA global EPS consist of the extratropical regions of the Northern and Southern Hemispheres and the tropical cyclone region. This approach not only captures the development of baroclinically unstable perturbations in the mid-to-high latitudes of both hemispheres but also accounts for the evolution of tropical cyclones when they are present.

Based on the initial SV computed for the Northern and Southern Hemispheres as well as tropical cyclones, Gaussian sampling techniques are employed to construct the initial perturbation fields for the CMA global EPS. Additionally, a scaling amplification factor $\sigma$ is applied to normalize the amplitudes of all initial perturbations, ensuring that their magnitudes are consistent with the actual analysis error levels. Finally, by adding and subtracting the initial perturbations to/from the assimilation analysis initial conditions in pairs, the perturbed initial conditions for the CMA global EPS are generated.

3.2.2   MSB initial perturbation techniques inshort-range regional ensembles

Due to the resolution limitations of global EPSs, the dynamical downscaling method fails to fully capture uncertainties at smaller scales that can be resolved by regional models. However, meso- and small-scale initial perturbations are essential for accurately capturing uncertainties in local severe weather. Therefore, the CMA regional EPSs have developed an MSB initial perturbation technique based on the ETKF scheme.

3.2.2.1   ETKF initial perturbation method

The ETKF method is an initial perturbation scheme developed based on Kalman filter theory. It rapidly estimates analysis errors from ensemble forecast perturbations and observational error variances (Long et al., 2011), thereby effectively capturing uncertainties in meso- and small-scale initial conditions.

In the CMA regional EPS, the forecast perturbation vector ${\mathit{X}}^{\mathit{f}}$ is transformed into the analysis perturbation vector ${\mathit{X}}^{\mathit{a}}$ through the transformation matrix T, as shown in Eq. (3). Specifically, the analysis perturbations at the current time are derived by by linearly combining the forecast perturbation vector using the transformation matrix T.

The core of the ETKF initial perturbation method lies in obtaining the transformation matrix T. In the ETKF scheme of CMA regional EPS, the transformation matrix T is defined using the method proposed by Wang and Bishop (2003). To centralize the perturbed members relative to the ensemble mean, we employ a spherical simplex centralization scheme. In the operational version of the CMA regional EPS, simulated observations are derived by interpolating the analysis fields from the global EPS to the observational space (Wang et al., 2018). The variables used to compute the transformation matrix T include meridional wind v and zonal wind u, and the variables for calculating the amplification factor are u, v, and specific humidity q. The system applies a 6-h cycling perturbation scheme to provide 6-h forecast perturbation fields for the ETKF cycle. The model runs 4 times per day at 0000, 0600, 1200, and 1800 UTC, among which the 0000 and 1200 UTC cycles are integrated for 84-h forecasts, whereas the 0600 and 1800 UTC cycles are for 6-h forecasts.

3.2.2.2   MSB initial perturbation method

Practical experience has shown that the initial perturbation fields in regional ensemble prediction, derived from the dynamical downscaling of global ensemble prediction, contain more large-scale perturbation information. Since 2016, the CMA regional EPS has adopted an MSB initial perturbation technique (Zhang et al., 2015; Xia et al., 2019). This method employs the two-dimensional discrete cosine transform (2D-DCT; Denis et al., 2002) to filter the initial perturbations from both the ETKF method in regional ensemble prediction and the dynamical downscaling of the T639 global ensemble prediction. Subsequently, through the 2D-DCT inverse transformation, the initial perturbation fields are reconstructed by integrating large-scale perturbations from the global ensemble prediction with meso- and small-scale perturbations from the regional ensemble prediction, with more multiscale initial condition uncertainties.

3.2.2.3   Conditional typhoon vortex relocation techniquein ensemble forecasting

To more reasonably describe the uncertainty in the positioning of typhoon vortex centers in ensemble forecasting, Wu et al. (2020) applied the best track data of typhoons from the CMA and the JMA during 2009 and 2018. They analyzed the uncertainties of typhoon vortex center positioning in these best tracks and developed a conditional typhoon vortex relocation method. Additionally, they established procedures for threshold discrimination of typhoon vortex relocation for ensemble members, as well as mathematical treatment processes for typhoon vortex separation and vortex relocation.

By analyzing and comparing the best typhoon vortex center locations from the CMA and the JMA during the study period, it was found that the maximum annual average positioning error was 17.18 km (in 2009), while the minimum was 11.98 km (in 2015). Given the 10-km horizontal resolution of the CMA-REPS, a critical threshold of 15 km was established to determine whether conditional typhoon vortex relocation should be performed for ensemble members. Specifically, if the distance between the typhoon vortex center in an ensemble member’s initial field and the best track exceeds 15 km, the typhoon vortex positioning in that member’s initial condition is deemed beyond the range of analyzed uncertainties, requiring conditional vortex relocation. Otherwise, the error is considered within the range of analyzed uncertainties, making relocation unnecessary. The relocation method utilizes the typhoon vortex separation technique developed by the U.S. Geophysical Fluid Dynamics Laboratory to isolate the typhoon vortex from the initial fields of ensemble members requiring relocation. The separated typhoon is then shifted to the observed position through interpolation. This scheme was implemented in the CMA regional EPS v3.0 in 2019.

3.2.3   Stochastic model perturbation technique forsubgrid physical processes

The EPS relying solely on initial perturbations suffers from insufficient divergence among ensemble members and inadequate reliability. In a complete system, incorporating model perturbation techniques is essential to effectively characterize model uncertainty and account for model errors (Palmer et al., 2009). Therefore, CMA’s global/regional EPS should introduce such techniques to represent uncertainties arising from model imperfections. At present, SPPT and SKEB model perturbation techniques are mainly used in the CMA’s forecast system, which are briefly introduced below.

3.2.3.1   SPPT scheme

The SPPT scheme in CMA’s global/regional EPS shares the same fundamental design principles as ECMWF’simplementation (Buizza et al., 1999). Specifically, during numerical integration, the model integration term can be decomposed into the integral term of non-parametric (dynamic) processes, and the integral tendency term of parametric physical processes. The model integral tendency term perturbed by the SPPT scheme can then be expressed as:

where $\varPsi (\lambda ,\phi ,t)$ is a three-dimensional random field with spatiotemporal correlation features. A stretch function is introduced to implement a range of custom perturbations (Li et al., 2008), which is defined as:

where β is a constant with a value of −1.27 and $\mu =({\varPsi }_{\rm max}+{\varPsi}_{\rm min})/2$; $\varPsi_{\rm max}$ and ${\varPsi }_{\rm min}$ represent the upper and lower bounds of the random field $\varPsi \left(\lambda ,\phi ,t\right)$, respectively; and $\psi (\lambda ,\phi ,t)$ is a three-dimensional random field with spatiotemporal correlated characteristics, defined as:

where $\lambda$, $\phi$, and t represent the longitude, latitude, and time, respectively; and ${Y}_{l,m}$ is the spherical harmonic function, L is the horizontal truncation scale of the random field, l and m are the total wavenumber and latitude wavenumber in the horizontal direction, and ${\alpha }_{l,m}\left(t\right)$ is the spectral coefficient of the random field. The correlation of the stochastic spectral coefficient ${\alpha }_{l,m}\left(t\right)$ in the time dimension is realized by the first-order Markov chain stochastic process (also known as the first-order autoregressive stochastic process), as shown in the following formula:

where t is the specific time interval that can correspond to the integration step of the model in the CMA’s EPS; τ is the time scale of the random field decorrelation; and ${R}_{l,m}\left(t\right)$ is a Gaussian distribution stochastic process that obeys a variance of 1 and a mean of 0. The stochastic type constructed by the above formula features time–space scale correlation and controllable perturbation amount. In the CMA global/regional EPSs, perturbations are applied only to the net tendency terms of potential temperature, horizontal wind components, and humidity variables.

Considering the different integration times and forecast objects of global and regional models, the parameter settings of the SPPT scheme in the CMA global and regional EPSs are different. In the CMA regional EPS, the SPPT perturbation amplitude ranges from 0.2 to 1.8 with a mean of 1.0. The scheme uses a decorrelation time scale of 6 h, a maximum wavenumber of 24, and a standard deviation σ of 0.27. In the CMA global EPS, to ensure operational stability, the stochastic function’s perturbation amplitude is constrained to 0.5–1.5. The vertical profile of tendency perturbations is introduced at both the near-surface layer and atmospheric top, where physical tendencies either remain unperturbed or receive only minimal-amplitude perturbations. The SPPT scheme can significantly improve the ensemble spread and missing rate of temperature and wind speed forecasts, improve the probabilistic forecasting skills of heavy precipitation, and significantly increase the ensemble spread of various elements (altitude field, temperature field, and wind speed) forecasts in tropical areas, which makes up for the lack of perturbation in tropical areas by SV initial perturbation technique (Peng et al., 2020).

3.2.3.2   SKEB scheme

The main purpose of the SKEB scheme is to represent the stochastic process and uncertainty of subgrid-scale energy upscaling transitions in the model (Shutts, 2005). The SKEB scheme in the CMA global EPS addresses excessive energy dissipation by employing a random flow function forcing with specific spatiotemporal correlations and local dynamic dissipation rates (Peng et al., 2019). The random stream function forcing ${F}_{\psi }$ is defined as follows.

In Eq. (8), $\varPsi \left(\lambda ,\varphi ,t\right)$ is a random type, and its generation method is the same as the definition of the random type of SPPT, as shown in Eq. (5); and $D\left(\lambda ,\varphi ,\eta ,t\right)$ is the dissipation rate of local dynamic energy. At present, the SKEB scheme in the CMA global EPS constructs the local dynamic energy dissipation rate based on the explicit horizontal scheme.

where k is a constant factor greater than 1, $\boldsymbol{u}$ is the horizontal wind speed, and ${\boldsymbol{u}}'$ is the change in the horizontal wind speed before and after the application of the horizontal diffusion scheme. Based on the relationship between the flow function and the rotational component of the horizontal wind field, the horizontal wind perturbation forcing terms of the CMA global EPS are constructed as ${S}_{u}$ and ${S}_{v}$. These stochastic forcing terms are then added to the tendency terms of the model’s horizontal wind field.

In the CMA global EPS, the minimum and maximum truncation wavenumbers L_min and L_max of the $\varPsi$ random type of SKEB scheme are 10 and 80, the decorrelation time scale τ is 6 h, the mean μ is 0, the standard deviation σ is 0.27, and the maximum ${\varPsi }_{\rm max}$ (minimum ${\varPsi }_{\rm min}$) is 0.8 (−0.8). The application of the SKEB scheme can improve the simulation ability of the CMA global EPS on the atmospheric kinetic energy spectrum, and significantly enhance the relationship between the ensemble mean error and the ensemble spread of the wind field forecast in the tropics.

3.2.4   Systematic bias and extreme informationextraction

In addition to the aforementioned initial perturbation and model perturbation techniques, the CMA global/regional EPSs have innovatively developed a suite of ensemble post-processing techniques and extreme weather forecasting applications. These include an ensemble dynamical method for correcting bias, the radar reflectivity for subgrid-scale precipitation, and the model of extreme weather index. These techniques addressed the challengesof extracting extreme forecast signals from the tails of probabilistic forecast distributions and correcting systematic model biases, leading to notable improvements in probability density distributions and enhanced extreme weather prediction capabilities in CMA-GEPS/REPS. A brief overview of these technologies is as follows.

3.2.4.1   Ensemble forecasting dynamic correction methodbased on subtracting model systematic bias intendency term

Current ensemble prediction techniques primarily address model stochastic errors. However, under pronounced model systematic biases, sole dependence on either initial perturbation or model perturbation methods is inadequate to optimize the error–spread relationship in ensemble forecasts. Following the research on the impacts of CMA global/regional model systematic biases on ensemble probability density distributions, an ensemble forecasting method was developed, applying dynamic correction through subtraction of model systematic bias in the tendency term. This approach removes bias tendencies from both dynamical and physical tendency terms during model integration, as shown in Eq. (12) (Chen et al., 2020; Han et al., 2023).

In the equation, ${e}_{j}\left(t\right)$ is the total tendency, $A({e}_{j},t)$ is the dynamic tendency, $P({e}_{j},t)$ is a physical tendency, and $\stackrel{\wedge }{{B}_{l}}\left({e}_{0}\right)$ is the subtracted bias. A linear bias coefficient is derived by linear regression and systematically removed from tendency terms at each model integration step. This method significantly enhances both the first- (ensemble mean) and second-order moment (ensemble spread) of the ensemble forecast probability density distribution, thereby improving the probabilistic prediction skill for surface meteorological elements.

3.2.4.2   A new calculated method of radar reflectivityfor subgrid-scale precipitation

Since radar reflectivity can reveal the characteristics of strong convective weather events, forecasters pay close attention to uncertainty in model forecasts. To better quantify radar reflectivity uncertainties, we developed and implemented a new radar reflectivity calculation method for subgrid-scale precipitation in the CMA regional EPS, as presented in Eq. (13). This addresses the current limitation where simulated radar reflectivity fails to represent precipitation from the Kain-Fritsch cumulus parameterization scheme (Chen et al., 2021).

In the equation, $Z\mathrm{_{total}}$ is the new radar reflectivity, $Z\mathrm{_{micro}}$ is the radar reflectivity from the cloud microphysics parameterization scheme, ${R}_{cu}$ is the subgrid precipitation, and A and b are empirical parameters of the Z–R relationship for the radar quantitative precipitation estimation. The methodology operates by subtracting the downdraft evaporation rate from the subgrid precipitation rate, followed by layer-wise estimation of radar reflectivity using the radar-derived reflectivity–rainfall (i.e., Z–R) relationship. This new reflectivity is then integrated with microphysics-simulated radar echoes to generate a novel three-dimensional radar reflectivity field. The new reflectivity fields improved diagnostic capability for precipitation generated by the cumulus parameterization scheme, particularly under conditions dominated by subgrid convective processes. It can better simulate the radar reflectivity related to subgrid precipitation of a single ensemble member and then improve the probabilistic prediction technique of radar reflectivity in CMA regional EPS.

3.2.4.3   Extreme weather prediction product

Based on the forecast, model climate, and historical climate data of 31 ensemble members of the CMA global EPS, extreme weather prediction products are developed, including EFI of ground elements, 2-m temperature anomaly probability forecast product, and medium anomaly probability forecast product based on Kalman filter bias correction technology.

(1) The ground elements EFI product based on the climate probability of CMA global model

With the development of deterministic and ensemble forecasting for each forecasting element (precipitation, wind, temperature, and cloud cover), forecasters aim to utilize EPS for generating early warning signals of extreme events. However, it is difficult to directly compare the difference between the observed meteorological elements and the forecast output of the model. Therefore, based on the CMA global EPS, a CMA ensemble forecast method of extreme weather called EFI was developed. The principle is to calculate the difference between the cumulative model climate probability distribution and the ensemble forecast probability distribution. The calculation of the model climate percentile is a crucial step in the calculation of the EFI, but the CMA global EPS has less historical data. Thus, the model climate percentile was calculated by using a 15-day sliding time window (centered on the target date, ±7 days) and a spatial window comprising the target grid point and its 9 nearest neighbors. The daily model climate sequence was constructed from this spatiotemporal sampling, and percentiles (0.01, 0.02, ..., 0.99, 1.00) of the forecast field were extracted to form the model’s climate percentile distribution (Wang et al., 2014).

(2) Medium anomaly probability forecast productbased on Kalman filter bias correction technology

The medium anomaly probability forecast products based on Kalman filter bias correction technology are developed, including atmospheric circulation anomaly probability products of 500-hPa geopotential height (Z500), 850-hPa temperature, and 850-hPa wind. These products help to better grasp the characteristics of extreme disaster weather and improve the probability forecasting skills of extreme disaster weather.

3.3   China’s global/regional EPSs and their comparisonwith international systems

3.3.1   China’s global/regional EPSs

In response to the characteristics of initial and model errors in CMA’s GEPS and REPS, Chinese technicians have independently developed key technologies for ensemble forecast, thereby overcoming the bottlenecks in initial perturbation, model perturbation, and ensemble forecast application technologies. They have developed an independent and controllable GEPS with a horizontal resolution of 50 km for 0–15-day forecasts, along with a REPS featuring 10-km resolution and a convective-scale EPS with 3-km resolution, both for 0–3-day forecasts. These systems provide robust scientific and technological support for global ensemble prediction, propelling China’s numerical prediction operations to achieve remarkable progress and reach an advanced level. The operational implementation of the CMA’s GEPS and REPS marks a significant milestone in advancing China’s NWP capabilities. Beyond daily forecasting services and major meteorological support, the products from the CMA’s GEPS and REPS continue to supply a wide range of forecast data to the TIGGE Center. They also provide core support for the World Meteorological Centre (WMC) Beijing and the WMO Severe Weather Forecasting Demonstration Programme for Southeast Asia (SWFDP-SeA) website, further elevating China’s role in global meteorological science. Table 2 illustrates the parameter settings of the CMA’s GEPS, REPS, and global ensemble forecast (CAEF) in detail.

Table  2.  The parameter configuration for CMA’s GEPS, REPS, and CAEF

Download as CSV

	CMA-GEPS	CMA-REPS	CMA-CAEF
Model version	CMA-GFS v3.1	CMA-MESO v4.3	CMA-MESO v5.1
Horizontal resolution	0.5°	0.1°	0.03°
Vertical level	87	50	50
Domain	90°S–90°N, 180°W–180°E	10°–60°N, 70°–145°E	15°–46°N, 108°–124°E
Initial perturbation of control member	Upscaled from analysis field via CMA-GFS 4DVAR	Downscaled from CMA-GEPS	Generated from analysis field via 3DVAR
Initial perturbation scheme	SVs	ETKF	Multi-scale blending based on SVs
Model perturbation scheme	SPPT, SKEB	SPPT	SPPT
Lateral boundary perturbation scheme	/	From CMA-GEPS	Blended dynamic perturbations
Number of members	31	15	15
Forecast length	360 h (0000/1200 UTC)	84 h (0000/1200 UTC) 6 h (0600/1800 UTC)	84 h (0000/1200 UTC)
Output interval	0–84 h (3 h), 84–360 h (6 h)	1 h	1 h
Postprocessing product	GRIB2 data, ensemble prediction products for normal, extreme weather and typhoon	GRIB2 data, ensemble prediction products for normal, emergency service and typhoon	GRIB2 data, ensemble prediction products for special events service
Note: GRIB2 refers to the Gridded Binary Format Edition 2, which is a standard encoding format sponsored by the WMO for the transmission of gridded data between the national meteorological centers.

 | Show Table

DownLoad: CSV

The CMA-GEPS integrates SV initial perturbation, model random physical perturbation, and generation technologies for extreme medium-range probabilistic forecast products. It features a horizontal resolution of 0.5°, a forecast range of 0–15 days, and 31 ensemble members. Leveraging ecFlow technology (a client/server workflow package of ECMWF that enables users to run a large number of programs in a controlled environment), the operational process of the CMA-GEPS has been established, enabling daily operations at 0000 and 1200 UTC. The system provides 29 types of products, including conventional forecasts and probabilistic forecasts for extreme weather events over the 0–15-day period.

The CMA-REPS v3.0 integrates a suite of advanced techniques, including the ETKF, multi-scale blended initial perturbation, SPPT, blended dynamic perturbations of LBCs, CTVR, and a new sub-grid scale precipitation radar echo reflectivity algorithm. The system operates with 15 ensemble members and runs 4 times daily, with a forecast duration of 84 h for the 0000 and 1200 UTC cycles, and 6 h for the 0600 and 1800 UTC cycles. It primarily supports the ETKF forecast cycle, preparing data for ETKF calculations and generating forecast perturbations. The system offers 49 types of ensemble products, with a focus on probabilistic forecasts for localized severe weather and tropical cyclones.

3.3.2   Analysis of international major numericalprediction centers

3.3.2.1   Comparison on key techniques and systems ofensemble prediction

Table 3 is a comprehensive comparison of major international EPSs. The main technical characteristics of the CMA are demonstrated as follows.

Table  3.  Technical status of international major operational EPSs

Download as CSV

Country	GEPS/horizontal resolution (number of members)/technology	REPS/horizontal resolution (number of members)/technology
USA	Spectral model/25 km (31)/ETKF, SPPT, SKEB	WRF, NMMB/16 km (26, continental USA)/3 km (11, key areas)/Hybrid ensemble Kalman filter with ETR, dynamical downscaling
UK	Grid-point model MOGREPS-G/20 km (18)/EDA, SPPT, SKEB	Grid-point model MOGREPS-UK/2.2 km (3 perturbations, 18 lagged ensemble members)/Dynamical Downscaling
France	Spectral model/7.5–37 km (35)/SVs, EDA, MP	Spectral model AROME-EPS/2.5 km (16)/EDA, SPP
Germany	Grid-point model GEPS/40 km (40)/LETKF, SPP	Grid-point model COSMO-DE/2.2 km (20)/LETKF, SPP
Canada	Grid-point model/39 km (21)/LETKF, SPP, SKEB	Limited area model HREPS/15 km (21)/Dynamical downscaling, physics perturbation
Japan	Spectral model/27–40 km (27)/SVs, LETKF, SPPT	Grid-point model/5 km (21)/Multi-scale SVs
EU	Spectral model/18 km (51)/SVs, EDA, SPPT	/
China	Grid-point model GRAPES/50 km (31)/SVs, SPPT, SKEB	Grid-point model GRAPES/10 (15, China) and 3 km (15, Beijing Winter Olympics area)/ETKF, MSB, SPP
Note: MOGREPS refers to the Met Office Global and Regional Ensemble Prediction System, and “G” in MOGREPS-G stands for “global”; NMMB indicates the Nonhydrostatic Mesoscale Model on the B-grid developed by NOAA; AROME denotes the Application of Research to Operations at Mesoscale; and COSMO-DE is a short for the Consortium for Small-scale Modeling Deutschland (“Germany” in English) Edition operated by the Deutscher Wetterdienst (DWD) of Germany.

 | Show Table

DownLoad: CSV

(1) The SV initial perturbation technique, a core component of the CMA-GEPS, has been innovatively enhanced with advanced algorithms for perturbation energy modulus, parallel computing, scaling, and linear sampling of initial perturbations based on a Gaussian distribution in a three-dimensional structure. This technique effectively captures the growth of large-scale baroclinic instability perturbation energy, significantly improving the accuracy and reliability of global ensemble forecasts. Given its high technical complexity, the SV initial perturbation technique is employed by only a few leading GEPSs, including those of ECMWF, the French Meteorological Service, and the JMA. The successful implementation of this technique marks CMA’s entry into the global forefront of meteorological forecasting, with its technical level matching that of the world’s most advanced centers.

(2) The CMA-REPS incorporates several cutting-edge techniques, including the ETKF for initial perturbations, the MSB, and the CTVR. These advancements have significantly improved the accuracy of probabilistic precipitation forecasts in China. The ETKF method is based on the same principles and framework as the local ensemble transform Kalman filter (LETKF) used operationally by the Meteorological Service of Canada and the JMA. Meanwhile, the MSB and CTVR techniques are originally developed by CMA, showcasing the agency’s innovation in REPS.

(3) The stochastic perturbation methods for physical processes, such as SPPT and SKEB, utilize a first-order autoregressive stochastic process and spherical harmonic function expansion. These methods produce random functions and key parameters that are temporally and spatially correlated and normally distributed. By introducing stochastic perturbations to the tendencies of physical processes and the dissipation of small-scale kinetic energy, these techniques provide a realistic characterization of the stochastic error growth in sub-grid physical process parameterizations. This approach is consistent with the stochastic perturbation schemes employed by leading meteorological centers in Europe (e.g., ECMWF and the UK) and North America (e.g., the USA and Canada).

(4) Aiming at ensemble forecast applications, we have designed techniques for extracting extreme information, dynamically correcting systematic biases, and developing algorithms that consider radar reflectivity for model convective precipitation. Additionally, we have developed technologies for generating probabilistic forecast products such as extreme weather forecast indices. These advancements provide valuable reference information for extreme weather forecasts.

(5) An independent and proprietary operational system for ensemble forecast has been established in China, encompassing global medium-range forecasts (0–15 days), regional short-range forecasts (0–72 h), and convective-scale EPS. These advancements have significantly propelled progress in China’s numerical weather prediction capabilities. A comparison of ensemble forecast techniques and system parameters across various international numerical centers reveals that global ensemble prediction resolutions typically range from 18 to 50 km, with ensemble members varying from 18 to 51. China’s global EPS aligns closely with these international standards, featuring a horizontal resolution of 50 km and 31 ensemble members. In the realm of regional ensemble forecast, countries with smaller land areas, such as the UK, France, and Germany, have achieved resolutions of approximately 2.5 km, supported by 10–20 ensemble members. In contrast, countries with vast territories, like Canada, the U.S., and China, have generally established both regional meso- and convective-scale EPSs. China’s regional EPS has a horizontal resolution of 10 km, while its convective-scale ensemble forecasting system achieves an even finer resolution of 3 km. Such a dual-system approach highlights China’s commitments to enhancing forecast accuracy and resolution.

In summary, unique technical methods and products with independent intellectual property rights have been developed in China in the aspects of system parameters and construction schemes for GEPS and REPS.

3.3.2.2   Forecast capability comparison between CMAglobal/regional EPSs and major internationalNWP centers

Based on data provided by the JMA from the International Global Ensemble Forecast Verification Center, Fig. 1 presents a comparison of the CRPS for the 7-day ensemble forecast and the average forecast error RMSE for the 10-day ensemble forecast of the 500-hPa geopotential height field from multiple centers over the period of 2003–2022. The verification data mainly include the GEPSs of several numerical centers, such as the ECMWF, the NCEP, the United Kingdom Met Office (UKMO), the Canadian Meteorological Centre (CMC), the Bureau of Meteorology, Australia (BOM), the JMA, the Korea Meteorological Administration (KMA), and the CMA. It can be seen that since the operational implementation of the CMA-GEPS in December 2018, its forecast skill for medium-range forecasts has been comparable to that of the GEPSs of Europe and the USA.

Download

Larger image

PowerPoint slide

Fig  1.  The ensemble mean continuous ranked probability score (CRPS; top) and RMSE (bottom) of the 500-hPa geopotential height (Z500) forecasts at 168 (day 7) and 240 h (day 10) from major international global EPSs, averaged monthly from January 2003 to December 2024. Different colors represent different NWP centers. RUMS refers to the regional unified model system. NH represents the Northern Hemisphere.

4.   Challenges and prospects in thedevelopment of ensemble forecast

4.1   Scientific challenges: Spatiotemporal variability ofpredictability in convective-scale models

The predictability of convective-scale models exhibits significant spatiotemporal variability. First, the predictability of convective-scale models is influenced by the interactions of circulations across different spatiotemporal scales. Even if the magnitude of initial errors is sufficiently small, these errors can rapidly saturate and undergo upscale growth, thereby affecting the predictability of large-scale weather systems (Zhang et al., 2007). The growth of initial errors in convective-scale models has the characteristics of strong nonlinearity, with the error growth rate being approximately 10 times that of synoptic-scale models (Hohenegger and Schar, 2007). Second, convective-scale models typically employ more sophisticated and complex parameterization schemes such as cloud microphysics and turbulent diffusion. However, these approaches are subject to numerous assumptions and empirical parameters, as well as errors arising from finite-difference methods and computational truncation, all of which contribute to the strongly nonlinear growth of initial errors in convective-scale models (Mu et al., 2011; Yano et al., 2018). Third, the growth of initial errors at different scales and their impact on precipitation are closely related to environmental forcing conditions (Johnson, 2014; Zhuang et al., 2021). The relative importance of initial errors at different scales depends on the quantity and type of moist convection, with more moist convection leading to greater growth of small-scale error perturbation energy (Nielsen, 2016). Consequently, convective-scale ensemble forecast techniques are facing significant challenges.

4.2   Frontier exploration of a combination ofensemble forecast and AI

In recent years, AI has achieved groundbreaking progress in meteorological forecast applications. The computational cost of deep learning-based large-scale meteorological models is only 1 × 10⁻⁴ of that of traditional numerical models, offering a new potential pathway for the development of large-sample ensemble forecast. However, this also introduces new challenges. First, the sources of uncertainty in AI models differ from those in traditional numerical models, and the impact of these uncertainties on the forecast skill and performance of AI models requires further in-depth investigation. Second, the uncertainty characteristics of deep learning models differ from those of traditional numerical forecast, and how to appropriately describe the uncertainty features of deep learning models remains an urgent issue to be addressed. Third, the physical consistency and interpretability of AI-based ensemble forecast models remain debated, necessitating corresponding research and improvements alongside the development of AI models. Furthermore, the advancement of AI has brought new challenges and opportunities to ensemble forecast. The framework and concepts of ensemble forecast may also need to undergo transformations in the future. For instance, (1) ensemble forecast may revert to probabilistic forecast, with the core objective of obtaining the probability density function (PDF) of variables, which conceptually derives from seeking the probability distribution of forecast errors. (2) The significance of representative perturbations may decrease. As machine learning models can rapidly generate tens of thousands of ensemble forecasts, a broader spectrum of forecast possibilities can be explored, which provides better probabilistic and statistical characteristics. (3) Partial replacement of dynamical models may occur, as smaller-scale structures in the atmosphere require high-resolution dynamical models for identification. Here, dynamical models essentially act as effective but time-consuming random perturbation generators. Whether machine learning models can learn and generate such perturbations is a topic worthy of further research.

4.3   Operational challenges: Inapplicability of traditionalensemble forecast techniques

In recent years, operational ensemble forecasting development has encountered three fundamental limitations. First, global and regional operational ensemble forecast is moving toward convective-scale resolution. Under the convective-scale modeling framework, the multi-scale structures of initial condition and model errors, along with their nonlinear evolution and interactions, remain key scientific challenges in understanding the mechanisms of forecast error growth. How to construct ensemble perturbation methods is a technical bottleneck in the development of convective-scale ensemble forecast. Second, earth system models are evolving toward seamless weather–climate integrated EPSs, while traditional ensemble forecast techniques have been developed separately for temporal and spatial scales, making them inadequate for seamless spatiotemporal ensemble forecasting requirements. Third, ensemble forecast is gradually transitioning toward earth system ensemble forecast. Yet, traditional ensemble forecast techniques primarily address the uncertainty in forecasting derived from the error growth of dynamically unstable atmospheric systems. The understanding of error characteristics from coupled component models and coupled data assimilation remains insufficient, requiring the development of ensemble forecast techniques suitable for multi-component coupled models.

5.   Conclusion

In recent years, ensemble forecasting in China has achieved remarkable development and progress, from theoretical methods to operational applications, making China a world leader in operational ensemble forecasting. However, there is still a certain gap compared to the performance of EPSs at the most advanced operational centers such as ECMWF. Currently, ensemble forecasting has been intentionally developed to gradually replace deterministic forecasting, emerging as the primary form and important trend in operational forecasting. The concept of probabilistic forecasting is also gradually being popularized and accepted among researchers, forecasters, and the public. Below, we briefly summarize and discuss the challenges facing the future development of ensemble forecasting, which are also key issues that need to be addressed in the scientific research and operational development of ensemble forecasting in China.

(1) The ensemble members of global ensemble forecasting systems will advance towards convective-scale resolution in the future. This will facilitate the sampling of fine-scale features from mesoscale to convective scale and the estimation of uncertainties in ensemble forecasting on a global scale. ECMWF is at the forefront of this development, targeting a global ensemble forecast at a resolution of 5 km by 2025. Higher resolution not only poses a significant challenge to the high-performance computing resources of major operational centers but also raises highly challenging scientific and engineering questions regarding how to effectively generate multi-scale initial ensemble perturbations and model physics perturbations, while ensuring the physical consistency of these perturbations.

(2) It is important to develop integrated ensemble forecasting for global medium-range to subseasonal-seasonal (S2S) predictions. As the sources of uncertainty in subseasonal and seasonal predictions become increasingly diverse and complex, conducting ensemble forecasting for S2S is an international frontier issue. The current trend among major international operational centers is to use integrated atmosphere–ocean coupled models for ensemble forecasting of short- to medium-range and subseasonal predictions. The main frontier scientific questions in this area are how to construct initial ensemble perturbations for the atmosphere–ocean coupled system and how to account for the impact of external forcing uncertainties on S2S ensemble forecasts.

(3) With the increasing resolution of global ensemble forecasting, it is necessary to reconsider the developmental positioning and key scientific issues of regional EPSs. Regional ensemble forecasts require enhanced spatial resolutions than global ensemble forecasts, potentially advancing towards sub-kilometer resolution in the future. Its forecasting targets may also gradually encompass smaller urban-scale meteorological elements. Moreover, how to more effectively consider the coordination and nonlinear interactions of initial, lateral boundary, and model process multi-source perturbations remains a critical scientific issue in regional ensemble forecasting.

(4) The rise of AI technology presents unprecedented challenges and opportunities for ensemble forecasting. The essence of ensemble forecasting lies in statistical sampling, while AI tools excel in constructing statistical characteristics of variables and their nonlinear relationships. How to effectively apply AI technology to ensemble forecasting is a vast and cutting-edge scientific and engineering issue. Future developments in this area hold the promise of significantly advancing the level of ensemble forecasting.