Asynchronous data assimilation with the EnKF in presence of additive model error

Pavel Sakov; Marc Bocquet

doi:10.1080/16000870.2017.1414545

Introduction

1.

The standard sequential data assimilation (DA) framework assumes that the time is discrete, and that observations assimilated at each time step are made at the analysis time. In practice, the time intervals between consecutive analyses can be substantial, and model states within the observation window can be quite different. It is often important in such cases for a DA method to take into account, in one way or another, the timing of observations. The method is then referred to as asynchronous, in contrast to its standard sequential formulation, which is referred to as synchronous.

The ensemble Kalman filter (EnKF, Evensen, 1994) is based on the Kalman filter (KF), which is a sequential method; however, in the linear perfect model case, the EnKF can be easily adapted to accommodate the asynchronous DA (ADA). Consider the associated minimisation problem:

(1a)

\begin{matrix} x_{0}^{⋆} = arg min J (x_{0}), \end{matrix}

(1b)

\begin{matrix} J (x_{0}) = ‖ x_{0} - x_{0}^{a} ‖_{{(P_{0}^{a})}^{- 1}}^{2} + \sum_{i = 1}^{k} {‖ y_{i} - H_{i} (x_{i}) ‖}_{R_{i}^{- 1}}^{2}, \end{matrix}

(1c)

\begin{matrix} x_{i} = M_{i - 1 \to i} (x_{i - 1}), i = 1, \dots, k . \end{matrix}

Here the lower index i refers to time

t_{i}

x_{i}

is the model state,

x_{0}^{a}

– the initial state estimate,

P_{0}^{a}

– the initial state error covariance estimate,^¹

y_{i}

– observations,

H_{i}

and

R_{i}

– the corresponding observation operator and observation error covariance matrix,

M_{i \to j}

is the model resolvent that maps the state from time

t_{i}

t_{j}

; the norm notation

{‖ x ‖}_{B}^{2} \equiv x^{T} B x

is used.

In order to simplify the notations, we define $b l k d i a g (A_{1}, \dots, A_{k})$ as the block diagonal matrix (or operator) made of the square matrices (or square operators) ${A_{i}\}}_{i = 1, \dots, k}$ . We also define $v e c (A_{1}, \dots, A_{k})$ as the matrix (or operator) that stacks the matrices (or operators) ${A_{i}\}}_{i = 1, \dots, k}$ that are assumed to have the same column number. With these notations, the cost function (1b) can be written as follows:

(2)

\begin{matrix} J (x_{0}) = ‖ x_{0} - x_{0}^{a} ‖_{{(P_{0}^{a})}^{- 1}}^{2} + {‖ y - H \circ M (x_{0}) ‖}_{R^{- 1}}^{2}, \end{matrix}

where

(3)

\begin{matrix} y \equiv v e c (y_{1}, \dots, y_{k}), \end{matrix}

(4)

\begin{matrix} R \equiv b l k d i a g (R_{1}, \dots, R_{k}) \end{matrix}

are the concatenated observation vector and corresponding observation error covariance matrix, respectively, and the forward map $H \circ M (x_{0})$ relates the initial state $x_{0}$ to the observations:

(5)

\begin{matrix} H \circ M (x_{0}) \equiv v e c {H_{i} \circ M_{0 \to i} (x_{0})\}}_{i = 1, \dots, k}) . \end{matrix}

Apart from the forward map (5) that yields forecast observations at times $t_{1}, \dots, t_{k}$ , the cost function (2) has the same form as that for a single DA cycle with synchronous observations. Consequently, in the linear case, when the tangent linear forward operator $H \circ M$ (the sensitivity of observations to the initial state) does not depend on $x_{0}$ , the cost function (2) can be minimised using existing EnKF solutions, provided that the standard (synchronous) observation functions are replaced by the generic forward function (5):

(6)

\begin{matrix} H (E^{f}) \leftarrow v e c {H_{i} \circ M_{0 \to i} (E_{0}^{a})}_{i = 1, \dots, k}), \end{matrix}

where $E$ is the EnKF ensemble; upper “f” refers to forecast; $H (E^{f})$ denotes the forecast ensemble observations and $E_{0}^{a}$ is the initial ensemble. In practice, this means that to assimilate asynchronously with the EnKF one simply needs to calculate ensemble observations using forecast ensemble at observation time (Sakov et al., 2010). There are no specific restrictions on $R$ in (2) – in particular, it does not have to be block-diagonal – therefore in theory observation errors can be correlated in time within a DA cycle.

Before formulated in Sakov et al. (2010), the recipe (6) has been used in a scheme called 4D-LETKF (Hunt et al., 2007), where it was empirically derived from a method called 4D-EnKF (Hunt et al., 2004).

In practical terms, the ADA makes it possible to avoid storing the full ensemble at each observation time, storing only ensemble observations instead. Because in large-scale geophysical systems the model size is typically much larger than the number of observations, this yields a significant economy of resources.

Formulation () of the minimisation problem implies assimilation time $t_{0}$ ; however, one interesting feature of the EnKF in the linear perfect model case is the time invariance of the optimising ensemble transform. If the update is conducted at time $t_{i}$ by applying transform $X$ to a forecast ensemble $E_{i}^{f}$ ,

(7)

\begin{matrix} E_{i}^{a} = E_{i}^{f} X, \end{matrix}

then, to the same effect, it can be applied at other time $t_{j}$ :

(8)

\begin{matrix} E_{j}^{a} = M_{i \to j} E_{i}^{a} = M_{i \to j} (E_{i}^{f} X) = (M_{i \to j} E_{i}^{f}) X = E_{j}^{f} X . \end{matrix}

This time invariance of the ensemble transform has been used in the ensemble Kalman smoother (Evensen, 2003) and made it possible to interpret the ensemble transform as a linear combination of the ensemble trajectories in the 4D-EnKF (Hunt et al., 2004).

The above framework for the ADA with the EnKF is based on two important assumptions: linearity and perfect model. While the KF is a linear method, it can be (and indeed is) successfully applied to weakly nonlinear cases, when the Jacobians $M (x) = \nabla_{x} M (x)$ and $H (x) = \nabla_{x} H (x)$ do not change much through the DA cycle and over the characteristic uncertainty range of the state. The same principle motivates the applicability of the approach (6) to the ADA with the EnKF. This is why we kept the nonlinear notation for the forward operator.^² In strongly nonlinear cases, the method is, generally, no longer applicable, and iterative minimisation algorithms (e.g. Gu and Oliver, 2007; Bocquet and Sakov, 2014) must be used.

The situation in the case with the model error is somewhat similar. Various aspects of additive deterministic (possibly time-correlated) model error and how formulation () should be generalised have been considered in Lorenc (2003), Trémolet (2006) and Carrassi and Vannitsem (2010). However, because (1c) assumes a lossless transfer of information between DA system states in time by the model, it has not been clear whether the concept of the ADA can still be useful in cases with substantial model error. Recently, an iterative method for strongly nonlinear systems with additive stochastic model error called IEnKF-Q has been proposed (Sakov et al., in press) that successfully employs combining propagated ensemble anomalies and model noise ensemble anomalies into a single ensemble of anomalies to obtain an algebraically simple solution for the Gauss–Newton minimisation. The authors suggested that a similar approach can possibly be used to generalise the concept of the ADA with the EnKF in the case of additive model error. This study proposes such a generalisation.

The minimisation problem

2.

Below we mainly follow the formulation and approach to the solution of the minimisation problem from Sakov et al. (in press), with the aim to put up an analogue of the recipe (6) that would permit a simultaneous assimilation of observations made at different times in the case of additive model error without loss of optimality. Sakov et al. (in press) considered a single DA cycle starting at time $t_{0}$ with observations at $t_{1}$ , while we consider a cycle with an arbitrary number of model steps:

(9a)

\begin{matrix} {x_{0}^{⋆}, \dots, x_{k}^{⋆}} = arg min_{{x_{0}, \dots, x_{k}}} J (x_{0}, \dots, x_{k}), \end{matrix}

(9b)

\begin{matrix} J (x_{0}, \dots, x_{k}) = {x_{0} - x_{0}^{a}∥}_{{(P_{0}^{a})}^{- 1}}^{2} + \sum_{i = 1}^{k} {y_{i} - H_{i} (x_{i})∥}_{R_{i}^{- 1}}^{2} \\ + \sum_{i = 1}^{k} {x_{i} - M_{i - 1 \to i} (x_{i - 1})∥}_{Q_{i}^{- 1}}^{2} . \end{matrix}

Here

x_{0}^{a}

is the initial state (the analysis from the previous cycle),

P_{0}^{a}

– the initial state error covariance,

Q_{i}

– the model error covariance at

t_{i}

, and upper star refers to the solution of the minimisation problem.

From a statistical perspective, (9b) relates to the joint probability density function on state and observation vectors since $J (x_{0}, \dots, x_{k})$ coincides with $- 2 ln p (x_{0}, \dots, x_{k}, y_{1}, \dots, y_{k})$ up to a constant. This implies that the model noise, as random variable, is assumed Gaussian, unbiased, of covariance matrix $Q_{i}$ at $t_{i}$ and uncorrelated in time (see Carrassi and Vannitsem, 2010, for a generalisation).

Let the covariances $P_{0}^{a}$ and $Q_{i}$ be factorised by the corresponding ensemble anomalies $A_{0}^{a}$ and $A_{i}^{q}$ :

(10a)

\begin{matrix} A_{0}^{a} {(A_{0}^{a})}^{T} = P_{0}^{a}, A_{0}^{a} 1 = 0, \end{matrix}

(10b)

\begin{matrix} A_{i}^{q} {(A_{i}^{q})}^{T} = Q_{i}, A_{i}^{q} 1 = 0, i = 1, \dots, k, \end{matrix}

where

1

is a vector of entries 1. In many applications, the representation of model errors through the anomalies

A_{i}^{q}

is more natural and convenient than through the

Q_{i}

. A discussion on this specific topic can be found in Section 4 of Sakov et al. (in press). Let us seek the solution of the minimisation problem () in the ensemble space of these anomalies:

(11a)

\begin{matrix} x_{0} = x_{0}^{a} + A_{0}^{a} w_{0}, \end{matrix}

(11b)

\begin{matrix} x_{i} = M_{i - 1 \to i} (x_{i - 1}) + A_{i}^{q} w_{i}, i = 1, \dots, k . \end{matrix}

The weights

w_{i}

are similar to the control variable

η_{i}

in the ‘forcing’ formulation of the weak-constraint 4D-Var (see Equation (2) in Trémolet, 2006), although using

A_{i}^{q} w_{i}

instead of

η_{i}

has a convenience of diagonalising

Q_{i}

The problem () becomes equivalent to

(12a)

\begin{matrix} {w_{0}^{⋆}, \dots, w_{k}^{⋆}} = arg min_{{w_{0}, \dots, w_{k}}} \tilde{J} (w_{0}, \dots, w_{k}), \end{matrix}

(12b)

\begin{matrix} \tilde{J} (w_{0}, \dots, w_{k}) = \sum_{i = 0}^{k} w_{i}^{T} w_{i} + \sum_{i = 1}^{k} {‖ y_{i} - H_{i} (x_{i}) ‖}_{R_{i}^{- 1}}^{2} . \end{matrix}

Let us concatenate

w_{0}, \dots, w_{k}

into

w

(13)

\begin{matrix} w \equiv v e c w_{0}, \dots, w_{k}) ; \end{matrix}

then (12b) becomes

(14)

\begin{matrix} \tilde{J} (w) = w^{T} w + {‖ y - H (x) ‖}_{R^{- 1}}^{2}, \end{matrix}

where $y$ and $R$ are given by (3) and (4), and

(15)

\begin{matrix} x \equiv v e c (x_{0}, \dots, x_{k}), \end{matrix}

(16)

\begin{matrix} H (x) \equiv v e c {H_{i} (x_{i})\}}_{i = 1, \dots, k}) . \end{matrix}

The cost function (14) is similar to the formalism developed by Lorenc (2003) or to Equation (11) of Desroziers et al. (2014), but differs in that it is not incremental and can be obtained directly.

In (14), $x$ is a function of $w$ because of (). Let us expand $x$ in power series of $w$ and retain the linear term only:

(17)

\begin{matrix} x = x^{f} + A w + {O (‖ w ‖}^{2}), \end{matrix}

where^³

(18)

\begin{matrix} x^{f} \equiv v e c {M_{0 \to i} (x_{0}^{a})\}}_{i = 0, \dots, k}), \end{matrix}

and

(19a)

\begin{matrix} A \equiv v e c ({\tilde{A}}_{0}, \dots, {\tilde{A}}_{k}), \end{matrix}

(19b)

\begin{matrix} {\tilde{A}}_{i} \equiv [A_{i}, 0], i = 0, \dots, k, \end{matrix}

(19c)

\begin{matrix} A_{i} \equiv \begin{matrix} A_{0}^{a}, & i = 0 \\ [M_{i - 1 \to i} A_{i - 1}, A_{i}^{q}], & i = 1, \dots, k . \end{matrix} \end{matrix}

Here the ensemble anomalies

{\tilde{A}}_{i}

are completed by zeros so that they all have the same number of columns as

A_{k}

. Note that the augmentation of the propagated state error anomalies and model error anomalies has also been used in the reduced rank square root filter by Verlaan and Heemink (1997) in their Equation (28). Substituting (17) into (14) yields a form convenient for Gauss–Newton minimisation:

(20)

\begin{matrix} \tilde{J} (w) = & w^{T} w \\ + {y - H (x^{f}) - Y w + {O (‖ w ‖}^{2})∥}_{R^{- 1}}^{2}, \end{matrix}

where

(21)

\begin{matrix} Y \equiv v e c {H_{i} {\tilde{A}}_{i}}_{i = 1, \dots, k}) . \end{matrix}

In the nonlinear case, (20) can be minimised by Gauss–Newton iterations. This involves using the approximate Hessian obtained by neglecting nonlinear terms $O (‖ w ‖^{2})$ . In the linear case, the cost function becomes quadratic, and the solution is given by:

(22a)

\begin{matrix} x^{⋆} = x^{f} + A w^{⋆}, \end{matrix}

(22b)

\begin{matrix} A^{⋆} = A T, T = D^{- 1 / 2} U, \end{matrix}

(22c)

\begin{matrix} w^{⋆} = D^{- 1} Y^{T} R^{- 1} y - H (x^{f})], \end{matrix}

(22d)

\begin{matrix} D \equiv I + Y^{T} R^{- 1} Y, \end{matrix}

where

{(\cdot)}^{1 / 2}

denotes the unique positive (semi)definite square root of a positive (semi)definite matrix, and

U

is an arbitrary mean preserving orthonormal matrix^⁴:

U U^{T} = I, U 1 = 1

These equations provide the update of the ensemble at all time levels. Similarly to the perfect model case, the same coefficients $w^{⋆}$ and the same ensemble transform $T$ are used over the whole assimilation window. However, due to the augmentation of the propagated ensemble anomalies with model error anomalies in (19c), the updated trajectories of the initial ensemble $E_{0}^{⋆}$ are no longer continuous in time.

Asynchronous DA

3.

Based on the solution of the minimisation problem in the previous section, we can now formulate a framework (rules) for the ADA with the EnKF in presence of additive model error. According to equations (), to update the system state, one needs to calculate the forecast observations $H (x^{f})$ and forecast ensemble observation anomalies $Y$ ; they in turn allow one to calculate the Hessian $2 D$ (approximate Hessian in the nonlinear case), coefficients $w^{⋆}$ and transform $T$ . Therefore, the essence of the ADA with the EnKF can be formulated as follows:

(23)

\begin{matrix} H (x^{f}) = v e c {H_{i} (x_{i}^{f})}_{i = 1, \dots, k}), \end{matrix}

(21)

\begin{matrix} Y = v e c {H_{i} {\tilde{A}}_{i}}_{i = 1, \dots, k}), \end{matrix}

where ${\tilde{A}}_{i}$ is defined by (19b) and (19c).

Subject to the new evolution equation for ensemble anomalies (19c), the rules (21) and (23) can be interpreted in the same way as in the perfect model case: calculate ensemble observations using the forecast ensemble at observation time.

The sequence of the complete analysis cycle in both the linear (upper) and the nonlinear (lower) case. From left to right, the nodes indicate the equation to use as referred to in the paper and the yielded quantities. The second, lower chain in the nonlinear case is needed if the observation operator is nonlinear. Arrows mean that the computation of the target node content requires the content of the origin node.

The sequence of the complete linear analysis cycle is synthesised in Fig. 1 (upper panel).

Discussion

4.

Increasing ensemble size

4.1.

The update equations () have the same form as in the perfect model case; however, the propagation of the ensemble anomalies in (19c) is substantially different. According to (19c), in presence of model error, the forecast ensemble of anomalies is not only transformed with the tangent linear model, but is also augmented with the corresponding ensemble of model error anomalies. This procedure implicitly implements the forecast covariance equation of the KF

(24)

\begin{matrix} P_{i} = M_{i} P_{i - 1} M_{i}^{T} + Q_{i}, \end{matrix}

so that $A_{i} A_{i}^{T} = P_{i}$ , and manifests an effective increase of the ensemble size at every time step:

(25)

\begin{matrix} m_{i} = m_{i - 1} + m_{i}^{q}, i = 1, \dots, k, \end{matrix}

where $m_{i}$ is the size of $A_{i}$ ; $m_{i}^{q}$ is the size of $A_{i}^{q}$ and $m_{0} \equiv m$ is the size of $A_{0}^{a}$ , i.e. the initial ensemble size. The practical choice for the $m_{i}^{q}$ has been discussed in Section 4 of Sakov et al. (in press).

Note that this increase of the ensemble size may not be strictly necessary for implementing (24). For example, if the column space of $A_{i}^{q}$ is a subspace of the column space of $M_{i} A_{i - 1}$ , then adding model error covariance can be done without changing the ensemble size. Yet, the augmentation of the state error and model error ensemble anomalies in (19c) is an essential part of the described minimisation procedure, regardless of whether it is possible to factorise the forecast state error covariance with a smaller ensemble or not. It enables to concatenate the vectors of ensemble weights $w_{i}$ into a single vector $w$ in (13) so that each component $w_{i}$ of this vector acts only on the columns of the augmented ensemble associated with the model error added at time $t_{i}$ (except that $w_{0}$ is indeed associated with the initial state error). In this way, the minimisation uses the same vector of coefficients $w$ over the whole observation window, similarly to the perfect model case.

Nonlinearity and approximation of $M$ and $H$

4.2.

Although the EnKF is a linear method targeting weakly nonlinear systems, it is still desirable to optimise it for strongly nonlinear situations, for two reasons. Firstly, it is quite common for a weakly nonlinear system to become strongly nonlinear during relatively short and infrequent transition periods. Secondly, there are a variety of EnKF-based iterative methods for strongly nonlinear systems that would benefit from such an optimisation.

The EnKF is a derivative-less method: it approximates the Jacobians $M_{i}$ and $H_{i}$ based on application of the forward function $H_{i} \circ M_{i}$ to the ensemble of model states $E_{i - 1}$ . In the linear case, the results of these approximations do not depend on the ensemble spread; in the nonlinear case to achieve correct sampling of nonlinear effects, the magnitude of the ensemble anomalies should be of order of the standard deviation of the state error. For a system with additive model error, this means that the forecast error anomalies should have magnitude similar to that of samples from $N (0, P_{i})$ , where the forecast covariance $P_{i}$ is given by (24).

Such sampling can be associated with the statistical estimator of the covariance,

(26)

\begin{matrix} P = \frac{1}{m_{A} - 1} A A^{T}, \end{matrix}

where $A$ is an ensemble of $m_{A}$ anomalies. Using this factorisation yields the characteristic magnitude of the anomalies of the order of its standard deviation, regardless of the ensemble size.

The minimisation procedure in Section 2 uses a different factorisation,

(27)

\begin{matrix} P = A A^{T}, \end{matrix}

which simplifies the algebraic side. Using this factorisation reduces the characteristic magnitude of anomalies with increasing ensemble size. To obtain anomalies of physically ‘right’ magnitude, one needs to scale them up with a factor of $\sqrt{m_{A} - 1}$ , but only if all samples in the ensemble are drawn from the same pool. If, for example, we introduce a new ensemble

\begin{matrix} B = [A, 0], \end{matrix}

then we should still scale it with $\sqrt{m_{A} - 1}$ rather than with $\sqrt{m_{B} - 1}$ , where $m_{B}$ is the size of $B$ .

Because the magnitudes of the anomalies $M_{i} A_{i - 1}$ and $A_{i}^{q}$ in (19c) cannot be assumed to be similar, it seems technically challenging to put up a generic re-scaling procedure for switching between factorisations (26) and (27) in the course of minimisation. This leaves us the following two main options.

(1) Using factorisation (27) over the whole DA cycle, which means no rescaling before or after application of nonlinear functions $M$ and $H$ . This approach underestimates the ensemble spread in approximations and shifts the minimisation towards the Newton method from the more inherent for the EnKF (and generally more robust) secant method.
(2) Scaling components of $A_{i}$ appended at different time steps according to the size of these components before application of $M$ and $H$ , and scaling back after that. This approach still does not make all anomalies of $A_{i}$ being sampled from $N (0, P_{i})$ , but can be a good practical compromise. It is elaborated below.

First, we define the forecast ensembles

E_{i}

(28)

\begin{matrix} E_{i} = \begin{matrix} E_{0}^{a}, & i = 0, \\ M_{i - 1 \to i} (E_{i - 1}), \\ x_{i}^{f} 1^{T} + \sqrt{m_{i}^{q} - 1} A_{i}^{q}], & i = 1, \dots, k, \end{matrix} \end{matrix}

where $E_{0}^{a}$ is the initial ensemble and

(29)

\begin{matrix} x_{i}^{f} \equiv M_{i - 1 \to i} (E_{i - 1}) 1 / m_{i - 1} . \end{matrix}

We then define the new anomalies ${\tilde{A}}_{i}$ that yield $A$ via (). They are composed of sub-blocks $A_{i}^{j}$ , each of them corresponding to the ensemble appended at time $t_{j}$ :

(30)

\begin{matrix} {\tilde{A}}_{i} \equiv [A_{i}^{0}, \dots, A_{i}^{k}], \end{matrix}

where

(31)

\begin{matrix} A_{i}^{j} \equiv \begin{matrix} (E_{i}^{0} - x_{i}^{f} 1^{T}) / \sqrt{m_{0} - 1}, & j = 0 \\ (E_{i}^{j} - x_{i}^{f} 1^{T}) / \sqrt{m_{j}^{q} - 1}, & 1 \leq j \leq i \\ 0, & i < j \leq k . \end{matrix} \end{matrix}

The submatrix $E_{i}^{0}$ corresponds to the first $m_{0}$ columns of $E_{i}$ ; $E_{i}^{j}$ , for $1 \leq j \leq k$ , corresponds to the submatrix from column $m_{j - 1} + 1$ to column $m_{j}$ of $E_{i}$ .

The forecast observation anomalies (21) are then calculated as

(32)

\begin{matrix} Y = v e c Y_{1}, \dots, Y_{k}), \end{matrix}

where

(33)

\begin{matrix} Y_{i} \equiv [Y_{i}^{0}, \dots, Y_{i}^{k}], \end{matrix}

and

(34)

\begin{matrix} Y_{i}^{j} \equiv \begin{matrix} [H_{i} (E_{i}^{0}) - H_{i} (x_{i}^{f}) 1^{T}] / \sqrt{m_{0} - 1}, & j = 0 \\ [H_{i} (E_{i}^{j}) - H_{i} (x_{i}^{f}) 1^{T}] / \sqrt{m_{j}^{q} - 1}, & 1 \leq j \leq i \\ 0, & i < j \leq k . \end{matrix} \end{matrix}

To summarise, the sequence of the complete nonlinear analysis cycle is synthesised in Fig. 1 (lower panel).

Observations correlated in time

4.3.

The formulation of the minimisation problem () assumes that observations at different times are noncorrelated, so that the merged observation error covariance matrix $R$ in (4) is block-diagonal; however, we can see no problem in removing this restriction, that is assuming that the observations within a DA cycle can be time correlated.

Ensemble resizing

4.4.

One practical implication of the proposed framework for ADA in presence of additive model error is the need for ensemble resizing at the end of a cycle to cut back the increased ensemble size due to (19c). If the augmented ensemble $A_{k}$ has the same rank as $A_{0}$ , then the ensemble size reduction from $m_{k}$ to $m_{0}$ can be performed losslessly (that is, without changing the analysed state error covariance). In large-scale geophysical EnKF systems, the ensemble size is always much smaller than the model size, and realistic implementations of model noise can be expected, generally, to project onto the left null-space of the ensemble. In this case, the resizing will involve removing the principal components with the smallest variance (Verlaan and Heemink, 1997; Raanes et al., 2015; Sakov et al., in press).

Observation time binning

4.5.

The minimisation problem () associates each time step with observations $(y_{i}, R_{i})$ and model error $Q_{i}$ . In practice the observations can be distributed in time continuously. They will need then to be binned by some time intervals, and model error estimates will need to be prescribed for each of these intervals. Because of the overhead associated with augmenting the ensembles of model error anomalies into the ensemble anomalies by (19c) one may have to use a rather coarse time binning for the model error. It is then possible to assume that the forward function $H_{i} \circ M_{i}$ in (34) propagates ensemble to the time of each particular observation, as the overhead associated with the time binning of ensemble observations can be much smaller. This is the same observation time binning strategy as the one proposed by Trémolet (2006) with his Equation (9).

ADA: further extending the KF

4.6.

The EnKF yields a large number of advantages over the KF: scalability, derivative-less formulation, ease of implementing the covariance propagation equation, robustness, and so on. Most of these advantages can be viewed as concerning the technical side of the Kalman filtering. By contrast, the ADA aims at extending the very framework of the sequential DA. While in the perfect model case it is possible to modify the square root formulations of the KF to permit the ADA similarly to that in the EnKF, the technique proposed in this note involves using nonsquare covariance factorisations of varying size within a DA cycle. We see it as yet another step in the evolution of the KF-based methods towards greater universality, well beyond the original scope of the KF.

Summary

5.

Equations (23), (21) and (19c) provide a framework for the ADA with the EnKF in presence of additive model error. Its derivation in the course of solving the minimisation problem warrants its optimality in the linear case.

The most essential new component of the proposed framework is Equation (19c) for propagation of ensemble anomalies. This equation manifests an effective increase of the ensemble size at each time step by the size of the ensemble of model error anomalies. We argue that this increase is necessary for the optimal ADA in presence of additive model error regardless of whether the model error anomalies project onto the column space of the ensemble or not; in practice, it deems it necessary to resize the ensemble at the end of the cycle.

Disclosure statement

[1] Conflicts of interest No potential conflict of interest was reported by the authors.

Notes

[2] The upper label ‘a’ refers to analysis as, later on, these initial state and error covariance matrix will coincide with the outcome of a previous analysis.

[3] Another reason is that in the linear case functions $M$ and $H$ can be affine: $M (x) = m + M x, H (x) = h + H x$ .

[4] In the EnKF instead of explicit propagation of state, it is common to use the average of the propagated ensemble: $M (x) \leftarrow M (E) 1 / m$ , where m is the size of ensemble $E$ .

[5] The presence of $U$ here symbolically accounts for the variety of the deterministic EnKF schemes.

Asynchronous data assimilation with the EnKF in presence of additive model error

Full Article

Introduction

1.

(1a)

(1b)

(1c)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

The minimisation problem

2.

(9a)

(9b)

(10a)

(10b)

(11a)

(11b)

(12a)

(12b)

(13)

(14)

(15)

(16)

(17)

(18)

(19a)

(19b)

(19c)

(20)

(21)

(22a)

(22b)

(22c)

(22d)

Asynchronous DA

3.

(23)

(21)

Figure 1.

Discussion

4.

Increasing ensemble size

4.1.

(24)

(25)

Nonlinearity and approximation of M and H

4.2.

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

(34)

Observations correlated in time

4.3.

Ensemble resizing

4.4.

Observation time binning

4.5.

ADA: further extending the KF

4.6.

Summary

5.

Disclosure statement

Notes

Nonlinearity and approximation of $M$ and $H$