An Analytical Four‐Dimensional Ensemble‐Variational Data Assimilation Scheme

Abstract

The usage of four‐dimensional variational (4DVar) scheme is limited by the static background error covariance and the adjoint model. In a hybrid frame of the four‐dimensional ensemble‐variational data assimilation scheme (4DEnVar), being able to avoid the tangent linear and adjoint models in the 4DVar and nowadays developed into a cutting‐edge research topic of the next‐generation data assimilation methods, an analytical 4DEnVar (A‐4DEnVar) scheme is designed. First, an analytical expression for explicit evolution of the background error covariances is derived. The expression collects the innovation of observations over an assimilation window simultaneously and propagates information to the initial background field by temporal cross covariances. Second, to estimate the adjoint model, the temporal covariances are constructed with ensemble members being centralized with respect to the model states integrated from the initial condition. Third, an iterative linear search process is introduced to minimize the cost function to update the analysis field until convergence. Twin experiments based on the Lorenz chaos model with three variables are conducted for the validation of the A‐4DEnVar scheme. Comparisons to the conventional 4DVar show that the A‐4DEnVar is comparable in accuracy even with a long assimilation window and sparse observations. The assimilation results also show that the A‐4DEnVar scheme can be implemented with a very small ensemble size which means that under circumstances without the tangent linear and adjoint models it can be easily incorporated into data assimilation systems in use.

Introduction

A better model initial condition combining the information from observations can improve the skill of weather or ocean forecasts dramatically. To this end, data assimilation schemes are proposed to transport information from sparse observations to the gridded model fields. To connect observations with state variables, variational schemes, and Kalman filter schemes have been widely adopted. Under similar statistical interpretations, most data assimilation schemes rely on Gaussian processes and optimal linear estimation theory due to their simplicity in reality. Background error covariances not only indicate the relationship between state variables, but also control the propagations of information through each grid points, thus playing a significant role in the data assimilation schemes (C. Liu et al., 2008).

Despite its importance, and due to its huge dimensions and complex dynamics of state variables, covariances cannot be always easily estimated. In variational data assimilation schemes, the background error covariance is static. Traditionally, the error covariances are assumed to be a function of vertical and horizontal distance. The parameters of these functions can be inferred by using an innovation method (Xu et al. 2001, 2002), NMC method or ensemble methods (Buehner, 2005; Wang et al., 2014). Bannister (2008a, 2008b) discussed these numerical methods in details. As an efficient approximation of a Gaussian filter, recursive filters, and diffusion filters can be employed to set up inhomogeneous and anisotropic background error covariances (Derber et al., 2003; Zhang et al.2015). A more effective construction has been done through multiscale method (Li et al., 2010). For these schemes mentioned above, the estimated covariances disregard their flow‐dependent property (without evolutions in time) and/or are constructed empirically (covariances have been abstracted to functions, explicitly or implicitly).

In data assimilation, to make accurate estimations of the initial fields, one should also utilize all available information from current and past observations. The connection between different time steps is hence extremely important. In four dimensional variational (4DVar) schemes, the evolution of the errors of state variables are reflected through the tangent linear and adjoint models (Thépaut et al., 1993). The initial state variables would be integrated forward to the observation times to obtain the innovation, and then integrated backwards by using the adjoint model to get the optimal solution by minimizing the cost function. The error covariances of states at different time are determined implicitly by the adjoint models in 4DVar (Amezcua et al., 2017).

In contrast to variational methods, sequential methods such as Kalman filter and extended Kalman filter (EKF, Evensen, 1992) describe covariances explicitly. The forecast steps and analysis steps blend information from the dynamic model and the observations alternately. Because of their large dimensionality, it is too expensive to apply an EKF to operational systems. Monte‐Carlo approaches such as the well‐known nonlinear data assimilation scheme, i.e., the ensemble Kalman filter (EnKF, Evensen et al., 1994) has been recommended to alleviate the problem. Unlike other schemes, in EnKF, the flow‐dependent covariances can be estimated explicitly. With more ensemble members and better localization and inflation technologies, the covariance would be expected to be more accurate. The other most conspicuous advantage of EnKF schemes is that they do improve the background error covariances but do not rely on the tangent linear model and its adjoint, which are indispensable in 4DVar schemes and cost too much in computational and development resources.

Though the EnKF based schemes are competitive for the estimation of the covariances of nonlinear process, there is insufficient evidence that EnKF schemes can provide a better analysis or forecast than variational schemes. Currently, the variational method is still the preferred choice of leading operational numerical weather prediction (NWP) centers. Both schemes have merits and demerits, and comparisons have been discussed in Bannister (2017). An attractive issue is to combine the advantages on both sides. The hybrid EnKF‐three‐dimensional variational (3DVar) schemes have been introduced and improved by Hamill (2000), Oddo (2016), and Storto et al. (2018) respectively, who decomposed background covariances to the time‐invariant 3DVar static covariances and flow‐dependent covariances to improve the estimation. These methods provide a better background error covariance but do not avoid the adjoint model in 4DVar schemes.

The high computational cost of the inverse of the covariance matrix and the human resources needed for coding the tangent linear or adjoint tangent linear models result in heavy and complex operations, especially when the data assimilation system evolves or is transmitted to other systems with different dynamic models used. In the four‐dimensional ensemble‐variational data assimilation scheme (4DEnVar), the background error covariances in 4DVar can be replaced by those implied by ensemble members. While spatial covariances of state variables indicate information about structures at a specific static observation time, the temporal cross covariances between different observation times point out the information transformation structure through time. The form of the covariance estimated from ensemble members makes it possible to utilize an incremental method to produce an algorithm of 4DEnVar in which adjoint model involved in the traditional 4DVar can be avoided. Introduced by C. Liu et al. (2008, 2009), the 4DEnVar was first named as ensemble‐based four‐dimensional variational data assimilation scheme (En4DVar), and subsequently renamed as 4DEnVar (C. Liu & Xiao, 2013; C. S. Liu & Xue, 2016) to highlight its features that the tangent linear model and its adjoint model are omitted (Lorenc, 2013). By using real observations in the Canadian operational global model, Buehner et al. (2010a, 2010b) tested 4DEnVar. Furthermore, studies were also proposed to compare the performance of ensemble based variational data assimilation schemes through the Lorenz model (Goodliff et al., 2015) and an NWP model (Lorenc et al., 2015). Amezcua et al. (2017) demonstrates the detailed description of strong‐ and weak‐constraint 4DEnVar using a Korteweg‐de‐Vries equation.

In C. Liu et al. (2008), using a preconditioning variable w and ensemble perturbations divided by $\sqrt{N-1}$ (here, denoted as ${{\mathbf{X}}^{\prime }}_b$), the analysis state x_a was transformed to ensemble space, i.e., [Image Omitted. See PDF] [Image Omitted. See PDF]where x_b is the ensemble member state vector and ${\overline{\mathbf{x}}}_b$ is the ensemble mean, and [Image Omitted. See PDF]where H is the observation operator, O denotes the observational error covariance, M is the tangent linear model and d is the innovation. The use of a preconditioning variable is a brilliant option to avoid the calculation of the background error covariance and also to avoid the calculation of the adjoint model. However, Equation 17 of C. Liu et al. (2008) demonstrated that, the background error covariances was constructed using the ensemble perturbations being centralized with respect to the ensemble mean state (denoted with a bar, as shown below), as used in the traditional ensemble‐based filters, i.e., [Image Omitted. See PDF]

Subsequently, an optimal algorithm could be employed based on the gradient of the cost function [Image Omitted. See PDF]to achieve the analysis increment. On the one hand, the covariances in 4DEnVar are sample covariances that centralized to the solution evolved from ${\overline{\mathbf{x}}}_b$ but not the x_b itself. As an estimation of the covariances of error in state variables, the sample covariances are not suitable. A better construction method is to centralize ensemble members with respect to the model state evolved from the initial condition at each iteration step like the conventional 4DVar does. These two kinds of covariance matrix are similar for linear models, but are potentially great differences for nonlinear models. On the other hand, though the adjoint model in 4DEnVar has been dropped explicitly, it is in fact represented implicitly by ensemble members. As ensemble members were generated near the background state, the implicit adjoint model was not updated during iterations, hence the ability of 4DEnVar was further restricted when applied to nonlinear models. To overcome this disadvantage, the ensemble members should be regenerated near the model state evolved from the initial condition after each iteration to track the model nonlinearity.

Based on a linear approximation of the dynamic model, the widely used tangent linear and adjoint models in 4DVar can be rewritten through temporal cross covariances. In what follows, we will show that the minimization of the cost function through an analytical solution is in fact feasible if the covariances are perfectly known. Based on the solution, and with finite ensemble size, an adjoint model free data assimilation method named analytical 4DEnVar (A‐4DEnVar, henceforth) will be introduced and evaluated. Different from previous 4DEnVar schemes, the proposed A‐4DEnVar centralizes ensemble members with respect to the model state evolved from the initial condition, and updates temporal covariances (corresponding to the adjoint model as 4DVar does) and the guessed initial states iteratively, so that nonlinearity is well considered.

The study is organized as follows: in Section 2, we describe the methodology of the study, including the formulation of the cost function based on temporal cross covariances and the minimization algorithm of the cost function based on finite ensemble size. In Section 3, validation experiments are described in details. Comparison to the conventional 4DVar is presented through a simple perfectly identifiable model, i.e., the Lorenz three‐variable chaotic model (Lorenz, 1963, Lorenz‐63 hereafter). The last section gives discussions, conclusions and directions for our future work.

An A‐4DEnVar

Evolution of Error Covariance

4DVar can simultaneously assimilate observations over a given time window by using a tangent linear model and its adjoint. The A‐4DEnVar associates these observations with model states through temporal cross error covariances that are calculated to replace the adjoint model. As a beginning, we review the basic general context of data assimilation and discuss the evolution of error covariances to rebuild the cost function.

In strong‐constraint 4DVar, x_i the model state vector which contains the state variables as elements at time t_i, must be a solution of the nonlinear model equations, i.e., [Image Omitted. See PDF]where ${\mathrm{M}}_{i,i-1}\left({\mathbf{x}}_{i-1},{\mathbf{F}}_{i-1}\right)$ denotes the forecast model operator from time t_i−1 to time t_i. The model's arguments consist of two components, the state vector x_i−1, and external forcing component F_i−1. In numerical ocean models, F_i−1 may include wind stress, heat flux, fresh water flux and open boundary conditions, etc.

In practice, we do not know the truth of the state and external forcing. As an approximation to the truth, a known model state ${\mathbf{x}}_i^{\ast }$ (at the beginning of data assimilation process, this is the background state), is a nearby estimation. Thus x_i can be decomposed into the known state and difference or perturbation (here distinguished with a superscript star and a tilde, respectively), [Image Omitted. See PDF]

The separation above makes it easy to design an iteration process. In the same way, the external forcing component can also be decomposed into a known forcing and its perturbation [Image Omitted. See PDF]

Substitution of the above decompositions into Equation 1 and linearizing the model operator based on the Taylor expansion leads to [Image Omitted. See PDF]where operators ${{\mathbf{M}}_{i,i-1}|}_{{\mathbf{x}}_{i-1}^{\ast },{\mathbf{F}}_{i-1}^{\ast }}$ and ${{\mathbf{M}}_{i,i-1}^{\mathbf{F}}|}_{{\mathbf{x}}_{i-1}^{\ast },{\mathbf{F}}_{i-1}^{\ast }}$ are the tangent linear expansion of ${\mathrm{M}}_{i,i-1}$ with respect to the model state and forcing component, respectively. The higher order moments are denoted as ε and could be removed if the perturbations are small enough. It is clear that the perturbation satisfies [Image Omitted. See PDF]where forcing terms in the subscript are omitted and ${\stackrel{\sim }{\mathbf{f}}}_{i-1}\equiv {{\mathbf{M}}_{i,i-1}^{\mathbf{F}}|}_{{\mathbf{x}}_{i-1}^{\ast },{\mathbf{F}}_{i-1}^{\ast }}\left({\stackrel{\sim }{\mathbf{F}}}_{i-1}\right)$ is introduced to simplify the notation.

Equations 1 and 3 show that the evolution of the model state and perturbation have their own structures. The model state relies on the nonlinear model while the perturbation relies on the linearized model in data assimilation. However, the linearized model is sensitive to the known state ${\mathbf{x}}_i^{\ast }$. As we emphasized before, one must note that the Taylor expansion is based on the known approximation of truth states. When calculating the temporal covariance, the perturbations are the difference between ensemble member states and the state evolved from the known approximation of the true initial conditions after each iteration, but not the states evolved from the average of the ensemble initial conditions.

To simplify the above formulas, we introduce symbol notation ${{\mathbf{M}}_{j,i}|}_{{\mathbf{x}}_i^{\ast }}$ to illustrate a multistep evolution of the perturbation from time t_i to t_j, [Image Omitted. See PDF]here, the evolution operator ${{\mathbf{M}}_{j,i}|}_{{\mathbf{x}}_i^{\ast }}$ depends not only on the model states at time t_i and t_j, but also on any states over the time interval. Using the multistep tangent linear model, the perturbation state at any time t_i can be related to the initial state perturbation. It is easy to see that, the perturbation states in Equation 3 being evolutions of the initial state perturbation should be reformulated as [Image Omitted. See PDF]

Equation 5 holds for any small perturbation ${\stackrel{\sim }{\mathbf{x}}}_i^k$ (k = 1,2,⋯,N, where N is the ensemble size), i.e., [Image Omitted. See PDF]

Generally speaking, the initial state and the forcing fields can be considered uncorrelated to each other so that if the ensemble size is large enough, the adjoint model satisfies [Image Omitted. See PDF]

In this study, to obtain the initial condition of the ensemble member (that is, ${\mathbf{x}}_0^k={\mathbf{x}}_0^{\ast }+{\stackrel{\sim }{\mathbf{x}}}_0^k$), we add Gaussian noise to the updated initial condition. These initial conditions are available to ensure the validity of Equation 6 if they are near to the updated initial condition. Then the ensemble members integrated using these ensemble initial conditions are centralized to the background rather than the ensemble mean as 4DEnVar does. To alleviate the huge computational cost, the covariance of the Gaussian noise is set to be proportional to background error covariance (see more details below).

We introduce the subscripts to distinguish the background error covariance B with temporal error cross covariances of ensemble members appeared in Equation 6. That are [Image Omitted. See PDF]

The temporal error cross covariances with subscript zero reveal the interactions between the initial states and time trajectories of the model. Using the new notations, Equation 6 can be simplified as ${\mathbf{B}}_{i0}\approx {{\mathbf{M}}_{i,0}|}_{{\mathbf{x}}_0^{\ast }}{\mathbf{B}}_{00}$and it is easy to see ${\mathbf{B}}_{0j}\approx {\mathbf{B}}_{00}{{\mathbf{M}}_{j,0}^T|}_{{\mathbf{x}}_0^{\ast }}$. The temporal error cross covariances evolved from the background field contain the information of evolutions of the model states, while in a traditional 4DVar scheme, the relationships are represented by tangent linear and adjoint models. But in our study, the adjoint operator is implicated in the temporal error cross covariance matrix, which will be seen in the following derivations.

Cost Function and Analytical Minimization

Based on temporal error cross covariances it is possible to construct an A‐4DEnVar data assimilation scheme without tangent linear or adjoint models. The conventional cost function of a 4DVar assimilation scheme with respect to the initial state x₀ consists of the background and the observation terms, i.e., [Image Omitted. See PDF]where x_b is the background state, H is the observational operator, R is the observation error covariance matrix and y_i denotes the observation at time t_i.

After the forcing term is prescribed, the model state at time t_i can be expressed as [Image Omitted. See PDF]

Since ${\mathbf{x}}_i^{\ast }$ is already known, the cost function can be rewritten as a function of perturbation as follows [Image Omitted. See PDF]

To achieve the minimization of the cost function, in practice, optimization algorithms such as Newton iterative method and steepest descent method are required. In order to see analytically the solution that an iterative algorithm produces, setting the gradient of Equation 10 equal to zero gives, [Image Omitted. See PDF]

Furthermore, by using equations ${\mathbf{B}}_{i0}\approx {{\mathbf{M}}_{i,0}|}_{{\mathbf{x}}_0^{\ast }}{\mathbf{B}}_{00}$ and ${\mathbf{B}}_{0j}\approx {\mathbf{B}}_{00}{{\mathbf{M}}_{j,0}^T|}_{{\mathbf{x}}_0^{\ast }}$, Equation 11 is equivalent to [Image Omitted. See PDF]

Similar with the solution of 4DVar schemes, the perturbation states are combinations of the background and weighted observation innovations, but all weights can be represented by spatial covariances (that are covariances with two same subscripts) and temporal error cross covariances (that are covariances with two different subscripts). Hence the analysis state of the initial field is [Image Omitted. See PDF]

The analytical solution consists of component ${\mathbf{x}}_0^{\ast }$ as its iterative initial, a combination of weighted differences from the background state and innovation. We name the proposed algorithm as A‐4DEnVar data assimilation method due to the usage of the analytical solution and to distinguish it from conventional data assimilations in which the cost function is minimized through optimal algorithms such as gradient descent, Newton's method or steep descent algorithm.

If the ensemble size is larger than freedom of state, the inverse of B₀₀ can be achieved. In this case, Equation 13 also indicates that the analysis state is actually a solution of a quadric form, with the inverse term of its Hessian matrix term Hess and gradient direction G corresponding to the traditional 4DVar as follows, [Image Omitted. See PDF] [Image Omitted. See PDF]

Traditionally, the inverse of the background error covariance costs huge computational resources, but if the dimension of the model is not so large, the analytical solution can be easily calculated. Otherwise, it is convenient to sample the noise ensemble members from a Gaussian distribution with covariance proportional to the background error covariance. Suppose that we generate members using μB, where μ is small to make sure that the second and higher order term in Taylor expansion can be omitted. It is clear that if the ensemble size is large enough or infinite, B₀₀ equals to μB, otherwise, B₀₀ is an approximation of μB. Then Equation 13 is simplified as [Image Omitted. See PDF]

Using finite ensemble size, it is feasible to achieve the pseudo inverse. How to obtain the inverse of the background error covariance is discussed in the Appendix.

If the model is linear, with sufficient ensemble members, Equation 13 would achieve its optimal solution directly. In nonlinear system, though the structure of adjoint model is unchanged, the value of adjoint model is related to the known state variable ${\mathbf{x}}_0^{\ast }$. It is necessary to introduce an iteration method into the A‐4DEnVar to update the state variable and hence, the value of adjoint model. Because the traditional Newton iteration is sensitive to the initial conditions, the algorithm should be modified to improve its efficiency. To this end, we introduce an iterative linear search process to ensure the convergence of the estimation. The iteration procedure in the proposed A‐4DEnVar has similar effect with the iteration procedure in conventional 4DVar. The analysis increment found from Equation 12 is used as a search direction in the variational procedure. An update step with factor α (with initial value 1.0) is searched to minimize the cost function before next iteration. If the minimum is achieved under some factor, the iteration is finished. Otherwise, a smaller factor, say a half of the previous one, for example, is used to update the analysis state. Theoretically, the factor should be updated carefully but not be roughly halved. But through numerical experiments, we realize that the simple linear search we used is enough to achieve the minimization and reduce the computational cost at the same time. The final form of analysis state is [Image Omitted. See PDF]

The analysis state${\mathbf{x}}_0^a$ will be used to replace ${\mathbf{x}}_0^{\ast }$ in the following iteration. These steps such as generating ensemble members using the new ${\mathbf{x}}_0^{\ast }$ to update the background error covariances and calculating ${\mathbf{x}}_0^a$ in Equation 15 will be iterated until convergence.

Validation With the Low‐Order Model

Twin Experiment Setup

We employ the Lorenz‐63 system to examine the performance of the proposed A‐4DEnVar scheme. The model is widely used in numerical experiments (Hoteit 2008; Lei et al. 2012; Miller 1994; Van Leeuwen, 2009, to name just a few) because of its computational simplicity and strong nonlinear interactions among variables. The Lorenz‐63 model consists of three equations: [Image Omitted. See PDF]

The model is set with the parameter values of (σ,r,b) = (10,28,8/3) (Lorenz, 1963). By using a fourth‐order Runge‐Kutta numerical scheme with time step dt = 0.01, the truth trajectories in different experiments are obtained by integrating the Lorenz‐63 model over the assimilation window. As demonstrated in the study of Goodliff et al. (2015) the performance of the data assimilation schemes does not depend on the true initial conditions, and so the initial conditions (−3.12346395; −3.12529803; 20.69823159) are used here, and are the same as in Goodliff et al. (2015).

In all experiments, observations are produced by adding a Gaussian white noise (the mean value is set to be zero and variance is set to be one) to the true variable values at the corresponding observational intervals. Different observational frequencies are also adopted. These observational frequencies simulate the features of the real climate observing system.

The analysis is based on the verification of three experiments: (1) model control run (CTL) with no observational constraint; (2) data assimilation using conventional 4DVar, in which the background state variables' errors are set to be a Gaussian noise with estimated background error covariance; and (3) data assimilation using the proposed A‐4DEnVar with the same initial background error covariance as used in conventional 4DVar and the ensemble members are generated by adding Gaussian noise with covariance μB, as mentioned in Section 2, to the initial background state variables. Empirically, for a 120 time steps data assimilation window, it is enough to set the parameter μ as 0.0001, but for a very long window, say more than 500 time steps for example, the parameter is set to be 10⁻⁸ in following experiments. Besides, initial conditions in the three experiments are the same.

To rigorously compare to other studies, e.g., Kalnay et al. (2007) and Goodliff et al. (2015) and consider the importance of background error covariance, we conduct the experiments with a background error covariance which is the same as used in Goodliff et al. (2015). Originating from the value used in our study, the true state x^t is generated. Starting from an arbitrary background error covariance, the 3DVar is conducted over a 5,000 time‐step cycle including many assimilation windows. The length of every assimilation window is set to be the same as the observation period, that is, only one observation is at the end of the assimilation window. The observations are the same as used in our study. After the 5,000 time‐step cycle, the background error covariance B is estimated using the forecast states x^f and the truth states x^t, i.e., [Image Omitted. See PDF]from t = 500 to t = 5000. The estimation process is done for more than 10 iterations until the estimation converges.

Figure 1 shows the background solution and 50 ensemble member solutions used in the preliminary test below where the window length is 500 time steps. It is clear that in the nonlinear model, all the initial conditions of ensemble members are nearby the initial background, but as these members evolve over time, the background solution is not the ensemble mean or the solution evolved from the mean of ensemble initial conditions anymore. Thus, as we emphasized, to estimate the adjoint linear model reasonably, ensemble members should be updated near the model state evolved from the guessed initial condition after each iteration.

1 Figure. The 50 ensemble member solutions, truth solution, background solution and solution evolved from the average of ensemble initial conditions of the Lorenz‐63 model over the 500 time steps.

We adopt the widely used root‐mean‐square error (RMSE) as a criterion to evaluate the performance of the conventional 4DVar and the A‐4DEnVar schemes. Here, [Image Omitted. See PDF]where T is the window length, x^analysis and x^truth represent the analysis states and the truth states for variable x, respectively.

Validation Results with the Lorenz‐63 Model

In experiments below, all state variables are observed. With 50 ensemble members, the background error covariance is calculated by Equation A2. These initial ensemble members are integrated to the end of the assimilation window. After temporal cross error covariances are calculated, minimization of cost function can be achieved by calculating the inverse term of Equation A3 and solving Equation 15.

In the preliminary test, the experiment is done for 10 cycles. Each single window length is set to be 500 time steps, which roughly refers to the transition timescale from one attractor to another (Yu et al., 2018), and the observation period is 25 time steps over the assimilation window. The time window and observation period are long enough if the model is an analogy for Rossby wave components because one dimensionless time unit (1 time unit is equivalent to 100 time steps in the experiments) in the Lorenz‐63 model represents several days in atmosphere motions (Ahrens, 1999). Figure 2a shows the results of the assimilation experiments (the first 5 cycles). The analysis states are very close to or even cover the true states under the fact that trajectories are transited from one attractor to the other. For all three variables, it is shown in Figure 2b that the differences between the RMSEs obtained by A‐4DEnVar and conventional 4DVar are all less than 0.06. So the proposed A‐4DEnVar performs as well as the conventional 4DVar scheme once the estimation converges. The RMSEs of the analysis states are less than the standard deviation of observations (that is unit), which demonstrates effectiveness and feasibility for the A‐4DEnVar. In Figure 2c, it is seen that the RMSEs (for the first window) of every variable decrease fast and within 20 iterations, and the A‐4DEnVar obtains the best estimations.

2 Figure. The performances for the proposed analytical four‐dimensional ensemble‐variational (A‐4DEnVar) and conventional four‐dimensional variational (4DVar) over 10 cycles in which each assimilation window is 500 time steps long and with observation every 25 time steps. (a) The solutions of the truth, background, data assimilation results from A‐4DEnVar and conventional 4DVar schemes. (b) The RMSEs for data assimilation schemes. (c) The RMSEs of three variables as a function of iterations when dealing with the first window.

To further evaluate the A‐4DEnVar, we set several other experiments. The data assimilation window lengths are varied from 24 to 720 time steps and observation periods are set to be 6, 12, and 24 time steps for each window. Every experiment is repeated 1,000 times with the same initial state variables and background but different observations. We calculate RMSEs to evaluate the statistical performance of the A‐4DEnVar. As is shown in Figure 3, the RMSEs for all experiments are less than the standard deviation of observations (that is unit, too), so the proposed scheme achieves good performance for each setting. Basically, for a given window length, the RMSE decreases with a smaller observational period. While for a given observational period, a longer data assimilation window would obtain a smaller RMSE due to its richer observational constraints. However, because of the increase of nonlinearity, the performance of A‐4DEnVar decreases when the window length is larger than 480 time steps with observations every 24 time steps. For shorter observation period, i.e., 6 or 12 time steps, the smallest RMSEs are at 480 time steps if the window length is less than 600 time steps. With the growth of nonlinearity, the RMSEs increase after 720 time steps (not shown).

3 Figure. RMSEs of (a) variable x, (b) variable y, and (c) variable z that are calculated from A‐4DEnVar with different window lengths and observation periods.

In operational systems, the dimensions of state variables are often huge. Data assimilation systems are burdened by integrating ensemble members. We evaluate the influence of ensemble size on the A‐4DEnVar scheme. Using the same window length and observation frequency as in the preliminary test, eight ensemble size experiments (from 5 to 1,000 members) are carried out. As the ensemble size increases, the same ensemble members from the previous ensemble size are reused and the extra members are freshly created. The RMSEs are shown in Figure 4. For all experiments, it is also clear that the proposed scheme can achieve similar RMSEs that are less than standard deviation of observation errors (here is set to be 1.0 as before) and background error. Though some experiments with larger ensemble size achieve higher RMSEs, we demonstrate that it is mainly because of the uncertainty induced by random perturbations when generating the ensemble members and 1,000 times iteration are used in all situations to save the computational resources. But considering the RMSE changes little with different ensemble sizes, we demonstrate that the A‐4DEnVar is not sensitive to the ensemble size.

4 Figure. RMSEs of A‐4DEnVar for different ensemble size. The window length and observation period are 500 and 25 time steps, respectively.

Conclusions

Hybrid data assimilation schemes can combine the advantages of conventional 4DVar schemes and ensemble based schemes. The basic motivation of 4DEnVar schemes are estimating the flow‐dependent background error covariance from ensemble members. Other major advantage is that since the adjoint model has been removed, 4DEnVar methods avoid the heavy cost for manual coding or system maintenance. The data assimilation part is separated from the numerical model system, so that it can be easily transplanted to other systems.

Using temporal cross covariances instead of the tangent linear model or the adjoint model, we derived the optimal solution of 4DVar. The derivation provides new insights into the Taylor expansion in ensemble data assimilation schemes. For the temporal covariance used in A‐4DEnVar, the perturbations are the differences between ensemble member states and the state evolved from initial conditions after each iteration. For linear or weakly nonlinear models, if the ensemble size is finite, our solution is equivalent to 4DEnVar proposed by C. Liu et al. (2008) but without any iterations. Like the conventional 4DVar scheme, the A‐4DEnVar can also handle nonlinear problems with a long time window by introducing an optimal step factor. The preliminary evaluation experiments based on the Lorenz‐63 chaotic model indicated a similar performance between conventional 4DVar and the proposed A‐4DEnVar.

With the highly nonlinear Lorenz three variables model, basic issues of the A‐4DEnVar scheme can be clearly examined. However, the huge size of variable states in operational systems requires a feasible technique to alleviate the so‐called sample error problem. Localization methods such as Schur product or localization with Gaussian function (Whitaker & Hamill, 2002) have been proposed to solve the problem in conventional 4DEnVar. As C. Liu et al. (2009) and Amezcua (2017) demonstrated, the empirical orthogonal function decomposition is a more efficient choice. To further explore our method, implementation details with a two dimensional shallow water model and a real atmospheric or oceanic model will be discussed in future works.

Acknowledgments

The authors declare no conflicts of interest. This research is cosponsored by grants from the National Key Research and Development Program (2018YFC1406206, 2016YFC1401800, 2017YFC1404103) and the National Natural Science Foundation (41876014, 11801402) of China.

Data Availability Statement

Data were not used, nor created for this research.

AppendixEnsemble Based Estimation of Analysis Field

In the analytical solution, the optimal initial states are linear combinations of differences from background states and innovations from observations. To achieve the optimal estimation, one problem must be solved, that is, the calculation of the inverse of background error covariance. First of all, because of the large dimensions, it is hard to seek the inverse of background error covariance. However, in operational systems, with about 20 to 50 ensemble members and proper localization technology, the background error covariance can be approximated very well. Especially, the Lorenz model only consists of three variables, so it is enough to use 20 to 50 members to well estimate background error covariance. The fact that the dimension of the ensemble member space is small makes it possible to propose a simple ensemble‐based method to calculate the analysis field.

Assuming there are N ensemble members, the background error is estimated by [Image Omitted. See PDF]where x is a state vector, and $\stackrel{\sim }{\mathbf{X}}$ is a matrix with each column representing an ensemble member perturbation.

If there are infinite members or ensemble size is large enough, it is clear that ${\mathbf{B}}_{00}=\frac{1}{N}{\stackrel{\sim }{\mathbf{X}}}_0{\stackrel{\sim }{\mathbf{X}}}_0^T=\mu \mathbf{B}$. But with finite members, the B₀₀ is just an approximation of μB, so it should be recalculated using ensembles. The covariance of ensemble members is therefore approximated by [Image Omitted. See PDF]

If the dimension of the background error covariance matrix is as small as the Lorenz‐63 model, the inverse term in analysis field can be directly solved. [Image Omitted. See PDF]

While if the dimension is as large as real‐world numerical ocean models have, the dimension of background error covariance matrix is very large and ensemble size is practically not enough, then the matrix is usually singular. So that the pseudo‐inverse but not the inverse is used in following derivation. The Jacobian decomposition of ${\stackrel{\sim }{\mathbf{X}}}_0^T{\stackrel{\sim }{\mathbf{X}}}_0$ is [Image Omitted. See PDF]where V is an orthogonal matrix and Λ is a diagonal matrix with its nonzero elements the eigenvalue of matrix ${\stackrel{\sim }{\mathbf{X}}}_0^T{\stackrel{\sim }{\mathbf{X}}}_0$. It can be seen that ${\stackrel{\sim }{\mathbf{X}}}_0\mathbf{V}$ and Λ are the eigenvector and the eigenvalues matrices of ${\stackrel{\sim }{\mathbf{X}}}_0{\stackrel{\sim }{\mathbf{X}}}_0^T$, respectively. It is easy to verify that the inverse of background error covariance is estimated through [Image Omitted. See PDF]

The generalized inverse of ${\stackrel{\sim }{\mathbf{X}}}_0$ is [Image Omitted. See PDF]

Furthermore, the summation term containing the temporal covariances can be decomposed as [Image Omitted. See PDF]

If we denote $\mathrm{\Delta }\equiv {\mathbf{V}}^T\left({\displaystyle \sum _i\frac{1}{N}{\stackrel{\sim }{\mathbf{X}}}_i^T{\mathbf{H}}_i^T{\mathbf{R}}_i^{-1}{\mathbf{H}}_i{\stackrel{\sim }{\mathbf{X}}}_i}\right)\mathbf{V}$, and using the generalized inverse of ${\stackrel{\sim }{\mathbf{X}}}_0$ in A6, the inverse term in analytical solution is [Image Omitted. See PDF]where I is the identity matrix. From this equation we can see that only ensemble‐size matrix inversion is needed. Other terms in solution can be calculated as in conventional schemes. So, the analytical solution can be obtained.