Comparison of Nonlinear Local Lyapunov Vectors with Bred Vectors, Random Perturbations and Ensemble Transform Kalman Filter Strategies in a Barotropic Model

In the past two decades, ensemble forecasting has been substantially developed to become a powerful approach that improves numerical weather prediction (NWP). The basic principle of the generation of initial ensemble members is to sample the uncertainties related to the initial analysis (Epstein, 1969; Leith, 1974). Various ensemble generation schemes based on dynamical error growth theory have been tested and used in weather prediction centers; for example, the bred vector (BV) method (Toth and Kalnay, 1993, 1997) used at NCEP, and the singular vector (SV) method (Lorenz, 1965; Molteni et al., 1996) at ECMWF. Afterwards, data assimilation (DA) schemes were further combined with the dynamical methods to better sample the analysis uncertainties, such as in the ensemble transform Kalman filter (ETKF) scheme (Bishop et al., 2001; Wang and Bishop, 2003; Wei et al., 2006, 2008) which is an improved version of the BV method, and in the introduction of advanced DA schemes in the SV method (Lawrence et al., 2009).

Numerous studies have been carried out to compare the BV and ETKF approaches (Wang and Bishop, 2003; Bowler, 2006; Wei et al., 2006, 2008; Descamps and Talagrand, 2007), from which is has been found that the better forecast skill of the ETKF scheme than the breeding scheme is mainly associated with two specific features. One is that the ensemble perturbations of the ETKF scheme can better reflect the geographical variations of analysis error variance than those of the breeding scheme. This may be because the BV perturbations are tuned by a time-invariant estimate of analysis error variance ("mask"; Toth and Kalnay, 1997), while the ETKF perturbations at each time are updated through the DA algorithm. The other is that the ETKF perturbations are more independent and uncorrelated than the BVs. The BVs are naturally evolved perturbations that represent the unstable modes associated with the day-to-day dynamical flow. However, dynamical flow would drive the BVs to project on the growing type of directions and thus increase the similarity of BVs, especially in local regions. Such similarities are reduced by the nonlinear error growth in NWP models, but can also be seen by the eigenvalue spectrum for BVs (Wang and Bishop, 2003; Bowler, 2006). In contrast, the ETKF ensemble perturbations are orthogonal in observation space (Wang and Bishop, 2003; Wei et al., 2006).

In view of the interdependence among BVs, (Feng et al., 2014) developed the nonlinear local Lyapunov vectors (NLLVs), which are the nonlinear extension of the linear Lyapunov vectors (LVs). The traditional LVs define a set of time-dependent orthogonal perturbation structures (Kalnay, 2002; Trevisan and Palatella, 2011). These perturbations point to the directions with different growth rates in phase space specified by positive, null and negative Lyapunov exponents (Benettin et al., 1980; Wolf et al., 1985; Fraedrich, 1987; Trevisan and Legnani, 1995; Vannitsem and Nicolis, 1997). Therefore, the LVs can be used to identify the unstable, neutral and stable subspaces. In reality, on the one hand, the unstable perturbation subspace has a significantly smaller dimension than the stable subspace, and thus is much easier to be sampled (Toth and Kalnay, 1997); while on the other hand, any perturbation, through model integration, will eliminate the non-growing components in it, and project onto the fast-growing subspace with relatively small dimension. Therefore, the identification of unstable perturbations is critical to the improvement of NWP (Zhang et al., 2015), like the accurate estimation of error variances in the background forecasts of DA (Trevisan and Uboldi, 2004; Carrassi et al., 2007), and the generation of ensemble perturbations in ensemble forecasting (Toth and Kalnay, 1997), etc. This is also why the LVs are gradually developed and applied in DA, ensemble generation and other areas in numerical forecasting. The NLLVs bear some similarities to the LVs; for example, they are orthogonal to each other and have different error growth rates from the largest to the smallest. However, they also have significant differences. The LVs are generated through a long-term integration using the tangent linear model and the direction of the leading LV is independent of the initial perturbation (Kalnay, 2002). In comparison, the NLLVs are obtained by integrating the nonlinear model from both the perturbed and unperturbed initial conditions, and represent different types of instabilities with different rescaling amplitudes of perturbations.

Both the BVs and the NLLVs are nonlinear extensions of LLVs, and can be dynamically produced by a sequence of two nonlinear model integrations with a periodic rescaling of their difference. BVs are random samples in the small fastest growing space, and have consistently large projection on the leading LLV (Magnusson et al., 2008). However, NWP models usually have a multidimensional growing subspace. Various independent fast growing directions should be comprehensively considered to span the unstable subspace (Duan and Huo, 2016). The fastest growing directions captured by the BVs have limited diversity and may not contain sufficient information to represent the development of analysis errors. In contrast, apart from the leading NLLV (LNLLV) being similar to the BVs, which tend to the fastest growing direction, other successively orthogonalized NLLVs point to additional fast growing directions in phase space. The NLLVs have larger diversity than the BVs and may capture a larger amount of components in analysis errors (Feng et al., 2014).

The NLLV method is related to the ETKF method, both of which are designed to orthogonalize perturbations. In the ETKF approach, the ensemble perturbation matrix is transformed by the matrix T associated with the covariance matrix of forecast members and observations (Wang and Bishop, 2003). In comparison, the NLLVs are orthogonalized using the Gram-Schmidt reorthonormalization (GSR) technique. Given that the ETKF scheme has been successfully implemented in operational forecast centers, it is worthwhile comparing the NLLV approach with the ETKF approach. The results of the BV and the random perturbation (RP) approaches are also given as reference, and all the comparisons are presented in a barotropic model. To clarify the differences between these four initialization methods, the experiments are implemented in a perfect model environment, and by using the same analysis and the same forecast model. An assimilation scheme independent of each ensemble initialization scheme —— namely, the ensemble Kalman filter (EnKF; Houtekamer and Mitchell, 1998; Evensen, 2003, 2004) —— is used to produce the initial analysis states.

This paper is organized as follows: Section 2 describes the BV, ETKF and NLLV ensemble generation schemes. An introduction to the barotropic model and the basic setup of the numerical experiments is provided in section 3. Section 4 presents the results of the quantitative comparison of the four ensemble methods. Section 5 provides a discussion and conclusions.

The breeding method proposed by Toth and Kalnay (1993, 1997) is used to simulate the error evolution in each analysis cycle. RPs are superposed on the analysis state to grow as a perturbation between two model integrations. At the end of each breeding cycle (12 hours in our experiments), the perturbations are rescaled by an appropriate factor. This, without the insertion of further random noise, defines a new set of perturbations that are used to generate the next perturbed initial states. This process is repeated and, as the effect of the particular initial RP diminishes, the resulting perturbations become the BVs. Initially, (Toth and Kalnay, 1993) used a global rescaling factor that is close to the empirically estimated global analysis RMSE. However, to ensure that the ensemble perturbation magnitude reflects the regionally varying uncertainties in the analysis, the breeding method used at NCEP includes an estimated mask as a constraint of BVs (Toth and Kalnay, 1997). Owing to our experiments having observations on every grid over the whole field and having the same precision, the distribution of analysis error variances is mainly determined by the forecast error variances. Therefore, only a simple global rescaling factor is applied in this study. The BVs are calculated using a Euclidean norm. Through a five-day breeding process, the generated BVs are tuned to the amplitude of the global analysis RMSE to use as the initial ensemble perturbations.

The ETKF perturbations are generated based on an ensemble DA algorithm (Wang and Bishop, 2003). We denote X' _f and X' _a as N forecast and analysis perturbations superposed on the control analysis. Forecast perturbations are listed as columns in the matrix X' _f and analysis perturbations are listed as columns in the matrix X' _a. The scheme is created to calculate a transformation matrix, T, which transforms the forecast perturbations according to \begin{equation} \label{eq1} {X}'_{\rm a}={X}'_{\rm f}{T} . (1)\end{equation} In this expression, T is an N× N matrix, and can be calculated as \begin{equation} \label{eq2} {T}={C}({\pmb{\it\Gamma}}+{I}_N)^{-1/2} ,(2) \end{equation} where columns of the matrix C contain the eigenvectors of (X' _f) ^TH ^TR^-1HX' _f/(N-1), and the diagonal matrix \(\pmb\it\Gamma\) contains the corresponding eigenvalues. H denotes the mapping from the model space to the observation space. R is the observation error covariance matrix. The ETKF ensemble perturbations can be proved orthogonal to each other in the space of the observations (Wang and Bishop, 2003).

The forecast perturbations in each cycle are multiplied by an inflation factor of 7% to maintain sufficient spread. The updated ensemble perturbations will be added to the control analysis to generate the initial conditions for the next filter step. The length of a cycle is 12 hours and the process is continued over 5 days to the initial analysis for forecasting. Each model grid has the same observation error variance of 0.02, i.e. the observation error covariance matrix R is a diagonal matrix with 0.02 entries.

Let us consider an n-dimensional nonlinear dynamical system: \begin{equation} \label{eq3} \frac{d}{dt}{x}(t)={F}[{x}(t)] , (3)\end{equation} where x(t)=(x₁(t),x₂(t),……,x_n(t)) ^T represents the state vector at time t governed by the dynamics F. Let δ(t₀) be a small error superposed on x(t₀) at an initial time t=t₀. The differential equation of the small error δ(t) is written as follows: \begin{equation} \label{eq4} \frac{d}{dt}{\delta}(t)={J}({x}(t)){\delta}(t)+{G}({x}(t),{\delta}(t)) , (4)\end{equation} where J(x(t))δ(t) and G(x(t),δ(t)) are the tangent linear terms and the high order nonlinear terms of the errors δ(t), respectively. By integrating Eq. (5) from t=t₀ to t₀+τ, the nonlinear error evolution can be obtained: \begin{equation} \label{eq5} {\delta}(t_0+\tau)={\eta}({x}(t_0),{\delta}(t_0),\tau){\delta}(t_0) ,(5) \end{equation} where η(x(t₀),δ(t₀),τ) is the nonlinear propagator (Ding and Li, 2007). In a chaotic system, each initial error vector tends to fall along the fastest growing direction. Therefore, for any initial error vector δ'(t₀), after a short time τ, δ'(t₀+τ) will capture the fastest growing direction. If the short time τ is not taken into consideration, taking δ'(t₀+τ) as the initial error δ₁(t₀), δ₁(t₀) is defined as the LNLLV here, which is similar to the BV. The first nonlinear local Lyapunov exponent (NLLE) along the fastest growing direction LNLLV can be approximately defined as \begin{equation} \label{eq6} \lambda_1({x}(t_0),{\delta}_1(t_0),\tau)=\frac{1}{\tau}\ln\frac{\|{\delta}_1(t_0+\tau)\|}{\|{\delta}_1(t_0)\|} , (6)\end{equation} where Λ₁(x(t₀),δ₁(t₀),τ) is the function of the initial state x(t₀) in phase space, the initial error δ₁(t₀), and evolution time τ. In contrast, the linear Lyapunov exponent relies on the assumptions of infinite evolution time and infinitesimal initial error (Lorenz, 1965; Yoden and Nomura, 1993; Ziehmann et al., 2000), and thus has limitations when applied to nonlinear error growth. The NLLE has been widely used to investigate the predictability of weather and climate (Li et al., 2006; Ding and Li, 2007; Ding et al., 2008, 2010, 2011; Li and Ding, 2011, 2013).

While the LNLLV tends to the fastest growing direction, the rest of the NLLVs can be derived via a periodic reorthogonalization by the GSR process (Wolf et al., 1985; Li and Wang, 2008, Feng et al., 2014). The orthogonalized perturbations are then globally rescaled to the initial size and again superposed on the new analysis to integrate (see Fig. 1). The process is performed for five days with a 12-hour cycle to acquire the NLLVs denoted by δ₁(t₀),δ₂(t₀),…,δ_N(t₀). This 5-day period is basically close to that used in (Vannitsem and Nicolis, 1997) to generate LVs in a quasigeostrophic model. The NLLVs are orthogonal to each other, representing the vectors along the directions from the fastest growing direction to the fastest shrinking direction in local phase space. The resulting NLLV perturbations are tuned to the amplitude of the global analysis RMSE and superposed on the analysis to generate ensemble members.

Once the NLLVs are obtained, the growth rate of perturbations along each direction can be estimated by the NLLEs, which are defined by \begin{equation} \label{eq7} \lambda_i({x}(t_0),{\delta}_i(t_0),\tau)=\frac{1}{\tau}\ln\frac{\|{\delta}_i(t_0+\tau)\|}{\|{\delta}_i(t_0)\|} , (7)\end{equation} where δ_i(t₀) is the i-th NLLV at time t₀. If we want to study the dynamic characteristics of the whole system, the ensemble average of the i-th NLLE, \(\overline \lambda_i(\delta_i(t_0),\tau)=\langle\lambda_i(x(t_0),\delta_i(t_0),\tau)\rangle\) (the symbol \(\langle\ \ \rangle\) denotes the ensemble average over a great number of different initial states) should be introduced.

The background forecast errors inserted in the analysis errors are the combined effects of the multidimensional growing subspace (Palatella and Trevisan, 2015) spanning the NLLVs. Therefore, to better sample the analysis errors, the initial ensemble perturbations are constructed by a simple random linear combination of the first N NLLVs (N=10). The NLLVs are denoted by \(\pmb\it\Delta=\delta_1,\delta_2,\ldots,\delta_N\), then the orthogonal ensemble perturbations \(\pmb\it\Delta\) can be calculated by \begin{equation} \label{eq8} {\pmb{\it\Delta}}_{\rm e}={\pmb{\it\Delta}}{\beta} , (8)\end{equation} where β is an N× N matrix. The column vectors of β are constructed by random numbers of a standard normal distribution and are independent from each other. The resulting ensemble perturbations are rescaled to the amplitude of the global analysis error and superposed on the analysis to generate ensemble members.

Figure 1.  Creation of NLLVs. A group of RPs is initially introduced on the analysis state, and integrated to the end of a breeding cycle. The orthogonalization processes are conducted regularly. The red solid curve represents the leading perturbed forecast during each cycle, and the black solid curves represent other perturbed forecasts from the second to the $n$th fastest. The direction of the fastest growing mode (LNLLV, red dashed line) is kept, and the other directions (black dashed lines) are orthogonalized successively to acquire NLLV$_{2}$ to NLLV$_{n}$ (blue dashed line). These perturbations are then rescaled to the same size as the initial perturbations, and added to the new analysis. After several days of cycling, the perturbations are dominated by growing components [modified from Feng et al. (2014)]

The two-dimensional quasi-geostrophic model is governed by the following system of dimensionless differential equations: \begin{equation} \label{eq9} \left\{ \begin{array}{l} \dfrac{\partial q}{\partial t}+J(\psi,q)=0\\[3mm] q=\nabla^2\psi-F\psi+f_0+f_0h_s,\quad {\rm in}\quad \Omega\times[0,T]\\[1mm] \psi|_{t=0}=\psi_0 \end{array} \right. ,(9) \end{equation} where the variables q and ψ are the potential vorticity and the stream function, respectively; ∇²=∂²/∂ k²+∂²/∂ l² and J(ψ,q)=ψ_kq_l-ψ_kq_l are the Laplacian operator and the horizontal Jacobian operator, respectively; k and l represent the zonal and meridional coordinates; and t denotes time. The planetary Froude number F=0.102, and the Coriolis parameter f₀=10.0. A double periodical boundary condition is used for \(\Omega=[0,K]\times[0,L]\). The barotropic model has been used in various studies of predictability and ensemble forecasts (Barkmeijer, 1992; Cheung and Chan, 1999; Mu and Zhang, 2006; Mu and Jiang, 2008; Durran and Gingrich, 2014).

All the experiments below are performed on a 32× 16 latitude-longitude horizontal grid. The grid spacing d=0.2 corresponds to a dimensional length of 200 km, and the time step dt=0.006 corresponds to 10 min. The basic flow is obtained by integrating (10) with an initial state \(\psi_0=0.5\sin(\pi k/6.4)+\sin(\pi l/3.2)+0.5\), as used by (Mu and Jiang, 2008). For simplicity, the topography is only a function of the l-direction \(h_s=h_0\times[\sin(2\pi l/3.2)+1]\).

The quasi-geostrophic model is integrated for 800 days as a "truth run", after a first "spin-up" run of 100 days. The quality of the forecasts are assessed by comparing with the truth run. For a fair comparison, the EnKF assimilation method, which is irrelevant to the ensemble generation schemes, is used to generate the control analysis. The observations are distributed at each grid point, and are produced by adding a random noise to the true state. The assimilation process is implemented every 12 hours on the stream function. The details of the EnKF assimilation steps are described in Appendix A.

Once the analysis states have been obtained, 10-day ensemble forecasts initialized from the four ensemble generation methods, respectively, are implemented. All the initial ensemble perturbations are produced using the Euclidean norm. The RP perturbations are simply generated by random numbers with a standard normal distribution, and then rescaled to the amplitude of initial analysis errors. The ensemble perturbations of the RP scheme are flow-independent and thus can provide a proper reference for the other theory-based schemes (Magnusson et al., 2009). To clarify the differences of the four initialization schemes in a simple way, in this study, the N=10 ensemble perturbations for the four schemes are directly superposed on the analysis state, without being centered, to generate ensemble members. The statistical quality of the forecasts is assessed by averaging a set of 100 different experiments, each corresponding to a different initial state randomly chosen from the 800-day truth run. All experiments are "perfect model" experiments in which the assimilation and forecasts are performed with the same model used for producing the truth run.

Before using the barotropic model to compute the NLLEs and NLLVs, it should be ensured that this model is a nonlinear dynamical system and has a chaotic attractor (Li and Chou, 1997). Figure 2 shows the temporal evolution of the stream function at one grid point for 800 days. It can be seen that the trajectory presents a chaotic behavior, which indicates the model is suitable for the research on nonlinear error growth.

To see if different growing directions can be effectively captured by the NLLVs, the transient error growths of NLLVs are measured by the NLLEs. Figure 3 shows the first 50 NLLEs averaged over 800 initial conditions with a one-day interval. The NLLEs are calculated for the error growth during the initial τ=12 h and the initial amplitude of perturbations is 0.01. It can be seen that the first 50 NLLEs are consistently positive and continuously arranged in a decreasing order. The results demonstrate that the barotropic model has a high-dimensional growing subspace. The development of the analysis error is a combined effect of various independent growing directions, and the NLLVs may be capable of capturing the unstable subspace effectively.

Figure 2.  The evolution of the stream function at one grid point of the barotropic model for 800 days.

Figure 3.  The first 50 NLLEs for the barotropic model with the evolution time $\tau$ equal to 12 hours and a finite initial amplitude of 0.01

Since analysis errors are poorly known, the initial ensemble perturbations should be as independent as possible to sample the possible direction of analysis errors. To see the independence of the global perturbations, the average variances explained by the eigenvectors of the ensemble covariance matrix are calculated for the four schemes, respectively (Wang and Bishop, 2003; Bowler, 2006). The assessments for the ensemble perturbations superposed on the control forecasts at the different lead times are shown in Fig. 4. A flatter distribution of the explained variance indicates better perturbation independence. At day 0, the explained variances for the NLLV and ETKF perturbations are strictly evenly distributed, which can be induced by their theoretical concepts. The RP perturbations are very close to being orthogonal to each other. However, the first eigenvector of the BVs explains almost 30% of the error variance when 10 perturbations are used, indicating that there is certain similarity among the perturbations. That is to say, the effective number of ensemble perturbations for the BV method is lower than its actual number. With the increase in lead time, the similarities among the ensemble perturbations for the four schemes gradually increase due to the dynamical evolution on the same basic flows. The performances of the four initialization schemes tend to be similar after about day 6.

Figure 4.  Mean explained variances of the eigenvectors of forecast covariance matrices measured from day 0 to 10. The BV, NLLV, RP and ETKF methods are denoted by triangles, open circles, squares and crosses, respectively.

Figure 5.  Evolution of mean error variance of the control run explained by the ensemble perturbations as a function of lead time.

Initial ensemble perturbations are created to sample the possible forecast errors. Therefore, the variance of forecast errors explained by the subspace of ensemble perturbations is an important assessment of the quality of perturbations. For a fair comparison of the four schemes, the ensemble forecast errors are defined as the difference between each ensemble forecast member and the control forecast, while the control forecast error is the difference between the true state and the control forecast. Let us assume the normalized ensemble forecast errors are denoted by E=e₁,e₂,…,e_N, and the control forecast error is e₀. The projection of the control forecast error on the subspace of ensemble forecast errors would be derived by p=EE ^Te₀. Then, the error variance of the control forecast explained by the ensemble perturbation subspace can be calculated through the correlation between p and e₀.

Figure 5 illustrates the variance of analysis errors explained by the ensemble perturbations superposed on the control forecasts as a function of lead time for the four methods. At the initial time, the RP approach performs the worst in terms of capturing the analysis error direction, as the RPs are independent of the flow at initial time. The NLLV and ETKF perturbations have similar explained variance, which is larger than that of the BV perturbations. This may be due to the more independent perturbations of the former two approaches than the latter one (see Fig. 4). The initial explained variances of the four schemes may also cause their relative performance afterwards. The ETKF and NLLV techniques basically remain similar and the best from day 0 to 10. The BV approach always has a lower score compared to the ETKF and NLLV approaches during the 10 days. The explained variance for the RP, despite exhibiting the fastest increase among the four schemes, has the worst performance throughout.

The RMSE and the pattern anomaly correlation (PAC) of the ensemble mean are often used to measure the overall forecast performance of ensemble forecasts. Figure 6 illustrates the RMSE and PAC as a function of time for the four schemes. The ensemble spread, calculated as the standard deviation of ensemble members, is also given as a reference. The ensemble spread of a statistically reliable ensemble system is usually close to the RMSE of the ensemble mean for all lead times (Buizza et al., 2005). From Fig. 6a, it can be seen that for the first 4 days the NLLV, RP and ETKF ensembles have similar RMSE scores, which are slightly smaller than that of the BV scheme. This may because the BVs are constrained to a smaller perturbation subspace compared to the other three types of initial perturbations (see Fig. 4). Note that the ensemble spread of the RP approach is smaller than the other three for 0-4 days. This is due to the RPs pointing to non-growing directions and will shrink initially (Toth and Kalnay, 1997). From day 4 to 10, the NLLV has slightly smaller RMSE and similar ensemble spread compared with ETKF, while the BV scheme has somewhat larger RMSE and smaller spread than both the NLLV and ETKF schemes. This is because the NLLV perturbations, which consider the nonlinear interactions of growing directions, have a better sampling of the development of the analysis errors than the ETKF perturbations, and both of them span more orthogonal and uncorrelated subspaces than do the BVs. Although the RP approach has the largest ensemble spread, its forecast skill decreases rapidly and remains the lowest among the four schemes. This can be explained by the RPs having very limited ability to capture the evolution of initial analysis errors. The ensemble spread for both the NLLV and ETKF methods are closer to their respective ensemble mean errors than the BV method, which implies that the former two have better statistical reliability than the latter one. The PAC score in Fig. 6b presents similar relative performance of the ensemble mean skill as the RMSE.

Figure 7 shows the ratio of the highest ensemble mean skill among the four initialization schemes for each of them as a function of lead time. It can be seen that the ETKF and NLLV approaches basically have similar ratios to reach the highest forecast skill, which is evidently higher than that of the BV and RP approaches. The RP approach performs better than the BV approach from day 1 to 4. Thereafter, the ratio of the RP scheme remains the lowest among the four schemes to the end of the lead time. The statistical results for each case basically correspond to the sample mean results.

Figure 6.  Evolution of the (a) mean RMSE (solid lines) and ensemble spread (dashed lines) and (b) PAC (solid lines) for the BV, NLLV, RP and ETKF schemes as a function of lead time.

Figure 7.  Evolution of the ratio of the (a) smallest RMSE and (b) highest PAC among the BV, NLLV, RP and ETKF initialization schemes for each of them as a function of lead time.

A commonly used measure in probabilistic forecasts is the Brier score (BS), which can be used to evaluate the reliability (BS _rel) and resolution (BS _res) of an ensemble forecast system (Murphy, 1973; Descamps and Talagrand, 2007; Stephenson et al., 2008). The score assesses the performance of the probability forecast of the occurrence of an event Φ. Further details are given in, for example, (Descamps and Talagrand, 2007) and (Stephenson et al., 2008). In this paper, Φ is defined as the event ψ>ψ _c, where ψ _c is the climatological mean at each grid point. A smaller value of BS and BS _rel indicates a better forecast skill; while the BS _res acts in an opposite way. Another event interval (ψ>ψ _c+σ, where σ is the standard deviation of ψ) is also introduced for diagnostics, but does not appear to affect the relative performance of the ensemble forecasts. Only the results of the event Φ are shown in this paper.

Figure 8a shows the basic BS evaluated for the four methods as a function of lead time. The four methods have similar basic BS scores during days 0-4. The NLLV scheme shows slightly better performance than the ETKF scheme. The BV scheme has a lower basic BS than the RP scheme for days 4-10, and both of them perform worse than the ETKF and NLLV schemes. Their differences in BS performance mainly come from the BS _res and the BS _rel components, as shown in Fig. 8b. After 4 days, the higher BS _res of the NLLV and the ETKF schemes than the BV scheme may be attributable to the worse sampling of the analysis errors by the BV approach. Meanwhile, the former two approaches have better independence and larger ensemble spread than the BV approach, which may result in a lower BS _rel relative to the BV approach. The RP scheme performs worst among the four schemes generally. This is because the RP ensemble is randomly generated, and thus is less likely to describe the probability distribution of the truth relative to the other three schemes.

Figure 8.  Evolution of (a) average basic BS as a function of lead time for the BV, NLLV, RP and ETKF schemes; and (b) resolution (dotted lines) and reliability (solid lines) contributions to the BS.

The NLLV ensemble generation scheme retains the advantage of the BV scheme, which uses the breeding cycles to capture the development of analysis errors in the assimilation cycles and is very simple and cheap to implement. Different from the BVs, the NLLV perturbations are periodically orthogonalized by the GSR process to identify various distinguished growing perturbations, and thus have more diversity and independence. In this paper, the application of the NLLV method introduced by (Feng et al., 2014) for the generation of ensemble perturbations is further extended to a moderately complex barotropic model. Except for the BV and RP approaches, the performance of the NLLV scheme for ensemble forecasting is further compared with that of the ETKF scheme, which is one of the most advanced ensemble initialization schemes. It is found that the NLLV scheme performs slightly better overall than the ETKF scheme for various scores when the observations to be assimilated are evenly distributed at each grid. This may because the development of analysis errors can be effectively sampled by the unstable perturbations acquired by the NLLVs. The BV approach performs worse than the ETKF and NLLV approaches due to the dependence among perturbations, but still dramatically exceeds the RP approach.

This study simply uses a barotropic model without model error, and the model size is less than that of operational forecasts. However, the quasi-geostrophic model is sufficiently large to demonstrate the relative performance of the ensemble schemes accurately. Despite the similar forecast skills of the NLLV and ETKF schemes, the generation of the NLLVs is significantly time-saving and easy to implement, compared to the ETKF scheme. The computational expense of the former is only about one third of the latter in the experimental environment of this paper. Meanwhile, the quality of ETKF perturbations largely depends on the estimation of the background and observation covariance matrices, while an accurate estimate of the forecast error covariance is very challenging due to such problems as ensemble "collapse" (Anderson and Anderson, 1999) and spurious correlations at large distances (Houtekamer and Mitchell, 2001). In a real world situation where the observations are irregularly distributed, a prior knowledge of the regional rescaling factors could be used to modulate the NLLV perturbations to better sample the analysis error variances. Several current approaches can provide the estimated masks, like the error growth statistical method (Pe\ na and Toth, 2014). Introducing the estimation of analysis error variances in tuning the NLLV perturbations may further improve the ensemble forecast skill of the NLLV scheme.

The Ensemble Kalman Filter

The true state is denoted by ${x}_{\rm true}$, which is derived from a long-term model integration. The simulated observations ${y}$ are generated based on the true state using \begin{equation} {y}={Hx}_{\rm true}+{\varepsilon}~,(A1)\end{equation} where ${H}$ provides a mapping from the model space to the observation space, and ${\varepsilon}$ represents independent realizations of the noise with a Gaussian distribution $N(0,0.02)$. The ensemble of forecast states is adopted as background states. The ensemble matrix is then defined as \begin{equation} {X}_{\rm f}=({x}_{\rm f,1},{x}_{\rm f,2},\ldots,{x}_{\rm f,{N}})~. (A2)\end{equation} where $\overline{{x}}_{\rm f}$ represents the mean of the ensemble. An ensemble perturbation matrix can then be written as \begin{equation} {X}'_{\rm f}=({x}_{\rm f,1}-\overline{{x}}_{\rm f},{x}_{\rm f,2}-\overline{{x}}_{\rm f},\ldots,{x}_{\rm f,{N}}-\overline{{x}}_{\rm f})~. (A3)\end{equation} The covariance matrix of the ensemble ${X}_{\rm f}$ is: \begin{equation} {P}_{\rm f}=\frac{1}{N-1}{X}'_{\rm f}{X}'^{\rm T}_{\rm f}~. (A4)\end{equation} The background forecast ensemble is to be updated by the observations. The set ${y}_i(i=1,2,\ldots,N)$ is a set of perturbed observations that are associated with each previous forecast of the ensemble. They are defined as: \begin{equation} {y}_i={y}+{\varepsilon}_i~. (A5)\end{equation} The RPs ${\varepsilon}_{i}$ follow the same distribution as ${\varepsilon}$. The observations are assimilated to produce a new analysis of the state: \begin{equation} {x}_{{\rm a},i}={x}_{{\rm f},i}+{K}({y}_i-{Hx}_{{\rm f},i}),\quad i=1,2,\ldots,N~. (A6)\end{equation} The Kalman gain ${K}$ is calculated by: \begin{equation} {K}={P}_{\rm f}{H}^{\rm T}({HP}_{\rm f}{H}^{\rm T}+{R})^{-1}~. (A7)\end{equation} This is actually a weight measuring the ratio of the forecast and observational error covariance. ${R}$ is the observational error covariance matrix.

To overcome the problem of undersampling, a 7{\%} variance inflation factor is applied to ${X}'_{\rm f}$ in this barotropic model. Moreover, the number of the ensemble is 100, which is much smaller than the model dimension. Therefore, the localization technique is applied to the matrix ${P}_{\rm f}$ to prevent spurious correlations at large distances. This is realized by the fifth-order function of \cite{Gaspari1999} with a distance of zero correlation equal to 800 km (four grids). The assimilation cycles are repeated for 30 days in each case to generate the analysis ensemble. The mean $\overline{{x}}_{\rm a}$ of the analyzed ensemble ${x}_{{\rm a},i}$($i=1,2,\ldots,N$) is regarded as the initial analysis state when performing the forecasts.