The ensemble Kalman filter combines the above introduced notion of an ensemble of model states to describe spatial and temporal correlations with the well-known Kalman Filter [27]. The pending task of generating an analysis ensemble obeying an obtained analysis correlation matrix can be completed by a square root filter (SRF), originating from the Kalman filter update for the correlation matrix applied to the ensemble representation

with the Kalman gain matrix

and the transformation matrix S given by

Taking the square root of the symmetric matrix U results in

which is the transformation matrix of the update for the centered ensemble in Equation (8). Note, the notation in Equations (2, 6) is the one used by Nakamura and Potthast [3] and differs from the one introduced by Hunt et al. [16] (see Equation 19). However, by multiplying S in Equation (30) with the inverse of W^a defined by Hunt et al, the identity W^a = S can be easily shown.

The update of the mean is obtained along the lines of the classical Kalman filter by

with K given in Equation (28). Comparing Equation (14, 31) leads to

in case of the SRF. The update for the full ensemble using the ensemble Kalman square root filter is therefore given by applying Equations (30, 31) to Equation (11, 15).

Here, we can now confirm the validity of Equation (17). From the definition of Q^b we know that the sum of the rows of Q is zero, such that the sum of the column of s̄=(Qb)TA with any matrix A ∈ ℝ^n×L is zero, and the sum of the columns of I − (Q^b)^TA is one. If I − (Q^b)^T A is symmetric, this means that the vector equal to 1 in each component is an eigenvector of I − (Q^b)^T A with eigenvalue 1. But then it will also be an eigenvector with eigenvalue 1 for each power of I − (Q^b)^TA, such that (17) is satisfied.

Before we investigate ultra-rapid data assimilation based on reduced data, we need to recall how a standard ensemble Kalman square root filter will react when we base our analysis on a reduced set of model variables. Let us study the calculation of the ensemble analysis for the ensemble Kalman filter with reduced data. The basic formula for the ensemble Kalman filter can be expressed as W=S+s̄ with S and s̄ given in Equations (30, 32)

Now, assume we observe y ∈ ℝ^m which depends on some subset x₁, …, x_ñ of the full set of variables x₁, …, x_n only. Given these reduced spaces the operator H will be of the form

In this case, the terms HQ^b ∈ ℝ^m×L and (HQ^b)^T ∈ ℝ^L×m will be a linear combination of the variables 1, …, ñ of the ensemble members. If we are given the variables x₁, …, x_ñ of the ensemble members only, the matrix W will not change. Also, y − Hx^b will depend only on the variables x1b,…,xñb of x^b. The solution z ∈ ℝ^m of

is calculated based on the variables x₁, …, x_ñ of Q^b and x1b,…,xñb of x^b only. We summarize the result of these arguments in the following lemma.

Lemma 2.1. If we have observations y dependent only on the variables x₁, …, x_ñ for ñ ∈ {1, …, n} of the full state x ∈ ℝⁿ of the state space of our dynamical system, the transformation matrix W of the ensemble Kalman square root filter update x^a − x^b depends on these variables of the centered ensemble Q and the mean first guess x̄b only.

A consequence of the above Lemma 0.0.1 is that, if we have reduced observations, the ensemble Kalman square root filter will give us an update matrix W which depends only on the variables under consideration.

But we need to pay attention to the update and propagation step. The update Equation (11) clearly updates all variables, since W ∈ ℝ^L×L, and thus all variables of x are updated by the ensemble Kalman filter. If the model M is based on all variables, in general we expect model propagation to be dependent on all variables as well. In general, an update based on the transformation matrix W will change all variables of the initial state. This means that the first guess of the next assimilation step highly depends on the application of the matrix W to all variables, not only to the variables x₁, …, x_ñ.

Clearly, in general we cannot run the full ensemble Kalman filter on a reduced set of variables, just because you need all prognostic variables to run the numerical model. We will see later, that this limitation does no longer apply when we are in the framework of ultra-rapid data assimilation.

Assume that we are given some ensemble xka,(1),…,xka,(L) of L states of our dynamical system at time t_k ∈ ℝ, which could be an analysis or a first guess from somewhere. Further, we assume that we have applied our model M to calculate forecasts based on xka,(ℓ) at times t_k+1, …, t_N for N > k. The corresponding forecasts are denoted by xk,ξf,(ℓ) for ξ = k + 1, …, N, analogous to (35).

We employ the following matrix notation. The matrix F is the matrix of the full forecast ensemble members in its columns, i.e.,

of forecasts xk,ξf,(ℓ) from t_k to t_ξ. The matrix W^(k) is the matrix of linear ensemble transform coefficients calculated based on the observations y_k at time t_k and the first guess ensemble at time t_k, i.e.,

When the analysis ensemble at time t_k is given by a generic ensemble data assimilation approach, we know that

with the matrix Wj,ℓ(k), j, ℓ = 1, …, L given by Equation (16) with the two quantities s̄ and S being dictated by the specific ensemble data assimilation system (e.g., Equations 30, 32), where the time index k refers to the analysis time t_k for which W_{j, ℓ} is calculated. Also, we note that the background xkb,(ℓ) is given by

Lemma 3.1. Here, we assume that the model M is a linear model M. In this case, the forecast ensemble xk,ξf,(ℓ) at time t_ξ when observations at time k are assimilated by a linear data assimilation method as in Equation (12), the forecast ensemble can be calculated by

Proof. In this case we have

for ℓ = 1, …, L and ξ ∈ {1, …, N}, where we used M_{k − 1, ξ} = M_{k, ξ}M_{k − 1, k}. □

Before we continue with our introduction of ultra-rapid data assimilation, we would like to study the reduced variable case in the above Lemma 0.0.2. Clearly, to apply M_{k, ξ} to a state x^(a) or x^(b), we need to know the full state. If only a part of the state x is available, starting the model is no longer possible. However, the Equation (40) is still valid for each of its components, i.e., if W is known, the variable xk,ξ,if,(ℓ) of xk,ξf,(ℓ) can be calculated from the knowledge of xk-1,ξ,if,(ℓ) for all ℓ = 1, …, L, representing the i-th variable of the state vector of the l-th ensemble member obtained by a forecast from time t_{k − 1} to time t_ξ.

Corollary 3.2 (Reduced Set of Model Variables). If the observation operator H depends on the variables x₁, …, x_ñ of the state x only, then the transformation matrix W^(k) for the assimilation of y_k can be calculated from a) the first guess ensemble data x1(b),…,xñ(b) and b) the observation y_k. For a linear model M, for the variables with index i we have

for i = 1, …, ñ, i.e., the formula (40) is valid and the ensemble forecast based on the analysis with observation y_k can be calculated from the knowledge of the reduced set of variables only.

The consequence of Equation (41) is that for linear models we can calculate the forecast based on the analysis at time t_k by a superposition of the forecast from time t_{k − 1}. The weight matrix Wℓ,j(k) is calculated from the ensemble analysis at time t_k given by the linear ensemble data assimilation scheme. We can also use Equation (41) recursively, which is formulated in the following Theorem.

Theorem 3.3. We assume we are given observations y_j, j = 1, …, k at times t₁, …, t_k. The goal is to calculate the forecasts xk,ξf,(ℓ) at time t_ξ based on the observations from t₁ to t_k and the initial ensemble x0a,(ℓ) at time t₀ with an ensemble data assimilation method as in Equation (11). If the model M is linear, we obtain

for ξ = k+1, …, N.

Proof. For a linear model, the generic step is given by Equation (41). Then, the same equation is applied to xk-1,ξf,(j), which leads to

and by the same step η times to

for η ≤ k assimilation steps. In matrix notation and for η = k this is Equation (43). □

Note that the recursive application of Equation (41) implies that any transformation matrix W⁽ⁱ⁾ is obtained using the observation y_i and the full ensemble F_{i − 1, ξ}.

The results for reduced data are also valid for the core formula (43). We collect the relevant statements into the following corollary. The matrix F_{k, ξ} contains the different state variables in its rows and the columns represent the ensemble under consideration. We employ the notation (_{Fk, ξ)i = 1, …, ñ} for the rows with the variable indices i = 1, …, ñ.

Corollary 3.4 (Reduced Set of Model Variables). If the observation operator H depends on the variables x₁, …, x_ñ of the state x only, then the transformation matrix W^(k) for the assimilation of y_k can be calculated from a) the first guess ensemble data (F_{0, k})i = 1, …, ñ , b) the observation y_k and c) the previous transformation matrices W⁽¹⁾⋯W^{(k − 1)} which depend on the corresponding observations y₁⋯y_{k − 1}. For a linear model M, for the variables with index i we have

i.e., the formula (43) is valid and the ensemble forecast based on the analysis with observation y_k can be calculated from the knowledge of the reduced set of variables only.

Smoothers are schemes which employ information from the future to improve the estimate about some present state. Alternatively, you could say that they use information now to update past states.

When we consider the scenario of ultra-rapid data assimilation, for the interval [t₀, t_N] we are given an ensemble of original states (35) over the full interval. When an observation is arriving at time t_k (ignoring delay usually needed for observation processing and transfer), we can employ the same techniques which are used for updating the analysis and forecast to the past interval [t₀, t_k].

Definition 3.5 (Ultra-Rapid Ensemble Smoother).

Given the original first guess ensemble F_{0, ξ} for ξ = 0, …, N on the time interval [t₀, t_N] we define the ensemble analysis given the data y₁, …, y_k by

This analysis ensemble is defined for the full time interval.

In general, a convergence analysis of an ensemble Kalman smoother and its comparison to a four-dimensional variational data assimilation (4D-VAR) scheme over the time window [t₀, t_n] can be found in Theorem 5.4.7 of Nakamura and Potthast [3]. For linear models and observation operators, the full Kalman smoother and 4D-VAR are equivalent.

Clearly, if we replace the full model M by the ensemble, this equivalence is no longer true. Also, if the numerical model M used to calculate the ensemble is different from the true model M_true, the temporal correlations, which are implicitly used when we employ the analysis matrix W^(k) to update the ensemble in the past or in the future, may not be correct with respect to the true ensemble correlations. In this case, the information y_k in the future of t₀ may not improve the state estimate at time t₀, but lead to additional errors in this state estimate. We will demonstrate this phenomenon in our numerical examples in section 4.

Here, we start our study with the Lorenz 63 model Lorenz [20]. It is a very well-known chaotic ODE system with three unknowns, compare for example Nakamura and Potthast [3].

The Lorenz 1963 model is a system of three non-linear ordinary differential equations

with constants σ, ρ, β known as Prandtl number, the Rayleigh number and a non-dimensional wave number. Here, for the constants we take the classical values σ = 10, β = 8/3 and ρ = 28. The implementation of the system is usually carried out by a higher-order integration scheme such as 4th-order Runge-Kutta, which we have employed for our numerical testing. The setup for our case study is shown in Figure 1A with 8 cycles for better visibility and Figure 1B with 40 cycles for studying the error evolution.

Here, we want to test the feasibility of ultra rapid data assimilation. The original curve is shown in black in Figure 1. The measurements are calculated by adding a random Gaussian error to this curve at the measurement times t₁, t₂, …, t_k with Δ_t = t_i+1 − t_i = 0.1 (without units). For the original, we have used the above ODE system with σ = 10 to generate the truth. To test data assimilation we have employed a modified system where σ = 12 was chosen. The mean of the original first guess ensemble for the full time period under consideration is shown in Figure 1 as a blue curve, with blue dots as the original first guess.

We have now followed two tracks. First, we have implemented an ensemble Kalman square root filter. We start with a first guess ensemble, which is generated at time t₀ by adding random Gaussian errors to the starting point of the original curve. Then we assimilate the observations (the black dots) using the Ensemble Kalman square root filter.

Second, the ultra-rapid data assimilation and forecasting cycle has been implemented. The ultra rapid data assimilation has been set up by first calculating the full first guess ensemble for the whole time interval under consideration. Then, a modified ensemble is calculated step by step following (43). We study N time steps (showing results for N = 8 and N = 40). In more detail, we have calculated the transformation matrix W^(k) based on the observations y_k at time t_k, k = 1, …, N and the transformed first guess ensemble xk-1,ξb,(ℓ). Here, ξ is the time index of the ensemble, i.e., ξ = 1, …, N. We carry out the assimilation for all time steps, changing the ensemble in the past as well as in the full future over the time interval under consideration.

The result of N time steps is shown in Figure 2. First, the example with N = 8 time steps is shown in Figure 2A,B, the first guess errors for N = 25 and N = 40 time steps in Figure 2C,D. Here, initially the ultra-rapid update is quite good, approximating well the full ensemble Kalman square root filter over 10 or 15 assimilation steps. Then, when the first guess ensemble and the true trajectory diverge further, the assimilation looses track and we obtain very large errors over time, as can be seen by the peak of the pink curve in Figure 2D at about t_k with k = 34.

Here, we also investigate the ultra rapid data assimilation tool as a smoother. We calculate the analysis ensemble Fk(a) defined in (47) based on the original first guess ensemble F₀.

In Figure 3 we study the filter and smoother results for a case with N = 32. For the latter we update both the future and the past. Errors of this with respect to the true curve are displayed in Figure 3C,D . Here, we need to note that we simulate a realistic setup in the sense that the true model M_true is different from the model M used to calculate the first guess ensemble. That has severe consequences for the convergence of the smoother. With the errors in the model, we obtain errors in the first guess ensemble and with this errors in the correlations and covariances which are exploited by the ultra-rapid Kalman filter forecast and the update of the past in the ensemble Kalman smoother.

Studying Figure 3C we see that the error is smallest on the diagonal, i.e., for the analysis and short-term forecast based on the ultra-rapid ensemble Kalman analysis. The errors for the analysis or forecast increase with distance to the current point in time. We expect that the errors are large for larger lead times. But in general we do not expect that the errors in the past, i.e., at the beginning of the interval [t₀, t_N] increase when we assimilate more and more data. When the ensemble reflects the correct correlations between the future and the past, the error should decrease. However, with a numerical model which is different from the true model, we also inherit errors into the temporal correlations. As a consequence we observe that the error at t₁ increases when we assimilate further data y_k for k in the second part of the interval [t₀, t_N].

In Figure 4 we evaluate the performance of URDA in a statistical manner by using different initialisations for the applied random number generator, which affects the observations drawn from a Gaussian distribution as well as the construction of the ensemble, and use different values of the starting point x₀ = x(t₀), which is used to obtain the truth as well as the ensemble. We evaluate differences of the corresponding mean from the truth and take appropriate ratios. In Figure 4A the mean-error of URDA is divided by the mean error of the free forecast (no data assimilation, also abbreviated by no-DA) for L = 5 ensemble member and N = 25 time steps on the trajectory. A clear positive impact is visible with only very few cases where the free forecast is better than URDA. Figure 4B shows the mean-error of the SRF divided by the one of URDA. As expected, the evaluation shows that in many cases the full SRF performs better compared to URDA. However, comparing with the results shown in Figure 3C this is what we expect due to the deviations after about N = 8 time steps. To test this, we show the result for the first 8 time steps of the run with N = 8 in Figure 4C and observe, that for a smaller N these compete indeed much better with the SRF. Figure 4D shows the result of the mean-error of URDA divided by the initial forecast with no data assimilation for L = 25 ensemble member and N = 25 time steps. We observe, that the improvement of URDA with L = 25 compared to L = 5 is not significant. This is no surprise since we deal with three prognostic variables where an ensemble of L = 5 is already sufficient to describe the relevant spread.

At the end of this section we highlight the impact of the time step in the model, which translates to the time the forecast from one point on the trajectory is performed. Note, this does not affect the performance of the Runge-Kutta-Scheme where the time step of the integration is kept fixed. We evaluate the ratio of the deviations from the mean error from the SRF to URDA. In Figure 5 we show results for different sizes of the time step dt. Again we used N_stat = 250 and the total number of time steps N = 25 with the number of ensemble members L = 5. We observe, that for ¹/₄ of the standard time step size dt = 0.100 the SRF and URDA perform almost equally. With increasing time step size we find more cases where the SRF outperforms URDA, which is still moderate for the standard time step size. For three times this step size we see a clear benefit for the SRF. Note, since we keep N = 25 fixed, the total length of the trajectory differs for the different time step sizes dt.

In the second part of our numerical study we would like to understand how ultra-rapid ensemble data assimilation can be applied to the case where only a reduced set of variables is passed down from the standard ensemble data assimilation framework.

In the framework of the Lorenz model, we have carried out a study the use of the observation operators

and study assimilation of observations of either the full state x, the first two variables of the state or the first variable of the state only.

We note that HQ will employ the corresponding selection of variables depending on the cases H₁, H₂ or H₃, i.e., H₁Q can be calculated from the knowledge of the first variable x₁ of x=(x1,x2,x3)T only. Similarly, H₂Q can be calculated based on the knowledge of (x₁, x₂) of x=(x1,x2,x3)T.

Here, we focus on the results for the use of H₂ in Figure 6. The effects are similar to the three-dimensional version. Figure 6A displays the first 8 steps, and we see that the SRF analysis and the URDA analysis are very close to each other. The error is shown in Figure 6B, here only for the two variables under consideration. Figure 6C,D display 30 assimilation setups. After 20 and 25 steps we observe first cases where URDA is worse than no-DA. In all other cases it is a big increase from the no-DA case and its quality becomes close to the quality of the full square-root filter with subsequent forecast.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.