HTML
--> --> -->Due to its effective probabilistic forecast information, the ensemble approach is believed necessary for the purpose of fully exploiting the potential of CAMs. Indeed, convection-allowing ensemble prediction systems (CAEPSs) play an important role in local high-resolution numerical weather prediction. Several leading forecast centers have started to run a CAEPS operationally in recent years, including COSMO-DE-EPS (Peralta et al., 2012), MOGREPS-UK (Hagelin et al., 2017), and AROME-EPS (Bouttier et al., 2012). CAEPSs have also been developed as part of NOAA’s Hazardous Weather Spring Experiment, such as SSEF (Clark et al., 2012) and NCAR’s EnKF-based EPS (Schwartz et al., 2015). See Table A1 for further details.
| Center | Model | Horizontal resolution (km) | Forecast range (h) | Members | Initial perturbations | Lateral boundary conditions | Model uncertainty |
| DWD | COSMO-DE-EPS | 2.8 | 27 | 20 | Downscaled and recentered perturbations from IFS, GFS, ICON, GSM | Downscaled from IFS, GFS,ICON, GSM | Multi-parameter |
| Met Office | MOGREPS-UK | 2.2 | 54 | 12 | Downscaled and recentered perturbations from MOGREPS-G | Downscaled from MOGREPS-G | RP |
| Météo-France | AROME-EPS | 2.5 | 45 | 12 | Downscaled and recentered perturbations from PEARP | Downscaled from PEARP | SPPT |
| RHMC | COSMO-Ru2-EPS | 2.2 | 48 | 10 | Downscaled from COSMO-S14-EPS | Downscaled from COSMO-S14-EPS | ? |
| KMA | LENS | 3 | 45 | 12 | Downscaled from EFSG | Downscaled from EFSG | RP |
| CAPS | SSEF | 3 | 60 | 20 | Downscaled from SREF | Downscaled from SREF | Multi- modelmulti- physics |
| NCAR | EnKF-based EPS | 3 | 48 | 10 | EAKF analysis | Downscaled from GFS+ random draws | ? |
TableA1. Summary of operational convection-allowing ensemble prediction systems.
However, the above CAEPSs have a common problem, known as under-dispersion. Typically, as the forecasts proceed, the ensemble spread among the ensemble members becomes too small to ensure that the observed weather pattern could be contained within the ensemble (Mascaro et al., 2010; Weyn and Durran, 2018).
Previous studies have pointed out that CAEPSs could benefit through improving the representation of model uncertainties (Johnson et al., 2014; Romine et al., 2014). To date, several approaches have been proposed in representing model uncertainties, including multi-model, multi-physics, multi-parameter, and stochastic physics. The multi-model approach considers the uncertainty of the dynamical core and parameterization by using different models (Iversen et al., 2011; Melhauser et al., 2017). Multi-physics and multi-parameter methods are used to describe the uncertainties embedded in the physical parameterizations by changing schemes (Charron et al., 2010; Hacker et al., 2011) or their empirical parameters (Gebhardt et al., 2011). Stochastic physics is more theoretical compared with the above approaches, in which the model uncertainties are described by introducing random perturbations into the physical processes, such as the well-utilized Stochastic Perturbed Parameterization Tendencies (SPPT) scheme (Buizza et al., 1999; Christensen et al., 2017), the Stochastic Total Tendency Perturbation scheme (Hou et al., 2006, 2008), and the Random Parameters/Stochastically Perturbed Parameterizations scheme (Bowler et al., 2008; McCabe et al., 2016; Ollinaho et al., 2017). In addition, considering the kinetic energy dissipation at the truncation scale, the Stochastic Kinetic Energy Backscatter scheme (Shutts, 2005; Berner et al., 2009) is another stochastic method that expresses the model structural uncertainty. Recently, a new approach based on the nonlinear forcing singular vector has been proposed to represent the combined effect of different sources of model uncertainties (Tao and Duan, 2019; Qin et al., 2020).
The above-mentioned model perturbation approaches have been widely utilized in global ensemble forecasts for a long time, and recently applied in CAEPSs for forecasting severe convective weather, tropical cyclone intensity, and so on. However, it is found that the above model perturbation approaches have limited effects in dealing with the under-dispersion problem of CAEPSs (Fresnay et al., 2012; Baker et al., 2014). Such limited effects may come from deficiencies in formulating the model uncertainties at the convection-allowing scales. Compared with global ensemble prediction systems, the representation of model uncertainty in CAEPSs is thought be lacking in systematic research and a theoretical basis, therefore making it an important issue worthy of further study.
As known, the main target of a CAEPS is to capture the probability of occurrence of rapidly evolving small- to mesoscale strong events. From the ensemble forecast viewpoint, the rapidly nonlinear growing perturbations should be described in the CAEPS’s formulations both for initial errors and model uncertainties. Mu et al. (2003) proposed the concept of conditional nonlinear optimal perturbation (CNOP). The CNOP represents the maximum nonlinear evolution of initial perturbations satisfying constraint conditions at the prediction time, which is thought to be greatly important for the predictability (Wang and Tan, 2010). Since Mu et al. (2003), CNOP has been applied in studying the various predictability issues, such as ENSO predictability, target observations, and so on. Duan and Huo (2016) also proposed the approach of orthogonal CNOPs, which extended the CNOP to ensemble forecasting. Furthermore, considering the predictability caused by the model error, Mu et al. (2010) extended the CNOP approach to search for the optimal model parameter perturbations, which resulted in the largest departure from a given reference state at a prediction time with physical constraints. This approach is referred to as CNOP-P, which has been extensively employed in understanding the climate model sensitivity to model parameters and investigating the impact of parameter errors on ENSO predictability, amongst other problems (Duan and Zhang, 2010; Yu et al., 2012; Yin et al., 2015).
As discussed above, a good CAEPS should reasonably describe the rapidly growing nonlinear errors or maximum growing perturbations at a short prediction time either due to initial errors or model uncertainties. In this sense, the CNOP-P approach may have potential to detect and describe the maximum perturbation growth caused by the model parameter uncertainties for CAEPSs. Compared with current formulations of model uncertainties in CAEPSs, the CNOP-P approach is expected to be more objective. However, the original computation method of CNOP-P requires integration of an adjoint model to calculate the cost function gradient. This limits the practical use of CNOP-P in operational applications. Wang and Tan (2010) proposed a pragmatic method to determine the CNOP by using an ensemble approach, which is easier to implement and more time-saving than the adjoint-based method. Thus, our research seeks to exploit a novel method to formulate the model uncertainties for CAEPSs by using the CNOP-P concept and the approach to determining the CNOP proposed by Wang and Tan (2010).
The rest of the paper is organized as follows. Section 2 briefly describes the CNOP-P approach utilized in this paper and the experimental design. The results of sensitivity tests and the impact of the CNOP-P approach on the ensemble forecasts are reported in section 3. Finally, conclusions and future work are summarized in section 4.
2.1. CNOP-P and its computation
According to Mu et al. (2010), in a nonlinear prediction model M,
The optimal parameter perturbations p' should correspond to maximum the prediction increments y' with a constraint of || p'||≤σ (σ is a small positive number, given here as 1). The p' can be obtained through solving the following equation:
In order to reduce the dependency of J(p') on initial values, Yin et al. (2015) introduced the ensemble-averaged type of cost function, which is constructed by using N short-term time series predictions with independent initial values U0,k (k=1,2,…,N) as follows:
The above cost function can be rewritten into a minimization problem for convenience of solution:
The key to this minimization problem is the calculation of its gradient with respect to the parameter perturbation, which can be expressed as:
As derived by Wang and Tan (2010) and Yin et al. (2015), the solution of Eq. (5) can be obtained by the following ensemble approach.
Firstly, for each initial value U0,k, n parameter perturbation samples p'i (i=1,2,…,n) are used to generate n corresponding prediction increment samples y'i,k (i=1,2,…,n; k=1,2,…,N). The parameter perturbations and prediction increment matrices are denoted as P' and Y'k, respectively:
The second step is to obtain a series of orthogonal basis functions of P' by using singular value decomposition:
where, Ap is the left singular vector of P', which is taken as the orthogonal basis function.
The third step is to estimate the gradient of the cost function based on ensemble projection. With Ap as the projection matrix, each parameter perturbation can be projected to the ensemble space:
where α is the parameter perturbation after projection.
According to the tangent linear relationship between p' and y'k, each prediction increment can also be projected to the same ensemble space:
where M' is the tangent linear model of M.
Then, a statistical model between p' and y'k can be established:
According to the orthogonality of Ap, we can obtain
The last step is to compute CNOP-P and detect the most sensitive parameters. The Spectral Projected Gradient 2 (SPG2) optimization algorithm of (Birgin et al., 2001) is used to obtain a result iteratively. To guarantee nonlinearity of the optimal solution, GRAPES-Meso (Global/Regional Assimilation and Prediction System, mesoscale version), with a horizontal grid size of 3 km, is used to calculate the cost function value. Meanwhile, the gradient of the cost function and constraint condition of parameter perturbations are required during iteration. The constraint condition is expressed as:
There are m parameters involved in the calculation of CNOP-P. Perturbations are standardized by dividing their corresponding default values. In this way, the greater the perturbation weight of a parameter in the CNOP-P, the greater the relative sensitivity of this parameter.
2
2.2. Experimental design
32.2.1. Model configuration
The model utilized here is the operational GRAPES-Meso, version 4.1, with a convection-allowing resolution setup. This model solves the fully compressible nonhydrostatic equations using a semi-implicit, semi-Lagrangian integration scheme on a spherical Arakawa-C grid horizontally, and Charney–Phillips staggered grids vertically. It has 50 vertical levels up to 10 hPa. For the 3 km resolution, the 30-s time step is adopted.Physical parameterization schemes used here include the Rapid Radiation Transfer Model longwave radiation scheme (Mlawer et al., 1997), the Dudhia shortwave radiation scheme (Dudhia, 1989), the Monin–Obukhov surface layer scheme, the Noah land surface scheme (Chen and Dudhia, 2001), the MRF boundary layer scheme (Hong and Pan, 1996), and the WSM6 cloud microphysics scheme (Hong and Lim, 2006). The cumulus convection scheme is not applied.
The initial conditions are generated by the GRAPES Cloud Analysis System (GCAS) (Qu et al., 2010) with a 3-km resolution, which has been tested in GRAPES-Meso_3km from 2015 to the present (Zhu et al., 2017). By using GCAS, less model spin-up time in the precipitation forecast was observed. NCEP Final Operational Global Analysis data at a spatial resolution of 1° are used as the background. Based on this background, GCAS integrates data sources from the Doppler weather radar 3D mosaic reflectivity data at a spatial resolution of 0.01°, blackbody temperature and total cloud amount data from the FY-2G geostationary satellite. The lateral boundary conditions are also provided by the NCEP Final Operational Global Analysis with a 1° horizontal resolution.
Because CAEPSs are mainly designed for obtaining the forecast probability of convection, a region with frequent occurrence of locally developed convection is preferred as the modeling area. As is known, most of the precipitation during the pre-summer rainy season in South China is of a convective nature with mesoscale organizational characteristics. Therefore, the model domain covers South China (18.5°–27.5°N, 105.0°–118.5°E) in this study. In order to reduce the influences from the boundaries on the CNOP-P computation, a smaller area (20°–26°N, 106.5°–117°E) is selected for calculating the CNOP-P, as shown in Fig. 1.
Figure1. Model forecast domain (D1: 18.5°–27.5°N, 105.0°–118.5°E) and computation domain for CNOP-P (D2: 20°–26°N, 106.5°–117°E).3
2.2.2. Experimental design for detecting the sensitive physical parameters based on CNOP-P
For the purpose of better isolating the physical processes influencing the initiation and development of convection, five locally developed rainstorm cases in 2017 that occurred in South China are chosen for constructing the ensemble to compute the CNOP-P, i.e., N=5 in Eq. (3). These five cases occurred under weak synoptic forcing. Here, the 12-h cumulative precipitation is used to calculate the cost function for detecting the most sensitive parameters, which means that CNOP-P could lead to the maximum variations in precipitation during the first 12 h of model integration. Details of the five cases are given as Cases 1–5 in Table 1. Then, a supplementary experiment is carried out to confirm the robustness of the above parameter detection. Five heavy rainfall cases under strong synoptic forcing are selected to compute the CNOP-P, which are given as Cases 6–10 in Table 1. It is worth mentioning that all the forecasts in Table 1 are initialized about 6 h before the onset of rainfall events. These cases are classified as well as documented by the forecasters in the Central Meteorological Bureau (Zhang et al., 2018).| Case | Initial forecast time (UTC) | Period of heavy rainfall (UTC) | Weather type | Trigger condition |
| 1 | 1200 6 May | 2100 6 May–0000 7 May | Rainstorm | Convective systems in the warm sector |
| 2 | 0000 9 May | 0600 9 May–0900 9 May | Rainstorm | Convective systems in the warm sector |
| 3 | 0000 15 Jun | 0900 15 Jun–1200 15 Jun | Wide-area rainstorm | Convective systems in the warm sector |
| 4 | 0000 20 Jun | 0900 20 Jun–1200 20 Jun | Wide-area rainstorm | Convective systems in the warm sector |
| 5 | 0000 26 Jun | 0900 26 Jun–1200 26 Jun | Rainstorm, thunderstorm gale | Convective systems in the warm sector |
| 6 | 0000 8 May | 0900 8 May–1200 8 May | Wide-area rainstorm, thunderstorm gale | Cold front |
| 7 | 1200 14 May | 1800 14 May–2100 14 May | Wide-area rainstorm | Cold front |
| 8 | 0000 23 May | 0900 23 May–1200 23 May | Wide-Wide-area rainstorm, thunderstorm gale gale | Cold front |
| 9 | 0000 6 Jun | 0600 6 Jun–0900 6 Jun | Wide-area rainstorm, thunderstorm gale | Cold front |
| 10 | 0000 16 Jun | 0700 16 Jun–1000 16 Jun | Wide-Wide-area rainstorm, thunderstorm gale gale | Cold front |
| 11 | 1200 20 Apr | 2100 20 Apr–0000 21Apr | Rainstorm, thunderstorm gale, hail | Cold front |
| 12 | 1200 3 May | 1800 3 May–0000 4 May | Wide-area rainstorm, thunderstorm gale, hail | Squall line |
| 13 | 1200 23 May | 1700 23 May–2000 23 May | Wide-area rainstorm | Cold front |
| 14 | 1200 1 Jul | 0000 2 Jul–0300 2 Jul | Rainstorm, thunderstorm gale | Convective systems in the warm sector |
| 15 | 0000 7 Jul | 0700 7 Jul–1000 7 Jul | Rainstorm | Convective systems in the warm sector |
| 16 | 0000 10 Jul | 0600 10 Jul–0900 10 Jul | Rainstorm, thunderstorm gale | Convective systems in the warm sector |
| 17 | 0000 19 Jul | 0600 19 Jul–0900 19 Jul | Rainstorm, thunderstorm gale | Convective systems on the edge of subtropical high and in the warm sector |
| 18 | 0000 21 Jul | 0800 21 Jul–1000 21 Jul | Rainstorm, thunderstorm gale | Convective systems on the edge of subtropical high and in the warm sectorarm sector |
| 19 | 0000 31 Jul | 0600 31 Jul–0900 31 Jul | Wide-area rainstorm, thunderstorm gale, hail | Cold front |
| 20 | 0000 1 Aug | 0600 1 Aug–0900 1 Aug | Rainstorm, thunderstorm gale | Convective systems in the warm sector |
| 21 | 0000 10 Aug | 0800 10 Aug–1000 10 Aug | Rainstorm, thunderstorm gale | Convective systems on the edge of subtropical high |
| 22 | 0000 22 Aug | 1200 22 Aug–1500 22 Aug | Rainstorm, thunderstorm gale | Typhoon |
| 23 | 0000 23 Aug | 0800 23 Aug–1000 23 Aug | Rainstorm | Typhoon |
Table1. Twenty-three typical heavy rainfall cases that occurred in South China in 2017.
Fifteen parameters as shown in Table 2 are selected from the boundary layer and microphysics, which are thought to be the most important two processes in influencing convective events. Theoretically, all the physical parameters should be selected to perturb for calculating the CNOP-P. However, due to the high computing costs, we only selected 15 parameters in this study. According to previous research on CAEPSs (Baker et al., 2014; Gebhardt et al., 2011; McCabe et al., 2016), it was believed that this selection would not lose the generality. When we generate the parameter perturbation samples, only one parameter is changed each time, and the perturbed value is ±50% of its default value.
| Index | Scheme | Parameter | Description | Default |
| P1 | MRF | rlam | Asymptotic mixing length (m) | 150 |
| P2 | MRF | brcr | Critical Richardson number | 0.5 |
| P3 | MRF | pfac | Profile shape exponent for calculating the momentum diffusion coefficient | 2 |
| P4 | MRF | cfac | Coefficient for Prandtl number at the top of the surface layer | 7.8 |
| P5 | WSM6 | avtr | A constant for terminal velocity of rain | 841.9 |
| P6 | WSM6 | avtg | A constant for terminal velocity of graupel | 330 |
| P7 | WSM6 | avts | A constant for terminal velocity of snow | 11.72 |
| P8 | WSM6 | n0r | Intercept parameter of rain (m?4) | 8.0 × 106 |
| P9 | WSM6 | n0g | Intercept parameter of graupel (m?4) | 4.0 × 106 |
| P10 | WSM6 | deng | Density of graupel (kgm?3) | 500 |
| P11 | WSM6 | peaut | Mean collection efficiency | 0.55 |
| P12 | WSM6 | xncr | Number concentration of cloud water droplet (m?3) | 3.0 × 108 |
| P13 | WSM6 | dimax | Maximum value for cloud ice diameter (m) | 5.0 × 10?4 |
| P14 | WSM6 | r0 | Critical mean droplet radius at which auto-conversion begins (m) | 8.0 × 10?6 |
| P15 | WSM6 | qs0 | Threshold amount for aggregation to occur (kg kg?1) | 6.0 × 10?4 |
Table2. The chosen physical parameters for computing CNOP-P.
3
2.2.3. Ensemble experiment design
By applying the CNOP-P, a group of the most sensitive parameters can be detected from the 15 parameters, which correspond to the maximum forecast error growth caused by the parameter uncertainties. Based on these most sensitive parameters, a novel formulation of model uncertainty is designed, and the ensemble experiment design is given below.As is known, currently existing systems typically use ensemble sizes of around 10 or more members, like MOGREPS-UK, AROME-EPS, NCAR’s EnKF-based EPS, etc. Due to the available computer resources, a 10-member ensemble run was designed in this study, which included one control member and nine perturbed members. For the purpose of investigating the advantages of this new approach in describing model uncertainty, we established the ensemble prediction system by only taking into account the model uncertainty. The well-utilized SPPT scheme is adopted as the benchmark, which represents the typical model uncertainty formulation in current operational global ensemble and regional ensemble systems (Bouttier et al., 2012; Leutbecher et al., 2017). Two sets of ensemble forecasts are carried out with different model perturbation methods. The first set of ensemble forecasts uses SPPT scheme, which is called EXP1 for short. The second uses the CNOP-P-based approach, which is called EXP2 for short. Next, the details of the two model uncertainty formulations are described.
In the SPPT scheme, a stochastic perturbation is added to the net tendency term of physical parameterizations:
where X0 is the net tendency term, Xp is the perturbed net tendency term, and φ is a 3D random field. The scheme firstly introduces a random field generator with certain space structures and time-correlated characteristics. This generator is an expansion of spherical harmonics along the horizontal directions based on first-order autoregressive random processes. Then, the random number is multiplied into the net tendency of all the physical processes. Aimed at maintaining the numerical stability, vertical profiles are introduced to reduce the perturbations near the surface and the model top. The following parameter configurations are based on Yuan et al. (2016): the random field follows a Gaussian distribution with a mean value of 0, a standard deviation of 0.27, a perturbation range of [0.2, 1.8], and temporal decorrelation scale of 6 h. Because of time limitations, this paper does not consider the tuning of the SPPT scheme, or its reformulation, e.g., by using the iSPPT formulation (Christensen et al., 2017). It should be noted that the SPPT scheme used here has been well tuned for the GRAPES mesoscale ensemble forecast system with a 15-km resolution, and has been operationally implemented since 2014 (Yuan et al., 2016).
For the CNOP-P-based approach, based on the detected sensitive parameters as described above, a random parameter algorithm (McCabe et al., 2016) is introduced to construct the CNOP-P-based formulation of the model uncertainty.
By giving a physically reasonable range [Pmin, Pmax] for varying to each parameter, the detected sensitive parameters are first mapped onto the range [?1, 1]:
where η is the mapped value, and Pmin and Pmax are the minimum and maximum values.
By using the first-order auto-regression model, the mapped parameters are stochastically perturbed and updated at time step t. The parameters are perturbed stochastically and independently at regular intervals throughout the relations:
where ηt is the value of mapped parameters at time t, r is the autocorrelation coefficient, and εt is the stochastic shock term. εt is formulated in terms of r and stochastic number At. r is expressed by the update time interval Δt and the characteristic time scale τp. According to McCabe et al. (2016), the update time interval is set to 5 min and the characteristic time scale is set to 7 d. Then, ηt will be mapped back to original range using Eq. (14).
3
2.2.4. Verification and metrics
For the purpose of systematically comparing the forecast skills using the two different approaches of model uncertainty formulation, 13 heavy rainfall cases that occurred during the flood season of South China in 2017 were selected to carry out reforecast experiments. These cases are given as Cases 11–23 in Table 1, including heavy rainfall events with strong and weak synoptic forcing.The near-surface variables, including 2-m specific humidity (Q2m), 10-m wind speed (V10m) and 2-m temperature (T2m), as well as hourly precipitation (R), are mainly evaluated in this study. Hourly observations from automatic weather stations are used in the verification.
The verification metrics include the ratio of spread to root-mean-square error (RMSE), outlier frequency, and continuous ranked probability score (CRPS). In an ideal ensemble, the ratio of spread to RMSE (hereafter, S/R) should be as close to 1 as possible; that is, the spread can represent the error evolution characteristics of the system. In reality, spread is often less than the RMSE owing to the observation error. Here, observation error is not taken into account in the verification in our study:
where (ˉ) denotes the mean value of all grid points, fi is the forecast of ensemble member i, n is the number of ensemble members, fmean is the ensemble mean, f0 is the control forecast, and o is the observation value.
Outlier frequency is used to evaluate the reliability of the CAEPS. Suppose there are N valid samples and n ensemble members. After sorting the forecast of ensemble members from small to large, n+1 intervals are generated. Si denotes the occurrence frequency of observation in the ith interval. The Talagrand distribution map can be drawn from the probability distribution:
In an ideal Talagrand distribution, there is not much difference between the probabilities of each interval, but in reality, the two ends are larger. Outlier frequency is a measure of the value at both ends of the distribution, which denotes the average probability that the observations fall outside the range of the forecasts. The ideal value is 1/(n+1).
CRPS is a probabilistic skill score to measure the differences in cumulative distribution functions between observations and forecasts. It can be expressed as the integral of the Brier score for all possible thresholds of the meteorological variables:
where F(x) and Fo(x) represent the cumulative distribution functions of the ensemble forecast and observation, respectively. The value for an ideal forecast of CRPS is 0. A lower CRPS indicates a more skillful forecast.
3.1. Sensitive parameter detection
As described in section 2.1, iteration is necessary in solving Eq. (5) for detecting the sensitive parameters. In this study, the iteration is stopped when the cost function becomes almost unchanged (< 1%), and the spectral projection gradient decreases by five to six orders of magnitude. Table 3 shows the results of the iterative process. It is found that the above criteria start being met from step 5 of iteration, and then the parameters remain stably ranked. Here, the top eight parameters are regarded as the most sensitive ones, denoted as P1, P13, P7, P5, P3, P10, P4 and P2 in Table 2. We attempt to provide a physical explanation of these eight sensitive parameters.| Iteration | Spectral projection gradient | Cost function | Sensitivity ranking corresponding to CNOP-P |
| 0 | 1.192 | ?923831.634 | 10 13 3 11 8 5 12 4 9 6 2 15 7 14 1 |
| 1 | 7.614 × 10?1 | ?7120544.451 | 1 5 3 2 14 6 11 10 12 13 8 9 15 7 4 |
| 2 | 2.052 × 10?2 | ?9873203.670 | 1 5 13 6 2 14 10 12 11 8 3 7 4 9 15 |
| 3 | 6.479 × 10?4 | ?9735736.603 | 1 13 5 7 10 2 12 4 3 11 9 15 6 14 8 |
| 4 | 4.630 × 10?5 | ?9745800.479 | 1 13 5 7 10 3 2 4 9 12 15 11 6 8 14 |
| 5 | 4.046 × 10?4 | ?9856915.210 | 1 13 5 7 10 3 2 4 9 12 15 11 6 8 14 |
| 6 | 2.211 × 10?4 | ?9615821.473 | 1 13 7 5 10 3 4 2 12 9 15 11 8 14 6 |
| 7 | 2.256 × 10?5 | ?9657034.976 | 1 13 7 5 10 3 2 4 12 9 15 11 8 6 14 |
| 8 | 2.773 × 10?5 | ?9727943.338 | 1 13 7 5 10 3 2 4 12 9 15 11 8 6 14 |
| 9 | 3.978 × 10?5 | ?9758126.346 | 1 13 7 5 3 10 4 2 12 9 15 11 8 6 14 |
| 10 | 1.212 × 10?5 | ?9686227.187 | 1 13 7 5 3 10 4 2 12 9 15 11 8 14 6 |
| 11 | 1.034 × 10?5 | ?9658981.241 | 1 13 7 5 3 10 4 2 12 9 15 11 8 14 6 |
Table3. The iteration of CNOP-P computation.
In the MRF scheme, the PBL and free atmosphere vertical diffusion is calculated based on the local K approach. The asymptotic mixing length (P1) is a physical quantity describing the size of the large energy-containing eddies in turbulence mixing, and is thus crucial in determining the vertical diffusion intensity. The mixing length scale is detected as the number one sensitive parameter, illustrating the rather important role of boundary layer processes in affecting locally developed convections. P2, the critical Richardson number, is a key parameter for calculating boundary layer height. P3 is a profile shape exponent for calculating the momentum diffusion coefficient, which determines the magnitude of the vertical diffusion as well as the height at which the diffusion reaches a maximum. P4, the Prandtl number coefficient, is involved in computing eddy diffusivity for temperature and moisture within the boundary layer. In the WSM6 microphysics scheme, the constants for the terminal velocity of rain and snow (P5, P7) play important roles in the vertical distribution of hydrometeors and auto-conversion in clouds. P10, the density of graupel, varies greatly, and can lead to obviously different characteristics of cloud microphysics. Setting P10 as a constant may underestimate or overestimate the precipitation efficiency of clouds. The maximum value of the cloud ice diameter, named P13 here, is related to the maximum mass of ice crystals, which is critical in auto-conversion from ice to snow, and further influences the rainfall intensity for the strong convection cases.
These eight sensitive parameters are applied to address the model uncertainty using Eqs. (14)–(17). Table 4 gives the range of perturbations for these parameters. It should be noted that this range is a little different from the perturbed range in the calculation of CNOP. To ensure the representativeness and independence of samples, we perturbed ±50% of the defaults to generate the parameter perturbation samples in computing CNOP-P. Here, however, the perturbations need to be closer to the actual values in the whole variation ranges, which should be large enough to add reliable variability into the forecasts, and small enough to make the forecasts reasonable. Although these exact thresholds are not clear, we define the values by referring to previous research (McCabe et al., 2016; Ollinaho et al., 2017).
| Index | Parameter | Min | Default | Max |
| P1 | rlam | 30 | 150 | 450 |
| P2 | brcr | 0.125 | 0.5 | 1 |
| P3 | pfac | 1 | 2 | 3 |
| P4 | cfac | 3.9 | 7.8 | 15.6 |
| P5 | avtr | 420 | 841.9 | 1263 |
| P7 | avts | 5 | 11.72 | 18 |
| P10 | deng | 200 | 500 | 900 |
| P13 | dimax | 2 × 10?4 | 5 × 10?4 | 8 × 10?4 |
Table4. Sensitive parameters and their ranges.
The following are the results of the supplementary experiment, which, using five heavy rainfall cases under strong synoptic forcing, was in sharp contrast to the above in which we selected five cases with weak synoptic forcing. The sensitive parameters detected by this experiment are P1, P13, P5, P7, P2, P3, P14 and P12, as shown in Table 5. Only two parameters changed compared with the experiments shown in Table 3. The result indicates that the detection and sensitivity ranking of parameters is relatively stable and robust.
| Iteration | Spectral projection gradient | Cost function | Sensitivity ranking corresponding to CNOP-P |
| 0 | 1.216 | ?1923239.880 | 10 13 3 11 8 5 12 4 9 6 2 15 7 14 1 |
| 1 | 3.401 × 10?1 | ?12148995.40 | 1 3 2 5 13 7 4 9 11 14 6 8 15 12 10 |
| 2 | 8.955 × 10?3 | ?12753984.43 | 1 5 3 14 12 13 2 4 7 9 11 8 10 6 15 |
| 3 | 3.236 × 10?4 | ?12522691.23 | 1 5 13 7 2 3 14 12 11 9 4 15 10 6 8 |
| 4 | 1.971 × 10?4 | ?12598208.52 | 1 13 5 7 2 3 14 12 11 4 9 10 15 6 8 |
| 5 | 1.201 × 10?4 | ?12565818.29 | 1 13 7 5 2 3 14 12 11 4 9 10 15 8 6 |
| 6 | 8.993 × 10?5 | ?12602431.25 | 1 13 7 5 2 3 14 12 11 9 4 10 15 6 8 |
| 7 | 6.820 × 10?5 | ?12546727.08 | 1 13 7 5 2 3 14 12 11 9 4 10 15 6 8 |
| 8 | 2.925 × 10?5 | ?12539894.94 | 1 13 7 5 2 3 14 12 11 9 4 10 15 6 8 |
| 9 | 1.211 × 10?4 | ?12602381.10 | 1 13 7 5 2 3 14 12 11 9 4 10 15 6 8 |
| 10 | 4.626 × 10?5 | ?12584054.48 | 1 13 7 5 2 3 14 12 11 9 4 10 15 6 8 |
| 11 | 2.772 × 10?5 | ?12547402.31 | 1 13 5 7 2 3 14 12 11 9 4 10 15 6 8 |
Table5. The iteration of CNOP-P computation using five cases under strong synoptic forcing as initial conditions.
2
3.2. Evaluation of the ensemble forecasts
The purpose of this paper is to introduce the CNOP-P-based model uncertainty formulation and understand its advantages for CAEPSs. For a relatively fair comparison, we chose a total of 13 heavy rainfall cases to understand the systematic effects, and calculated the averaged differences of metrics over a forecasting time of 6–24 h between EXP1 and EXP2, i.e., EXP2 minus EXP1, to investigate the added value of the CNOP-P-based approach. According to the definition of these metrics, EXP2 performs better when the difference of S/R is greater than 0, or the differences of outlier frequency and CRPS are less than 0, and vice versa. Figure 2 shows the 6–24-h averaged differences of S/R on four variables for the 13 cases, with T2m, V10m, Q2m and R shown in Figs. 2a–d, respectively. From all the positive values of S/R differences in Fig. 2, it is obvious that the CNOP-P-based method, i.e., EXP2, exhibits systematic improvement with respect to the S/R relationship for all the near-surface variables and precipitation. This result indicates that the CNOP-P-based method has great potential in improving the under-dispersive problem of current CAEPSs.
Figure2. The 6–24-h time-averaged differences of S/R (spread/RMSE) for 13 cases. Black lines represent EXP2’s S/R on four variables, and columns represents the difference between EXP2 and EXP1 (i.e. EXP2 minus EXP1). (a) 2-m temperature; (b) 10-m wind speed; (c) 2-m specific humidity; (d) hourly precipitation.
Figure3. As in Fig. 2 but for outlier frequency.Similar to Fig. 2, Fig. 3 gives the 6–24-h averaged differences of outlier frequency for each variable. As mentioned above, outlier frequency denotes the average probability that the observations fall outside the range of the forecasts. As shown in Fig. 3, the time-averaged differences of outlier frequency for T2m, V10m and Q2m exhibit negative values for all 13 cases. Only two exceptions are found in the precipitation outlier frequency difference (Fig. 3d). It can be concluded that, on average, the CNOP-P-based method can reduce the probability that the observations fall outside the range of the forecasts, thus enhancing the reliability of the CAEPS.
The CRPS is frequently used to assess the respective accuracy of two probabilistic forecasting models. The 6–24-h time-averaged differences of CRPS for the two probabilistic forecasts are shown in Fig. 4 for each variable. According to the definition of CRPS, a negative time-averaged difference means a better probabilistic forecast. From Fig. 4, it can be found that, for all four variables, about 11 out of 13 cases show negative differences, indicating a better forecast corresponding to the cumulative distribution function of the observations.
Figure4. As in Fig. 2 but for CRPS.As shown in Figs. 2–4, the above verification results indicate that the CNOP-P-based method can bring more spread, and has relatively reliable and skillful probabilistic forecasts compared with the SPPT scheme. That is, perturbing the most sensitive parameters or physics may have more potential to improve CAEPSs than perturbing the net physics tendency.
Next, we check the different effects of the two approaches in the vertical direction. Figure 5 shows the vertical profiles of ensemble spread for the 24-h forecast of zonal wind, temperature and specific humidity, the profiles of which are the averages of the 13 cases. For the zonal wind field in Fig. 5a, EXP2 can obtain larger spread within the boundary layer, while EXP1 seems to be more dispersive in the middle atmosphere. This may partly reflect the characteristics of the two approaches in formulating the model uncertainties, but on the whole the difference is not significant. For temperature and specific humidity in Fig. 5b and Fig. 5c, EXP2 has consistently and statistically significant larger spread over the troposphere compared with EXP1. This again confirms the merits of the CNOP-P-based approach, because through detecting the most sensitive parameters from CNOP-P, the process-specific and event-oriented model uncertainties can be realized. More specifically, focusing on the CAEPS of convective events in this study, the most sensitive parameters embedded in the PBL and microphysics are detected by using CNOP-P. Thus, the formulation of model uncertainty of these two key processes for convection, which has solid theoretical basis, can bring more spread for temperature and humidity over the troposphere.
Figure5. Vertical profiles of ensemble spread for the 24-h forecasts of 13 cases: (a) zonal wind (units: m s?1); (b) temperature (units: °C); (c) specific humidity (units: g kg?1).Finally, it should be noted that some ensembles might not necessarily be under-dispersive in terms of all parameters, so that our study may be more about creating “a right kind of spread” rather than increasing the spread in absolute terms. A sensitivity test was conducted to explain why we chose to perturb the first eight sensitive parameters in the ranking of Table 3. Five comparison experiments were implemented, each of which perturbed the first 4, 6, 8, 10 and 12 parameters in the sensitivity ranking of Table 3. We chose Case 1 in Table 1 for ensemble prediction, and calculated the differences of three metrics over the 6–24-h forecasting time (Figs. 6–8). From the results of the three verification metrics, we found that the perturbation scheme with eight parameters is more effective than that with four or six parameters. However, increasing to more than eight parameters does not obviously improve the forecasts. It can be seen that the uncertainty of parameters is mainly due to the important sensitive parameters.
Figure6. The 6–24-h time-averaged S/R (spread/RMSE) for the case initialized at 1200 UTC 6 May. EXP_4 (red), EXP_6 (blue), EXP_8 (black) EXP_10 (grey), and EXP_12 (green) are five ensemble forecasts that perturbed the first 4, 6, 8, 10 and 12 parameters, respectively: (a) 2-m temperature; (b) 10-m wind speed; (c) 2-m specific humidity; (d) hourly precipitation.
Figure7. As in Fig. 6 but for outlier frequency.
Figure8. As in Fig. 6 but for CRPS.The CNOP-P approach is very robust in detecting the sensitive parameters when applied to the convection-allowing-scale model.
The CNOP-P-based method has potential in improving the under-dispersive problem of current CAEPSs, and has more reliable forecast skill for T2m, V10m, Q2m and precipitation compared with the SPPT scheme.
The CNOP-P-based method brings more spread for humidity and temperature over the troposphere owing to its process-specific and event-oriented formulation of model uncertainties.
Although the cases we selected in this study are specific to South China and may therefore limit the conclusions, the experiments with these 13 heavy rainfall cases in the 2017 flood season strongly suggest a very bright future in terms of the application of the CNOP-P-based approach in CAEPSs.
There is, however, still a practical and critical issue to resolve in terms of how to select the parameters if two of the eight parameters are different when different cases are used. The CNOP-P approach has been proven theoretically in previous studies and in some climate model applications, and it is believed in practical NWP applications that a stable group comprising the most sensitive parameters could be obtained if the cases are chosen to share the similar physical features or processes. The supplementary experiment in our study illustrates this issue, i.e., the parameters based on the “strong forcing” cases were similar to those based on the “weak forcing” cases. Much research is needed in which the CNOP-P approach is further applied in operational CAEPSs—such as, its application in heavy rainfall events in the mid-to-high latitudes, the impact of parameter perturbation constraints, the use of spatiotemporal random perturbation approach when perturbing the parameters, and the choice of cost functions.
Acknowledgements. We sincerely appreciate the constructive comments and suggestions of the anonymous reviewers. This work was supported by the National Key Research and Development Program of China (Grant No. 2017YFC150 1904).
APPENDIX
