# Quasi-static ensemble variational data assimilation: a theoretical and numerical study with the iterative ensemble Kalman smoother

Anthony Fillion, Marc Bocquet, Serge Gratton

► To cite this version:

Anthony Fillion, Marc Bocquet, Serge Gratton. Quasi-static ensemble variational data assimilation: a theoretical and numerical study with the iterative ensemble Kalman smoother. Nonlinear Processes in Geophysics, 2018, 25, pp.315-334. 10.5194/npg-25-315-2018. hal-02147982

HAL Id: hal-02147982

https://hal.science/hal-02147982

Submitted on 5 Jun 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OATAO

Open Archive Toulouse Archive Ouverte

Open Archive Toulouse Archive Ouverte

OATAO is an open access repository that collects the work of Toulouse researchers and makes it freely available over the web where possible

This is a publisher's version published in:
http://oatao.univ-toulouse.fr/22606

Official URL

DOI: https://doi.org/10.5194/npg-25-315-2018

**To cite this version:** Fillion, Anthony and Bocquet, Marc and Gratton, Serge *Quasi-static ensemble variational data assimilation: a theoretical and numerical study with the iterative ensemble Kalman smoother.* (2018) Nonlinear Processes in Geophysics, 25. 315-334. ISSN 1023-5809
<figure><img src="image_3.png" /></figure>

# Quasi-static ensemble variational data assimilation: a theoretical and numerical study with the iterative ensemble Kalman smoother

Anthony Fillion<fnref n="1" />,<fnref n="2" />, Marc Bocquet<fnref n="1" />, and Serge Gratton<fnref n="3" />

**Correspondence:** Anthony Fillion (anthony.fillion@enpc.fr)

Received: 30 October 2017 – Discussion started: 16 November 2017

Revised: 14 March 2018 – Accepted: 15 March 2018 – Published: 27 April 2018

**Abstract.** The analysis in nonlinear variational data assimilation is the solution of a non-quadratic minimization. Thus, the analysis efficiency relies on its ability to locate a global minimum of the cost function. If this minimization uses a Gauss-Newton (GN) method, it is critical for the starting point to be in the attraction basin of a global minimum. Otherwise the method may converge to a *local* extremum, which degrades the analysis. With chaotic models, the number of local extrema often increases with the temporal extent of the data assimilation window, making the former condition harder to satisfy. This is unfortunate because the assimilation performance also increases with this temporal extent. However, a quasi-static (QS) minimization may overcome these local extrema. It accomplishes this by gradually injecting the observations in the cost function. This method was introduced by Pires et al. (1996) in a 4D-Var context.

We generalize this approach to four-dimensional strong-constraint nonlinear ensemble variational (EnVar) methods, which are based on both a nonlinear variational analysis and the propagation of dynamical error statistics via an ensemble. This forces one to consider the cost function minimizations in the broader context of cycled data assimilation algorithms. We adapt this QS approach to the iterative ensemble Kalman smoother (IEnKS), an exemplar of nonlinear deterministic four-dimensional EnVar methods. Using low-order models, we quantify the positive impact of the QS approach on the IEnKS, especially for long data assimilation windows. We also examine the computational cost of QS implementations and suggest cheaper algorithms.

# 1 Introduction

## 1.1 Context

Data assimilation (DA) aims at gathering knowledge about the state of a system from acquired observations. In the Bayesian framework, this knowledge is represented by the posterior probability density function (pdf) of the system state given the observations. A specificity of sequential DA is that observations are not directly available; they become available as time goes by. Thus, the posterior pdf should be regularly updated.

In order to do so, one usually proceeds in two steps: the analysis and the propagation (or forecast). During the analysis step, a background pdf is used as a prior together with the observation likelihood to construct the (often approximate) posterior pdf, following Bayes' theorem. During the propagation step, this posterior pdf is propagated in time with the model to yield the prior pdf of the next assimilation cycle.

In general these posterior and prior pdfs are not easily computable. In the Kalman filter, assumptions are notably made on the linearity of operators, to keep these pdfs Gaussian; this way, they are characterized by their mean and covariance matrix. Linear algebra is then sufficient to enforce both Bayes' theorem and the propagation step into operations on means and covariances.

However, with nonlinear models, the Kalman filter assumptions do not hold as the posterior and prior pdfs are not Gaussian anymore. A possibility in this case is to enforce Gaussianity with approximations. This requires the selection of mean and covariances intended for the Gaussian surrogate pdfs. With the Kullback-Leibler divergence, the best Gaus-

<footnote n="1">¹CEREA, Joint Laboratory École des Ponts ParisTech and EDF R&D, Université Paris-Est, Champs-sur-Marne, France</footnote>

<footnote n="2">²CERFACS, Toulouse, France</footnote>

<footnote n="3">³INPT-IRIT, Toulouse, France</footnote>
sian approximation of a pdf is achieved by equating the mean and covariances (see, e.g., Bishop, 2006). However, the integrations necessary to evaluate these moments are also prohibitive.

In the 4D-Var algorithm (see, e.g., Lorenc, 2014, and references therein), the Laplace approximation gives us a way to work around the problem by replacing the posterior mean with the presumed unique minimizer of the cost function over all input values. A model propagation is then sufficient to estimate the prior pdf mean. This global approach calls for efficient global optimization routines. However, in practice, solving a global optimization problem is challenging when the number of unknowns is large, and local methods focused on finding a minimizer over an open subset like Gauss-Newton are often preferred (see, e.g., Björck, 1996).

Unfortunately, the Gauss-Newton method's ability to locate the global minimum depends on the minimization starting point and the cost function properties. Furthermore, missing this global minimum is likely to cause a quick divergence (from the truth) of the sequential DA method. Thus, it is critical for the assimilation algorithm to keep the minimization starting point in a global minimum basin of attraction.

## 1.2 Quasi-static variational data assimilation

Keeping the minimization starting point in a global minimum basin of attraction is constraining because, with a chaotic model, the number of local minima may increase exponentially with the data assimilation window (DAW) time extent $L$ (Pires et al., 1996; Swanson et al., 1998). Unfortunately, assuming a perfect, chaotic – and hence unstable – model, the assimilation performs best for the longest time extents. Several strategies have been investigated to go beyond this restriction. Pires et al. (1996) propose the quasi-static (QS) minimization in a 4D-Var context: as the observations are progressively added to the cost function, the starting point (or first guess) of the 4D-Var minimization is also gradually updated. This led to the method known as QSVA for quasi-static variational assimilation. Ye et al. (2015) propose the gradual increase of the model error covariances in the weak-constraint 4D-Var cost function in a minimization over an entire trajectory; this way the model nonlinearity is gradually introduced into the cost function (see also Judd et al., 2004). They also propose to parallelize this minimization over multiple starting points to increase the chance to locate the global minimum.

On the one hand, 4D-Var benefits from the QS approach to approximate the posterior and prior means. On the other hand, with traditional 4D-Var, the prior covariance matrix is taken as static. This is appropriate when only one cycle of assimilation is considered. But this limits the dynamical transfer of error statistics from one cycle to the next, for instance when Pires et al. (1996) propose gradually moving their fixed-lag data assimilation window in order to build a sequential QSVA.

## 1.3 Ensemble variational methods

By contrast, four-dimensional ensemble variational (EnVar) schemes allow one to perform both a nonlinear variational analysis and a propagation of dynamical errors via the ensemble (see chap. 7 of Asch et al., 2016). The improvement brought by QS minimizations to these schemes has been suggested and numerically evaluated in Bocquet and Sakov (2013, 2014) and Goodliff et al. (2015). This motivates a more complete analytical and numerical investigation.

The iterative ensemble Kalman smoother (IEnKS) (Bocquet and Sakov, 2014; Bocquet, 2016) is the archetypic of such a four-dimensional nonlinear EnVar scheme, where the ensemble parts of the algorithm are deterministic. Using low-order models (usually toy-models), it was shown to significantly outperform 4D-Var, the EnKF or the ensemble Kalman smoother in terms of accuracy.

The IEnKS improves the DA cycling by keeping track of the pdfs' mean and covariance matrix. To do this, a Laplace approximation is used to replace the posterior mean and covariance matrix with the minimizer of the cost function and an approximation of the inverse Hessian at the minimizer, respectively. These moments are then used to update the ensemble statistics. The updated ensemble is next propagated to estimate the prior mean and covariance matrix. Hence, it is also critical for the IEnKS to locate the global minimum of the cost function.

Here, we are interested in the application of the QS minimization to the IEnKS. One of the variants of the IEnKS called the multiple data assimilation (MDA) IEnKS was shown (Bocquet and Sakov, 2014) to be quasi-static by design and can be seen as the EnVar generalization of the sequential QSVA by Pires et al. (1996). It was first tested in Bocquet and Sakov (2013). However, the MDA IEnKS is a specific variant of the IEnKS whereas we see here the IEnKS as an exemplar of deterministic nonlinear 4D EnVar methods.

Goodliff et al. (2015) have applied QSVA numerically to a collection of hybrid and EnVar techniques on the Lorenz 1963 model (Lorenz, 1963), where they vary the magnitude of nonlinearity. Nonetheless the focus of their study was not cycling and the transfer of information from one cycle to the next, which is critical to EnVar methods. They also showed that the ensemble transform Kalman smoother outperforms all of the benchmarked methods. We have furthermore shown that this smoother was systematically outperformed by the IEnKS by design, which was also demonstrated on numerics (Bocquet and Sakov, 2013). This strengthens our claim that the IEnKS can be used here as an exemplar of deterministic nonlinear 4D EnVar methods.

## 1.4 Outline

The rest of the paper is organized as follows. In Sect. 2, the performance dependency of 4D-Var and IEnKS algorithms
on the DAW parameters is investigated. The emphasis is on the transfer of information from one cycle to next, which distinguishes the IEnKS from 4D-Var. In order to do this, a brief presentation of 4D-Var and the IEnKS algorithms is given. We then define a measure of performance for assimilation algorithms. This definition is used to give analytic expressions for the accuracy of both algorithms with a linear, diagonal, autonomous model. This quantifies the impact of cycling on the algorithms. After these preliminaries, the nonlinear, chaotic case is studied. In Sect. 3, we provide and describe the algorithms of the quasi-static IEnKS (IEnKSQS). Section 4 is dedicated to numerical experiments with two Lorenz low-order models and to the improvement of the numerical efficiency. Conclusions are given in Sect. 5.

We emphasize that the algorithmic developments of this study are not meant to improve either high-dimensional or imperfect model data assimilation techniques. Even if Miller et al. (1994) show some similarities between the perfect and imperfect settings of a model of intermediate complexity, model error would generally forbid the use of very long DAWs as sometimes considered in this study. Instead, the objective of this paper is to better understand the interplay between chaotic dynamics, ensemble variational data assimilation schemes and their cycling, irrespective of whether they could be useful in high-dimensional systems.

## 2 The data assimilation window and assimilation performance

After reviewing 4D-Var (in a constant background matrix version) and the IEnKS algorithms, we investigate the dependency of assimilation performance on the DAW key parameters. This illustrates the cycling improvement brought in by the IEnKS compared to 4D-Var. We will see that, with chaotic models, the longer the DAW is, the better the accuracy of these algorithms. Which highlights the QSVA relevance for cycled data assimilation.

The evolution and observation equations of the system are assumed to be in the following form:

$$y_l = \mathcal{H}(x_l) + \varepsilon_l, \quad (1a)$$

$$x_{l+1} = \mathcal{M}(x_l), \quad (1b)$$

where the unknown state $x_l$ at time $t_l$ is propagated to $t_{l+1}$ with the model resolvent $\mathcal{M}: \mathbb{R}^m \to \mathbb{R}^m$. The model is assumed to be perfect (no errors in Eq. (1b)) and autonomous ($\mathcal{M}$ does not depend on time). The observation operator $\mathcal{H}: \mathbb{R}^m \to \mathbb{R}^d$ relates the state $x_l$ to the observation vector $y_l$. The observation errors $(\varepsilon_l)_{l \ge 0}$ are assumed to be Gaussian with mean $\mathbf{0} \in \mathbb{R}^d$ and covariance matrix $\mathbf{R} \in \mathbb{R}^{d \times d}$; they are uncorrelated in time.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic of a DAW. The state variable at $t_0$ is $x_0$, the observation vector $y_K$ at time $t_K$ is the first of the DAW to be assimilated and the observation $y_L$ at present time $t_L$ is the last to be assimilated. These observations, possibly observation vectors, are represented by black dots.</figcaption></figure>

### 2.1 4D-Var and IEnKS algorithms

Both 4D-Var and the IEnKS use a variational minimization in their analysis step. The objective of this minimization is to locate the global maximum of the posterior pdf $p(x_0|y_{L:K})$ of the system past state $x_0$ given the observations $y_{L:K} = [y_K, \dots, y_L]$ at times $t_{L:K} = [t_K, \dots, t_L] \in \mathbb{R}^S$. The system state and observations are seen as random vectors with values in $\mathbb{R}^m$ and $\mathbb{R}^d$, respectively. The posterior pdf quantifies how our knowledge on the state $x_0$ changes with realizations of $y_{L:K}$. Thus, its maximum is the most probable state after assimilating the observations. The DAW is displayed in Fig. 1. The parameters $K$ and $L$ are the time index of the DAW first and last assimilated observation batch, respectively. The number of observation vectors used within the DAW is $L-K+1$.

To specify this posterior pdf, we have to make further assumptions on $x_0$.

The initial state $x_0$ is assumed to be Gaussian with mean $\mathbf{x}_0^b \in \mathbb{R}^m$ and covariance matrix $\mathbf{B} \in \mathbb{R}^{m \times m}$:

$$p(x_0) = \mathcal{N}(\mathbf{x}_0|\mathbf{x}_0^b, \mathbf{B}). \quad (2)$$

Following these assumptions, an analytic expression can be obtained for the posterior smoothing pdf $p(\mathbf{x}_0|\mathbf{y}_{L:K})$ at the first cycle, or for the cost function associated with this pdf. The latter is defined as

$$G(\mathbf{x}_0|\mathbf{y}_{L:K}) = -\ln p(\mathbf{x}_0|\mathbf{y}_{L:K}). \quad (3)$$

The notation $G$ is used rather than the traditional $J$ to refer to an exact cost function, i.e., a cost function defined from an exact posterior pdf. Bayes' theorem consequently yields the following:

$$G(\mathbf{x}_0|\mathbf{y}_{L:K}) = -\ln p(\mathbf{x}_0) - \ln p(\mathbf{y}_{L:K}|\mathbf{x}_0) \\ + \ln p(\mathbf{y}_{L:K}), \quad (4)$$

and with the Gaussian assumption applied to the background and observation errors we find
$$
G(\mathbf{x}_0|\mathbf{y}_{L:K}) = \frac{1}{2} \|\mathbf{x}_0^b - \mathbf{x}_0\|_{\mathbf{B}^{-1}}^2 + c_0 \\
+ \frac{1}{2} \sum_{l=K}^{L} \|\mathbf{y}_l - \mathcal{H} \circ \mathcal{M}^l(\mathbf{x}_0)\|_{\mathbf{R}^{-1}}^2, \qquad (5)
$$

where $\|\mathbf{x}\|_{\mathrm{A}}^{2} = \mathbf{x}^{\mathrm{T}}\mathbf{A}\mathbf{x}$ is the norm of $\mathbf{x}$ associated with a sym-
metric positive definite matrix $\mathbf{A}$; $c_0$ is a normalization con-
stant ensuring $\int e^{-G(\mathbf{x}_0|\mathbf{y}_{L:K})} d\mathbf{x}_0 = 1$; $\mathcal{M}^l$ stands for $l$ com-
positions of $\mathcal{M}$.

The propagation corresponds to a time shift of S time
steps. Thus, at the kth assimilation cycle, the posterior pdf
is p(x_{KS}|y_{KS+L:K}). Using Bayes' theorem the kth cycle cost
function is

$$
G(\mathbf{x}_{KS}|\mathbf{y}_{KS+L:K}) = -\ln p(\mathbf{x}_{KS}|\mathbf{y}_{(k-1)S+L:K}) \\
+ \frac{1}{2} \sum_{l=K}^{L} \left\| \mathbf{y}_{KS+l} - \mathcal{H} \circ \mathcal{M}^{l}(\mathbf{x}_{KS}) \right\|_{\mathbf{R}^{-1}}^{2} + c_{KS}, \quad (6)
$$

where $-\ln p(\mathbf{x}_{KS}|\mathbf{y}_{(k-1)S+L:K})$ is the background term. If the model and observation operators are nonlinear, an analytical expression for this prior is not accessible and one needs to approximate it. The 4D-Var and the IEnKS algorithms are solutions based on distinct approximation strategies.

The 4D-Var cost function at the kth cycle, based on the
static error covariance matrix **B**, is defined by

$$
J(\mathbf{x}_{KS}; \mathbf{y}_{KS+L:kS+K}, \mathbf{x}_{KS}^b) = \frac{1}{2} \|\mathbf{x}_{KS}^b - \mathbf{x}_{KS}\|_{\mathbf{B}^{-1}}^2 \\
+ \frac{1}{2} \sum_{l=K}^{L} \left\| \mathbf{y}_{KS+l} - \mathcal{H} \circ \mathcal{M}^l(\mathbf{x}_{KS}) \right\|_{\mathbf{R}^{-1}}^2. \quad (7)
$$

The analysis of 4D-Var consists in by minimizing Eq. (7), yielding $\mathbf{x}_{KS}^a$ at $t_{KS}$. Because this cost function depends on realizations of the random observations, $\mathbf{x}_{KS}^a$ is also a random variable. This analysis is then propagated at time $t_{(k+1)S}$ with the resolvent of the model to produce the next cycle background state:

$$
\mathbf{x}_{(k+1)S}^{b} = \mathcal{M}^{S}(\mathbf{x}_{KS}^{a}). \qquad (8)
$$

In general, $G(\mathbf{x}_0|\mathbf{y}_{L:K})$ and $J(\mathbf{x}_0; \mathbf{y}_{L:K}, \mathbf{x}_0^b)$ only coincide at the first cycle of an assimilation because of the assumption outlined by Eq. (2). Subsequently, the background term of the 4D-Var cost function $\frac{1}{2}\|\mathbf{x}_{KS}^b - \mathbf{x}_{KS}\|_{\mathbf{B}^{-1}}^2$ is a Gaussian approximation of the exact background term $-\ln p(\mathbf{x}_{KS}|\mathbf{y}_{(k-1)S+L:K})$. By definition, the background error covariance matrix $\mathbf{B}$ of the traditional 4D-Var cost function is the same for each cycle – this is not the case for the IEnKS.

The IEnKS (Bocquet and Sakov, 2014) is an ensemble
method with a variational analysis. Two versions of the al-
gorithm exist: the singular data assimilation (SDA) version

where observations are assimilated only once, and the multi-
ple data assimilation (MDA) version where they are assimi-
lated several times. We focus on the SDA version in the the-
oretical development. The MDA version, which can be seen
as the first published quasi-static EnVar scheme is used in
the numerical experiments for comparison. Note that, for the
SDA IEnKS, the number of observations is L − K + 1 = S .

At the kth cycle of the IEnKS, the background ensemble at $t_{KS}$ is obtained by a propagation from the previous cycle. The ensemble members are the columns of the matrix $\mathbf{E}_{KS}^b$, which is seen as a random matrix with values in $\mathbb{R}^{m \times n}$. It is used to estimate the prior mean and covariance matrix as follows:

$$
\mathbb{E}[\mathbf{x}_{KS} | \mathbf{y}_{(k-1)S+L:K}] \approx \bar{\mathbf{x}}_{KS}^{b}, \qquad (9)
$$

$$
\mathbb{C}[\mathbf{x}_{KS} | \mathbf{y}_{(k-1)S+L:K}] \approx \mathbf{X}_{KS}^{b} \mathbf{X}_{KS}^{b\mathrm{T}}, \quad (10)
$$

where $\mathbb{E}$ and $\mathbb{C}$ are the expectation and covariance operators, respectively; $\bar{\mathbf{x}}_{KS}^b$ and $\mathbf{X}_{KS}^b$ are the empirical mean and normalized anomaly of $\mathbb{E}_{KS}^b$, respectively:

$$
\bar{\mathbf{x}}_{KS}^{b} = \mathbf{E}_{KS}^{b} \frac{\mathbf{1}_{n}}{n}, \qquad (11)
$$

$$
\mathbf{X}_{KS}^{b} = \mathbf{E}_{KS}^{b} \frac{\mathbf{I}_{n} - \frac{1}{n} \mathbf{I}_{n}^{\mathrm{T}}}{\sqrt{n-1}}, \qquad (12)
$$

with $\mathbf{1}_n = [1, \dots, 1]^T \in \mathbb{R}^n$ a vector of ones and $\mathbf{I}_n$ is the iden-
tity of $\mathbb{R}^n$. Note that Eqs. (9) and (10) are approximations
because of sampling errors. If the state vector is of the form

$$
\mathbf{x}_{KS} = \bar{\mathbf{x}}_{KS}^{b} + \mathbf{X}_{KS}^{b} \mathbf{w}_{KS}, \qquad (13)
$$

with $\mathbf{w}_{KS} \in \mathbb{R}^n$ the control variable in the ensemble space,
the IEnKS cost function is defined in the ensemble space by

$$
J(\mathbf{w}_{KS}; \mathbf{y}_{KS+L:kS+K}, \mathbf{E}_{KS}^{b}) = \frac{1}{2} \|\mathbf{w}_{KS}\|^2 \\
+ \frac{1}{2} \sum_{l=K}^{L} \left\| \mathbf{y}_{KS+l} - \mathcal{H} \circ \mathcal{M}^{l} (\bar{\mathbf{x}}_{KS}^{b} + \mathbf{X}_{KS}^{b} \mathbf{w}_{KS}) \right\|_{\mathbf{R}^{-1}}^2, \quad (14)
$$

where $\|\cdot\| = \|\cdot\|_{\mathbf{I}}$ with $\mathbf{I}$ the identity matrix. The analysis of the IEnKS consists in minimizing this cost function. It yields an updated ensemble $\mathbf{E}_{KS}^a$ at time $t_{KS}$ verifying

$$
\bar{\mathbf{x}}_{KS}^{a} = \bar{\mathbf{x}}_{KS}^{b} + \mathbf{X}_{KS}^{b} \mathbf{w}_{KS}^{a}, \qquad (15)
$$

$$
\mathbf{X}_{KS}^{a} = \mathbf{X}_{KS}^{b} \left[ \nabla^2 J \left( \mathbf{w}_{KS}^{a}; \mathbf{y}_{KS+L:kS+K}, \mathbf{E}_{KS}^{b} \right) \right]^{-1/2} \mathbf{U}, \quad (16)
$$

with $\bar{\mathbf{x}}_{KS}^a$ and $\mathbf{X}_{KS}^a$ representing the empirical mean and normalized anomalies of the analyzed ensemble, respectively, $\mathbf{w}_{KS}^a$ the cost function minimizer, $\nabla^2 J(\mathbf{w}_{KS}^a; \mathbf{y}_{KS+L:kS+K}, \mathbf{E}_{KS}^b)$ its Hessian at this minimum usually approximated in the Gauss–Newton method by

$$
\nabla^2 J(\mathbf{w}_{KS}^a; \mathbf{y}_{KS+L:kS+K}, \mathbf{E}_{KS}^b) \approx \mathbf{I}_n + \sum_{l=K}^{L} \mathbf{F}_l^T \mathbf{R}^{-1} \mathbf{F}_l, \quad (17)
$$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Chaining of the 4D-Var and IEnKS first two cycles with $S=3, L=4$. The first 4D-Var analysis uses the background $\mathbf{x}_0^b$ at $t_0$ and the observations $y_{L,K}$ to give the analysis $\mathbf{x}_0^a$ at $t_0$. It is propagated $S$ steps forward in time to produce the new background at time $t_S$ where another analysis can be performed. The IEnKS does the same but with an ensemble. The dashed rectangle symbolizes the current DAW, black dots represent the observations assimilated in the current cycle, gray dots represent already assimilated observations and white dots represent observations not assimilated.</figcaption></figure>

with $\mathbf{F}_l = \frac{d\mathbf{F}_l}{d\mathbf{w}_{kS}}(\mathbf{w}_{kS}^a)$, where $\mathcal{F}_l: \mathbf{w}_{kS} \mapsto \mathcal{H} \circ \mathcal{M}^l(\bar{\mathbf{x}}_{kS}^b + \mathbf{X}_{kS}^b \mathbf{w}_{kS})$, in order to avoid computing the model second derivatives. The exponent $-1/2$ refers to the unique symmetric definite positive inverse square root of a symmetric definite positive matrix and $\mathbf{U} \in \mathbb{R}^{n \times n}$ an orthogonal matrix such that $\mathbf{U} \mathbf{1}_n = \mathbf{1}_n$. This analyzed ensemble is also propagated to time $t_{(k+1)S}$ to produce the next cycle background ensemble

$$ \mathbf{E}_{(k+1)S}^{b} = \mathcal{M}^{S}(\mathbf{F}_{kS}^{a}). \qquad (18) $$

The cycle is completed by using the cost function Eq. (14) with time indexes incremented by $S$ in the next analysis. 4D-Var and the IEnKS are displayed in Fig. 2.

## 2.2 Performance of assimilation

In order to evaluate the efficiency of 4D-Var and the IEnKS with DAW parameters $S$ and $L$, two measures of accuracy are investigated here: the usual empirical RMSE, and a theoretical counterpart.

At the $k$th cycle, the algorithm generates at time $t_{kS}$ an analysis $\mathbf{x}_{kS}^a$ from the observations. This analysis is propagated with the model $l$ steps forward in time to yield the analysis $\mathbf{x}_{kS+l}^a$ meant to approximate the system's true state $\mathbf{x}_{kS+l}$. A traditional measure of assimilation performance is the root mean square error (RMSE). It is defined by

$$ \text{RMSE} = \frac{1}{\sqrt{m}} \| \mathbf{x}_{kS+l} - \mathbf{x}_{kS+l}^a \|. \qquad (19) $$

The RMSE takes different names depending on when it is computed. If $l=L$, it is called the filtering RMSE; if $0 \le l \le L-1$, it is the smoothing RMSE with lag $L-l$. In the following, the smoothing RMSE will correspond to the one with (maximum) lag $L$.

The RMSE rigorously depends on the random variable realizations, and thus it is also a random variable. In our numerical experiments, as is usually the case, the RMSE is averaged over the cycles to mitigate this variability:

$$ a\text{RMSE}_N = \frac{1}{N} \sum_{k=0}^{N-1} \frac{1}{\sqrt{m}} \| \mathbf{x}_{kS+l} - \mathbf{x}_{kS+l}^a \|. \qquad (20) $$

Let us assume that there is a random couple $(\mathbf{x}_{\infty S+l}, \mathbf{x}_{\infty S+l}^a)$ whose distribution is invariant and ergodic with respect to the shift transformation:

$$ T: (\mathbf{x}_{kS+l}, \mathbf{x}_{kS+l}^a) \mapsto (\mathbf{x}_{(k+1)S+l}, \mathbf{x}_{(k+1)S+l}^a). \qquad (21) $$

Then by Birkhoff's ergodic theorem (see Walters, 1982) the sequence $(a\text{RMSE}_N)_N$ converges when $N \to \infty$ and its limit aRMSE verifies

$$ a\text{RMSE} = \frac{1}{\sqrt{m}} \mathbb{E} \left[ \| \mathbf{x}_{\infty S+l} - \mathbf{x}_{\infty S+l}^a \| \right], \qquad (22) $$

where the expectation $\mathbb{E}$ is taken over $p(\mathbf{x}_{\infty S+l}, \mathbf{x}_{\infty S+l}^a)$. In this case, the aRMSE measures the long-term impact of the cycling on the assimilation accuracy. This limit is difficult to exploit algebraically. That is why, in the theoretical developments, we prefer the expected MSE (eMSE), denoted by $P$ as follows:

$$ P_{kS+l} = \mathbb{E} \left[ \| \mathbf{x}_{kS+l} - \mathbf{x}_{kS+l}^a \|^2 \right], $$

where the expectation is taken over $p(\mathbf{x}_{kS+l}, \mathbf{x}_{kS+l}^a)$. In the following subsection, we will focus on the long term impact of the cycling on $P$. Simplifying assumptions will be made to express $P_{\infty S+l} = \lim_{k \to \infty} P_{kS+l}$ as a function of $S$ and $L$.

## 2.3 Performance in the linear, diagonal, autonomous case

In order to obtain analytical expressions of the eMSE for 4D-Var and the IEnKS, we make drastic simplifying assumptions.

First, the model is assumed to be the resolvent of a linear, diagonal, autonomous ordinary differential equation. Thus, it can be expressed as

$$ \mathcal{M}^l(\mathbf{x} + \delta\mathbf{x}) = \mathcal{M}^l(\mathbf{x}) + \mathbf{M}^l \delta\mathbf{x}, \qquad (23) $$

where $\mathbf{M} = \text{diag}(\alpha_i)_{i=1..m}$ is diagonal and does not depend on $\mathbf{x}$. We further assume $\mathcal{H} = h\mathbf{I}_m, \mathbf{B} = b\mathbf{I}_m, \mathbf{R} = r\mathbf{I}_m$, where $\mathbf{I}_m$ is the identity matrix of $\mathbb{R}^m$ and $h, r, b > 0$. With these assumptions, Appendix A provides an expression for the 4D-Var asymptotic eMSE in the univariate case. The generalization to the diagonal multivariate case is obtained by summing
up the eMSEs of each direction:

$$P_{\infty S+l}^{\text{4D-Var}} = \sum_{i=1}^{m} \begin{cases} \infty & \text{if } \Delta_i \ge 1 \\ \frac{b^2 \Sigma_{K,i}^L}{\alpha_i^{2(S-l)}} \frac{\Delta_i}{1 - \Delta_i} & \text{otherwise,} \end{cases} \quad (24a)$$

$$\Sigma_{K,i}^{L} = \frac{h^2 \alpha_i^{2(L+1)} - \alpha_i^{2K}}{r \alpha_i^2 - 1}, \quad (24b)$$

$$\Delta_i = \frac{\alpha_i^{2S}}{(1 + b \Sigma_{K,i}^L)^2}. \quad (24c)$$

The case $\Delta_i \ge 1$ means that too much credit is given to the background variance, which is approximated for all cycles by the constant $b$ in our 4D-Var scheme. Therefore, the information carried by the observations is not sufficient to mitigate the exponential growth of errors in the propagation.

Concerning the IEnKS, the anomalies are assumed to be full rank to avoid any complication due to singular covariance matrices. Moreover, the linearity of the model is employed to express the background statistics:

$$\bar{\mathbf{x}}_{(k+1)S}^{\mathrm{b}} = \mathcal{M}^{S}(\bar{\mathbf{x}}_{kS}^{\mathrm{a}}), \quad (25a)$$

$$\mathbf{X}_{(k+1)S}^{\mathrm{b}} = \mathbf{M}^{S} \mathbf{X}_{kS}^{\mathrm{a}}, \quad (25b)$$

$$\bar{\mathbf{x}}_{0}^{\mathrm{b}} = \mathbf{x}_{0}^{\mathrm{b}}, \quad (25c)$$

$$\mathbf{X}_{0}^{\mathrm{b}} = \mathbf{B}^{1/2}. \quad (25d)$$

This way, sampling errors are avoided as they are not the focus of this study. The background ensemble $\mathbf{E}_{kS}^{\mathrm{b}}$ becomes a notational shortcut for the pair $(\bar{\mathbf{x}}_{kS}^{\mathrm{b}}, \mathbf{X}_{kS}^{\mathrm{b}})$. This simplified IEnKS is actually a Kalman smoother (Cosme et al., 2012; Bocquet and Carrassi, 2017). With these assumptions, Appendix B gives an expression for the IEnKS asymptotic eMSE in the univariate case. The optimality of the IEnKS eMSE is also proven. The generalization to the diagonal multivariate case is also obtained by summing up the eMSEs of each direction:

$$P_{\infty S+l}^{\text{IEnKS}} = \sum_{i=1}^{m} \begin{cases} 0 & \text{if } |\alpha_i| \le 1 \\ \frac{r}{h^2 \alpha_i^{2(L-l)}} \frac{\alpha_i^2 - 1}{\alpha_i^2} & \text{otherwise.} \end{cases} \quad (26)$$

This expression shows that the eMSE components on the stable directions are null. Indeed one expects the IEnKS to be at least more efficient than a free-run, whose errors in the stable directions tend to zero.<fnref n="1" /> This is not the case for 4D-Var since the static background covariance matrix introduces spurious variance in the stable directions as seen in Eq. (24). In Trevisan et al. (2010), 4D-Var error variances in the stable directions are forced to zero to improve the accuracy of the assimilation.

In the following, we study the eMSE dependency on the DAW parameters, $S$ and $L$, for both algorithms. We focus on a bivariate case with $\alpha_1 = 1.2$, $\alpha_2 = 0.8$ in order to have one stable and one unstable direction; $h = b = r = 1$. Using Eqs. (24) and (26), the asymptotic smoothing and filtering eMSEs are displayed as a function of $L, S$ (Fig. 3). The eMSE components on the stable and unstable directions are also shown. Those graphs are interpreted in the following.

Specifically, the eMSE expression for the IEnKS is of the form

$$P_{\infty S+l}^{\text{IEnKS}} = c_1 \alpha_1^{2(l-L)}, \quad (27)$$

where $c_1$ does not depend on $S, L, l$. The contribution to $P_{\infty S+l}^{\text{IEnKS}}$ on the stable direction is zero. It depends only on the lag $L-l$; it does not depend on $S$. Moreover, the filtering eMSE ($l=L$) is constant: the propagation compensates for the analysis. The smoothing eMSE ($l=0$) decreases exponentially with $L$.

Concerning 4D-Var, we assume $K$ to be fixed and $S \to \infty$ to give an asymptotic eMSE expression:

$$P_{\infty S+l}^{\text{4D-Var}} = P_{\infty S+l}^{\text{IEnKS}} (1+o(1)) + c_2 \alpha_2^{2l} (1+o(1)), \quad (28)$$

where $c_2$ is constant with $S, l$ and $o(1) \to 0$ when $S \to \infty$. The unstable component is close to the IEnKS overall eMSE. The biggest difference with it concerns the eMSE on the stable component. The inexact background variance modeling adds a detrimental term to the eMSE.

To qualify the long term impact of the cycling on the errors, the filtering eMSE is more instructive than the smoothing eMSE. Indeed, the smoothing eMSE is improved with $L$ as it adds future observations (with respect to the analysis time) in the DAW. This improvement dominates the potentially detrimental impact of the Gaussian background approximations. In the filtering eMSE, this improvement is balanced by the propagation of the analysis at the end of the DAW. In Fig. 3, the filtering eMSE stable component is mitigated such that it has little effect. However, the parameter $S$ improves the filtering eMSE on the unstable component. The bigger $S$ becomes, the closer $P_{\infty S+l}^{\text{4D-Var}}$ is to $P_{\infty S+l}^{\text{IEnKS}}$, which is optimal (cf. Appendix B). A qualitative explanation is that to assimilate the same numbers of observations, a 4D-Var using high values of $S$ needs fewer cycles; therefore, it less often relies on the background approximation, making the analysis more trustworthy.

Figure 4 displays the asymptotic eMSEs of both algorithms as a function of the lag for $S = L = 5$;

these curves are similar to those of Trevisan et al. (2010). Concerning 4D-Var, the eMSE can be written as a function of the lag $L-l$ as follows:

$$P_{\infty S+l}^{\text{4D-Var}} \approx c_1 \alpha_1^{-2\text{lag}} + c_2 \alpha_2^{2(L-\text{lag})}. \quad (29)$$

Thus, the unstable component of the eMSE is an exponentially decreasing function of the lag and the stable component

<footnote n="1">¹The fact that the errors lie in the unstable subspace is more general. (Bocquet and Carrassi, 2017)</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The asymptotic smoothing (a) and filtering (b) eMSEs of 4D-Var and the IEnKS. eMSEs are displayed as functions of *L*, *S* with their components on the stable and unstable subspaces. The *L* parameter is on the abscissa axis for 4D-Var and on the ordinate axis for the IEnKS. The *S* parameter is on the abscissa axis for the IEnKS and on the ordinate axis for 4D-Var.</figcaption></figure>

is an exponentially increasing function of the lag. The sum therefore decreases when the unstable component is dominant and increases when the stable component is dominant.

In this section, we studied the accuracy of both the cycled IEnKS and 4D-Var eMSEs as a function of the DAW parameters under linear, autonomous, diagonal assumptions. We found that the DAW parameter *L* improves the smoothing eMSE and *S* improves the filtering eMSE. These properties will be discussed and numerically investigated in a nonlinear context in the next section.

## 2.4 Performance in the nonlinear, chaotic case

The results of Sect. 2.3 foster the use of the largest possible *L* to improve the 4D-Var and IEnKS smoothing eMSEs. For filtering, the error propagation at the end of the DAW

balances the gain in eMSE due to the assimilation of future observations. Thus, the filtering eMSE is not improved by the assimilation of observations distant in time; it is affected by the background pdf approximation. To assimilate the same number of observations, algorithms using high values of *S* need fewer cycles. Therefore, they less often rely on the background approximation. That is why the 4D-Var filtering performance is improved with *S*. For the IEnKS as in Sect. 2.3, the prior pdf approximation is exact so that there is no dependence on *S*. This is no longer true in the nonlinear case. Hence, with a similar reasoning, one actually expects an improvement on the IEnKS filtering performance with *S*. As a matter of fact, it has been shown that with a nonlinear chaotic model, the filtering accuracy increases with *L* in most cases (see Bocquet and Sakov, 2014, and Sect. 4 of the present paper).
<figure><img src="image_7.png" /><figcaption>Figure 4. The asymptotic eMSEs as a function of the lag (lag = 0 is the filtering performance and lag = L is the smoothing performance). The superposing curves have been slightly translated for better readability, S = L = 5, h = b = r = 1.</figcaption></figure>

## 2.4.1 Multiple local minima

However, with a chaotic model, Pires et al. (1996) show that the 4D-Var cost function number of local extrema increases with L, making minimization problematic. We show in this section that the IEnKS cost function suffers from the same problem.

This behavior will be illustrated with the Lorenz 95 (L95) model (Lorenz and Emanuel, 1998). It represents a mid-latitude zonal circle of the atmosphere and is described by a set of m nonlinear differential equations:

$$ \frac{dx_j}{dt} = (x_{j+1} - x_{j-2})x_{j-1} - x_j + F, \quad (30) $$

where $x_j$ is the j-th modulo m component of $x$, $m = 40$ and $F = 8$. This equation is integrated using a fourth-order Runge-Kutta scheme with a time step of $\delta t = 0.05$. The dynamics of L95 are chaotic; the L95 model's largest Lyapunov exponent is $\lambda \simeq 1.7$.

Figure 5 shows a typical IEnKS cost function profile in one direction of the analyzed ensemble space for multiple values of the DAW parameters. The system is observed at every time step and $\mathcal{H} = \mathbf{B} = \mathbf{R} = \mathbf{I}_m$.

The curves have more and more local extrema as L increases. The curves with the highest amplitudes of the ripples are found for small values of S. Indeed, an averaging effect may settle in as the number of observations increases.

This hilly shape causes minimization problems. A possible minimization procedure for the IEnKS is the Gauss-Newton (GN) algorithm (e.g., Björck, 1996). GN is not a global procedure meaning that, depending on the starting point and the cost function properties, the algorithm can con-

verge towards any local extremum, take many iterations or even diverge. However, if the cost function is quadratic, the global minimum is reached in one iteration. The cost function non-quadraticity is induced by the model nonlinearity. In the following, we will give a heuristic argument yielding a bound on the S parameter beyond which the GN method probably misses the global minimum in the configuration $\mathcal{H} = \mathbf{B} = \mathbf{R} = \mathbf{I}_m$. To some extent, our argument can be seen as an improvement on the notion of useful DAW length introduced by Pires et al. (1996) beyond which the performance gain is negligible. In contrast with the useful DAW length, we account in the following for both the cycling and the non-linearity.

## 2.4.2 Effective data assimilation window length

First, the GN convergence properties are drastically simplified. We assume the method converges to the global minimum if, and only if, the minimization starting point is in a neighborhood of the global minimizer where the IEnKS cost function is almost quadratic.

Unfortunately, this minimizer is unknown because the cost function depends on realizations of many random variables. In order to eliminate this variability, Pires et al. (1996) introduced a so-called error-free cost function. We prefer an averaged cost function $J_{\infty S}$ in this particular study defined by

$$ J_{\infty S}(\mathbf{w}) = \lim_{N \to \infty} \frac{1}{N} \sum_{k=0}^{N-1} J(\mathbf{w}; y_{kS+L:kS+K}, \mathbf{E}_{kS}^{\mathbf{b}}). \quad (31) $$

Relying on an ergodicity assumption, Appendix C proves that this averaged cost function verifies

$$ J_{\infty S}(\mathbf{w}) = \frac{1}{2} \|\mathbf{w}\|^2 + \frac{dS}{2} + \frac{1}{2} \sum_{l=K}^{L} \mathbb{E} \left[ \left\| \delta \mathbf{x}_{\infty S+l}^{\mathbf{b}} \right\|^2 \right], \quad (32a) $$

$$ \delta \mathbf{x}_{\infty S+l}^{\mathbf{b}} = \mathcal{M}^l \left( \bar{\mathbf{x}}_{\infty S}^{\mathbf{b}} + \mathbf{X}_{\infty S}^{\mathbf{b}} \mathbf{w} \right) - \mathcal{M}^l (\mathbf{x}_{\infty S}), \quad (32b) $$

where the ergodic random variables $\bar{\mathbf{x}}_{\infty S}$, $\mathbf{X}_{\infty S}^{\mathbf{b}}$ and $\mathbf{x}_{\infty S}$ have been defined in Sect. 2.2 and Appendix C. As seen in Eq. (32), a sufficient condition for the starting point $\mathbf{w} = \mathbf{0}$ to be in a neighborhood of the global minimizer where the cost function is assumed almost quadratic requires that $\bar{\mathbf{x}}_{\infty S}^{\mathbf{b}}$ be in a neighborhood of $\mathbf{x}_{\infty S}$ where all the $(\mathcal{M}^l)_{K \le l \le L}$ are almost linear.

In the univariate case, if the model behavior is almost linear and unstable, we can use Eq. (26) to estimate the terms in the sum in Eq. (32) at the starting point $\mathbf{w} = \mathbf{0}$:

$$ \mathbb{E} \left[ \left\| \delta \mathbf{x}_{\infty S+l}^{\mathbf{b}} \right\|^2 \right] \approx \alpha^{2(S+l-L)} \frac{\alpha^2 - 1}{\alpha^2}, \quad (33) $$

where $\alpha$ is the linear part of the model and the extra S accounts for the propagation. But in a necessarily bounded
<figure><img src="image_8.png" /><figcaption>Figure 5. Cost functions of the IEnKS projected in one direction of the analyzed ensemble (hence centered and normalized) with various DAW parameters S, L. A quasi-static minimization can be visualized from these panels. It begins with the bottom-left cost function; the orange square is the starting point and the green dot is at the minimum. From the bottom-left cost function up to the top-right cost function, batches of 9 then 10 observation vectors are progressively added to the DAW, and the minimizer (green dot) is updated accordingly.</figcaption></figure>

physical system, the right-hand side of Eq. (33) cannot grow indefinitely with $l + S$. Such model saturation imposes

$$ \mathbb{E} \left[ \left\| \delta \mathbf{x}_{\infty S+l}^{\mathrm{b}} \right\|^2 \right] \leq B, \qquad (34) $$

with B representing a bound. Hence, Eqs. (33) and (34) yield the following inequality on S:

$$ S \leq S_{\max}, \qquad (35a) $$

$$ S_{\max} = \frac{\ln(B) - \ln(1 - \alpha^{-2})}{2\ln(\alpha)}. \qquad (35b) $$

We choose $l = L$ because it corresponds to the most constraining case. When Eq. (35a) is violated, $\bar{\mathbf{x}}_{\infty S}^{\mathrm{b}}$ departs from $\mathbf{x}_{\infty S}$ such that the nonlinearities of $\mathcal{M}^L$ are significant.

To apply this inequality to the L95 model we choose the following:

$$ \alpha = \lim_{N \to \infty} \frac{1}{N} \sum_{k=0}^{N-1} \sigma \left( \frac{d\mathcal{M}}{dx}(\mathbf{x}_k) \right), \qquad (36) $$

$\sigma$ being the mean of the singular values greater than 1. This corresponds to an average of the error amplification by $\frac{d\mathcal{M}}{dx}(\mathbf{x}_{\infty S})$ in the local unstable subspace. We also choose

the average squared norm between two long trajectories for the bound B:

$$ B = \lim_{N \to \infty} \frac{1}{mN} \sum_{k=0}^{N-1} \left\| \mathbf{x}_{kS} - \mathcal{M}^{kS} \left( \bar{\mathbf{x}}_0^{\mathrm{b}} \right) \right\|^2. \qquad (37) $$

This quantity is greater than $\mathbb{E} \left[ \left\| \delta \mathbf{x}_{\infty S+l}^{\mathrm{b}} \right\|^2 \right]$ because the IEnKS asymptotic performance is at least better than a free run. From the values of $\alpha$ and B we find $S_{\max} = 14$.

Figure 6 shows the filtering and smoothing aRMSEs of an IEnKS $L = S$ with L95 as a function of S; the performance strongly deteriorates for $S > 16$, which is remarkably consistent with our estimation.

Figure 6 also shows another difference with the linear case: the IEnKS filtering aRMSE depends on S. The former discussion on the local extrema explains this dependency for large values of S.

However, for small values of S, the decreasing aRMSE has not been explained. This is a consequence of the Gaussian background approximation. At each cycle, the IEnKS uses the background ensemble $\mathcal{E}_{kS}^{\mathrm{b}}$ to estimate the first two moments of the background pdf and makes the following ap-
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Smoothing aRMSE (a) and filtering aRMSE (b) of a Gauss–Newton IEnKS $L = S$ as a function of $S$ with the L95 model over $5 \times 10^5$ cycles (logarithmic scale). We used the finite-size version (Bocquet, 2011) to account for sampling errors and avoid the need for inflation, $\mathcal{H} = \mathbf{R} = \mathbf{B} = \mathbf{I}_m, n = 20$.</figcaption></figure>

proximation:

$$p(\mathbf{x}_{kS} | \mathbf{y}_{(k-1)S+L:K}) \approx \mathcal{N} \left( \mathbf{x}_{kS} | \bar{\mathbf{x}}_{kS}^{\mathrm{b}}, \mathbf{X}_{kS}^{\mathrm{b}} (\mathbf{X}_{kS}^{\mathrm{b}})^{\mathrm{T}} \right). \quad (38)$$

Because the model is nonlinear, this pdf is unlikely to be Gaussian. Therefore, this approximation results in a loss of information and the more it is used, the farther from $G$ the IEnKS cost function is. This is exactly what happens when $S$ is small: to assimilate the same number of observations, the IEnKS uses more cycles so that it relies more on the Gaussian background approximation.

### 3 Quasi-static algorithms

We have seen in the previous section that the effective DAW length is constrained by the cost function non-quadraticity. In this section we review and propose algorithms able to overcome these minimization issues and reach longer DAWs.

Quasi-static algorithms were introduced by Pires et al. (1996) in a 4D-Var context. The idea behind QSVA is to control the way observations are added to the cost function in order to keep the minimization starting point in the basin of attraction containing the global minimizer. The method is carried out by repeating the minimization with an increasing number of observations: the first minimization is performed using the cost function with a single observation vector, then the number of observation vectors is increased and another minimization can be performed with the former minimizer as the new starting point. The process is then repeated until all observations are accounted for.

This procedure is directly applicable to the IEnKS cost function minimization. The left panel in Fig. 7 is a schematic of a QS minimization and Algorithm 1 gives the pseudo-code of an IEnKS with a QS minimization. The new parameters $(L_q)_{q<N_Q}$ control the number of observations added at each minimization, where $N_Q$ is the total number of batches of observation vectors.

The first three lines initialize the minimization starting point, the ensemble mean and anomaly matrix. The **for** loop in lines 4–23 repeats the QS minimization. The **while** loop in lines 6–22 is the Gauss–Newton minimization. Lines 7 and 8 center the ensemble on the current minimizer. Lines 9 and 10 initialize the cost function gradient and the approximate Hessian. The **for** loop in lines 12–18 compute the observation terms of the cost function, the gradient and the approximate Hessian. Lines 13 and 14 use a finite difference formula to compute the tangent linear and adjoint of $\mathbf{w} \mapsto \mathcal{H} \circ \mathcal{M}^l (\bar{\mathbf{x}}_0^\mathrm{b} + \mathbf{X}_0^\mathrm{b} \mathbf{w})$. Lines 15 and 16 use this adjoint to update the gradient and approximate the Hessian. Lines 19 and 20 solve the linear system of the Gauss–Newton algorithm to update the current minimizer. When GN convergence is reached, this minimizer will be used as a starting point for the next QS minimization. Line 24 updates the ensemble. Line 25 propagates the updated ensemble to the next assimilation cycle. Figure 5 illustrates the QS scheme on a single analysis, as described in the caption.

To cycle the scheme, the DAW is then shifted with a small $S$. This ensures minimal cost function deformation, since few vectors of observation enter and few leave the DAW. This new cost function is then minimized using the forecast of the preceding minimizer as a starting point.

To keep this cost function deformation statistically consistent, Bocquet and Sakov (2014) advocate the use of $S = 1$ and assimilate one observation vector at the end of the DAW, which avoids multiple assimilation of the same observations. It is also the easiest way to ensure inequality (Eq. 35a). However, this method has been shown to be suboptimal because of the frequent Gaussian background approximations. Moreover, it ultimately fails to be QS since there is only one distant observation vector per analysis.

An alternative is to keep $S = 1$ but use all observations in the DAW and, consequently, relax the condition $K = L-S+1$. This way, observations are assimilated several times. This is done in the MDA IEnKS (Bocquet and Sakov, 2014). To keep the statistics consistent, at least in the linear/Gaussian case, the observations error covariances should be adequately altered. Hence, the MDA IEnKS is truly a QS scheme and makes a good reference scheme for our numerical experiments. However, the multiple assimilation of the observations introduces spurious correlations in the nonlinear/non-Gaussian case, which entail sub-optimality. A scheme similar to the IEnKSQS has also been successfully used in Carrassi et al. (2017) to compute model evidence. Indeed, the efficient computation of model evidence as an integral over the state space depended on the proper identification of a global
<figure><img src="image_10.png" /><figcaption>Figure 7. Schematics of an IEnKS QS minimization (a) and an IEnKS QC minimization (b). The rectangle contains the observations to be assimilated and the black dots represent the observations used in the current minimization. The “Gauss–Newton” surrounded by arrows represents the iterations of the Gauss–Newton procedure. The number of quasi-static steps is $N_Q = 3$. The flow of observations is controlled by the parameters ($L_0, L_1, L_2$) = (3, 5, 6). For the QC IEnKS, the number of GN iterations is controlled by the ($j_0, j_1, j_2$) parameters.</figcaption></figure>

maximum of the integrand. However, its implementation was based on the update of the ensemble whenever an observation batch is added to the cost function, which is not as numerically efficient as the scheme presented here.

The success of the QS minimization lies in the fact that, when an observation is successfully assimilated, the eMSE is reduced. Thus, the analysis probability mass concentrates around the true state. The analysis is then more likely to be in a neighborhood of the true state where the model is linear. The cost function non-quadraticity can then be increased by adding a new term in it. This is confirmed by the following argument: let $P^{(q)}$ be the IEnKS asymptotic eMSE at the qth step of a QS minimization. With the notations and assumptions of Sect. 2.3, i.e., in a linear context, we have the following recurrence relation:

$$ (P^{(q+1)})^{-1} = (P^{(q)})^{-1} + \sum_{L_q+1}^{L_{q+1}}, \qquad (39a) $$

$$ P^{(-1)} = \alpha^{2S} \frac{\alpha^2 - 1}{\alpha^2 (L_{N_Q-1} + 1)}. \qquad (39b) $$

Thus, $P^{(q+1)} < P^{(q)}$ and we can increase the cost function non-quadraticity by adding new terms in it as long as the

propagation of errors does not exceed the bound:

$$ \alpha^{2L_{q+1}} P^{(q)} \leq B. \qquad (40) $$

This implies

$$ L_{q+1} \leq L_q + S_{\max}, \qquad (41a) $$

$$ L_0 \leq L_{N_Q-1} + S_{\max} - S, \qquad (41b) $$

which yields $S \leq N_Q S_{\max}$. Therefore, the QS minimizations allow for a $N_Q$ times longer DAW.

Unfortunately, these QS minimizations are very expensive. Indeed, they add a third outer loop repeating $N_Q$ GN minimizations. The GN iterations used to compute the intermediate starting points give unnecessary precision; all that is required for these starting points is to be in a neighborhood of $x_{\infty S}$ where the model is almost linear. Thus, one can restrain the number of intermediate GN loops and save the full convergence to the last minimization. This is done in the quasi-convergent IEnKS (IEnKS QC) with the parameters ($j_q$)$_{q<N_Q}$. These parameters correspond to the numbers of GN loops in the intermediate QS minimizations. They are typically equal to 1 except for the last one. Algorithm 2 gives the pseudo-code of the IEnKS QC and the right panel in Fig. 7 is a schematic for it.
**Algorithm 1 One cycle of the IEnKS<sub>QS</sub>**

Require: $\mathbf{E}_0^b$ the background ensemble at $t_0$; $\lambda$ the inflation; $(L_q)_{q<N_Q}$ a list of DAW time indexes; $\epsilon$ the finite differences step; $\delta, j_{max}$ GN end of loop parameters; $\mathbf{1} = (1, \dots, 1)^T \in \mathbb{R}^n$ and $\mathbf{I}$ is the identity of $\mathbb{R}^n$.

$$1: \quad \mathbf{w}_0^a = \mathbf{0}$$

$$2: \quad \bar{\mathbf{x}}_0^b = \mathbf{E}_0^b \mathbf{1}/n$$

$$3: \quad \mathbf{X}_0^b = \lambda (\mathbf{E}_0^b - \bar{\mathbf{x}}_0^b \mathbf{1}^T)$$

$$4: \quad \text{for } q = 0, \dots, N_Q - 1 \text{ do}$$

$$5: \quad j = 0$$

6: repeat

$$7: \quad \bar{\mathbf{x}}_0^a = \bar{\mathbf{x}}_0^b + \mathbf{X}_0^b \mathbf{w}_0^a$$

$$8: \quad \mathbf{E}_0 = \bar{\mathbf{x}}_0^a \mathbf{1}^T + \epsilon \mathbf{X}_0^b$$

$$9: \quad \nabla J = (n - 1) \mathbf{w}_0^a$$

$$10: \quad \tilde{\nabla}^2 J = (n - 1) \mathbf{I}$$

$$11: \quad \mathbf{E}_K = \mathcal{M}^K (\mathbf{E}_0)$$

$$12: \quad \text{for } l = K, \dots, L_q \text{ do}$$

$$13: \quad \bar{y}_l = \mathcal{H}(\mathbf{E}_l) \mathbf{1}/n$$

$$14: \quad \mathbf{Y}_l = (\mathcal{H}(\mathbf{E}_l) - \bar{\mathbf{y}}_l \mathbf{1}^T) / \epsilon$$

$$15: \quad \nabla J = \nabla J - \mathbf{Y}_l^T \mathbf{R}^{-1} (\mathbf{y}_l - \bar{\mathbf{y}}_l)$$

$$16: \quad \tilde{\nabla}^2 J = \tilde{\nabla}^2 J - \mathbf{Y}_l^T \mathbf{R}^{-1} \mathbf{Y}_l$$

$$17: \quad \mathbf{E}_{l+1} = \mathcal{M}(\mathbf{E}_l)$$

18: end for

$$19: \quad \text{solve } \tilde{\nabla}^2 J \delta \mathbf{w} = \nabla J$$

$$20: \quad \mathbf{w}_0^a = \mathbf{w}_0^a - \delta \mathbf{w}$$

$$21: \quad j = j + 1$$

$$22: \quad \text{until } ||\delta \mathbf{w}|| \le \delta \text{ or } j \ge j_{\max}$$

23: end for

$$24: \quad \mathbf{E}_0^a = \mathbf{x}_0^a \mathbf{1}^T + \sqrt{n-1} \mathbf{X}_0^b \tilde{\nabla}^2 J^{-1/2}$$

$$25: \quad \mathbf{E}_S^b = \mathcal{M}^S (\mathbf{E}_0^a)$$

## 4 Numerical experiments with low-order models

In the following, we perform numerical experiments with the Lorenz 1963 (L63) and Lorenz 1995 (L95) models. L95 has already been presented in Sect. 2.4. L63 (Lorenz, 1963) is a simplified model for atmospheric convection. It is defined by the following ordinary differential equations:

$$\frac{dx}{dt} = \sigma (y - x), \qquad (42a)$$

$$\frac{dy}{dt} = \rho x - y - xz, \qquad (42b)$$

$$\frac{dz}{dt} = xy - \beta z. \qquad (42c)$$

These equations are integrated using a fourth-order Runge-Kutta scheme with a time step of $\delta t = 0.01$ and $(\sigma, \rho, \beta) = (10, 28, 8/3)$. The dynamics of the L63 model are chaotic, with a largest Lyapunov exponent given by $\lambda \approx 0.91$.

Both models are assumed perfect. The truth run is generated from a random state space point. The initial ensemble is generated from the truth with $\mathbf{B} = \mathbf{I}_m$ where $m = 40, 3$ and $n = 20, 3$ for L95 and L63, respectively. Observation vectors are generated from the truth with $\mathcal{H} = \mathbf{R} = \mathbf{I}_m$ every $\Delta t = 0.05$ for L95 and every $\Delta t = 0.02$ for L63. A burn-in period of $5 \times 10^3 \times \Delta t$ is enforced in both cases.

The IEnKS parameters are $\epsilon = 10^{-4}$, $\delta = 10^{-3}$, $j_{max} = 20$, $N_Q = 1$, $L_0 = L$. For the IEnKS<sub>QS</sub>, the QS parameters are $(L_0, L_1, \dots, L_{N_Q-1}) = K + (0, 1, \dots, N_Q - 1) \times \frac{s-1}{N_Q-1}$. For the IEnKS<sub>QC</sub>, the QS parameters are the same and, in addition, $(j_0, \dots, j_{N_Q-2}, j_{N_Q-1}) = (1, \dots, 1, 20)$. Sampling errors are systematically accounted for using the IEnKS finite-size version (Bocquet, 2011; Bocquet and Sakov, 2012; Bocquet et al., 2015) which avoids the need for inflation and its costly tuning. Finally the aRMSE is averaged over a number of cycles which is determined by the number of observations assimilated. We use $5 \times 10^5$ observation vectors for L95 and $5 \times 10^6$ observation vectors for L63.

Unlike Goodliff et al. (2015), our numerical experiments neither address increasing nonlinearity, nor do they address the use of climatological background error covariance matrices. Instead, we focus exclusively on the IEnKS performance dependence on the DAW key parameters $L, S$ and the number $N_Q$ of QS minimizations. Goodliff et al. (2015) numerically evaluates the QSVA approach with hybrid and EnVar techniques, similarly to Bocquet and Sakov (2013), and confirms the findings of Pires et al. (1996) albeit in a EnVar context. In terms of accuracy, Goodliff et al. (2015) show that the ensemble transform Kalman smoother (ETKS) outperforms all hybrid schemes in their numerical experiment. Since the IEnKS systematically outperforms the ETKS in all
<caption>**Algorithm 2** One cycle of the IEnKS<sub>QC</sub></caption>

**Require:** $\mathbf{E}_0^b$ the background ensemble at $t_0$; $\lambda$ the inflation; $(L_q)_{q<N_Q}$ a list of DAW time indexes; $(j_q)_{q<N_Q}$ the number of intermediate GN loops; $\epsilon$ the finite differences step; $\delta$ GN end of loop parameter; $\mathbf{1} = (1, \dots, 1)^T \in \mathbb{R}^n$ and $\mathbf{I}$ is the identity of $\mathbb{R}^n$;

1: $\mathbf{w}_0^a = \mathbf{0}$
2: $\bar{\mathbf{x}}_0^b = \mathbf{E}_0^b \mathbf{1}/n$
3: $\mathbf{X}_0^b = \lambda (\mathbf{E}_0^b - \bar{\mathbf{x}}_0^b \mathbf{1}^T)$
4: **for** $q = 0 \dots N_Q - 1$ **do**
5:     $j = 0$
6:     **repeat**
7:         $\bar{\mathbf{x}}_0^a = \bar{\mathbf{x}}_0^b + \mathbf{X}_0^b \mathbf{w}_0^a$
8:         $\mathbf{E}_0 = \bar{\mathbf{x}}_0^a \mathbf{1}^T + \epsilon \mathbf{X}_0^b$
9:         $\nabla J = (n - 1) \mathbf{w}_0^a$
10:        $\tilde{\nabla}^2 J = (n - 1) \mathbf{I}$
11:        $\mathbf{E}_K = \mathcal{M}^K (\mathbf{E}_0)$
12:        **for** $l = K, \dots, L_q$ **do**
13:            $\bar{\mathbf{y}}_l = \mathcal{H}(\mathbf{E}_l) \mathbf{1}/n$
14:            $\mathbf{Y}_l = (\mathcal{H}(\mathbf{E}_l) - \bar{\mathbf{y}}_l \mathbf{1}^T)/\epsilon$
15:            $\nabla J = \nabla J - \mathbf{Y}_l^T \mathbf{R}^{-1} (\mathbf{y}_l - \bar{\mathbf{y}}_l)$
16:            $\tilde{\nabla}^2 J = \tilde{\nabla}^2 J - \mathbf{Y}_l^T \mathbf{R}^{-1} \mathbf{Y}_l$
17:            $\mathbf{E}_{l+1} = \mathcal{M}(\mathbf{E}_l)$
18:        **end for**
19:        **solve** $\tilde{\nabla}^2 J \delta \mathbf{w} = \nabla J$
20:        $\mathbf{w}_0^a = \mathbf{w}_0^a - \delta \mathbf{w}$
21:        $j = j + 1$
22:    **until** $\|\delta \mathbf{w}\| \le \delta$ or $j \ge j_q$
23:    **end for**
24:    $\mathbf{E}_0^a = \mathbf{x}_0^a \mathbf{1}^T + \sqrt{n-1} \mathbf{X}_0^b \tilde{\nabla}^2 J^{-1/2}$
25:    $\mathbf{E}_S^b = \mathcal{M}^S (\mathbf{E}_0^a)$

conditions (Bocquet and Sakov, 2013, 2014) as long as the DAW length is not excessively long (for a chaotic model), then one concludes that our RMSEs would be systematically equal to or smaller than those reported for any hybrid scheme in Goodliff et al. (2015).

Figures 8 and 9 show the aRMSE of the IEnKS and IEnKS<sub>QS</sub> for both L95 and L63 as a function of S and L. The smoothing and filtering performance of the IEnKS increases for small values of L, S then decreases for high values of L, S. This is due to the appearances of local minima. As noted by Goodliff et al. (2015) with L63, the QS variant allows one to reach much longer DAWs, and improves the performance. However, some limits to the L95 model method are visible with when, for L = 50, best smoothing aRMSEs are reached for S < 50. However, the IEnKS<sub>QS</sub> filtering performance is invariant with L and improves with S as in the 4D-Var filtering performance shown in Fig. 3. As suggested by Pires et al. (1996), one may be tempted to estimate the useful DAW length from past observations beyond which the performance gain is negligible. However, their study estimated useful DAW length in a one-cycle 4D-Var context with a focus on the filtering RMSE. Hence, this estimation is not directly relevant for a cycled IEnKS<sub>QS</sub>. By contrast, in this study, a lot of observations have already been assimilated

and condensed in the background approximation. Thus, the performance gain with the DAW length comes from the precision of this Gaussian background approximation; a precision that the linearized theory is not able to provide, attested to by the filtering ($l = L$) performance independence of Eq. (26) with the DAW parameters.

Figure 10 compares the smoothing aRMSE (first column), the filtering aRMSE (second column) and the number of ensemble propagations (third column) of the IEnKS<sub>QS</sub> ($N_Q = S, S = L$), the IEnKS ($S = L$) and the IEnKS-MDA ($S = 1$), for both L95 and L63 as a function of L. The number of ensemble propagations is the total number of ensemble propagations in units of $\Delta t$ divided by the total number of assimilated observation vectors. For $L < 20$, all three algorithms show smoothing and filtering performance improvements with L. For $L > 20$, the IEnKS filtering and smoothing RMSEs increase because of the multiple local extrema. For $L < 40$, the IEnKS<sub>QS</sub> has smaller aRMSEs than the IEnKS-MDA. Because the IEnKS<sub>QS</sub> is SDA by design, it does not suffer from sub-optimality related to multiple assimilations and nonlinearity. Moreover, the IEnKS<sub>QS</sub> always requires less propagations of the ensemble, which improves the computational cost. However, for $L > 40$ with the L95 model, the quasi-static approach cannot sustain the nonlinearity any-
<figure><img src="image_11.png" /><figcaption>**Figure 8.** IEnKS (lower triangles) and IEnKS<sub>QS</sub> (upper triangles, $N_Q = S$) smoothing and filtering aRMSEs as a function of *L* and *S* with the L95 model. The *L* parameter is on the abscissa axis for the IEnKS and on the ordinate axis for the IEnKS<sub>QS</sub>. The *S* parameter is on the abscissa axis for the IEnKS<sub>QS</sub> and on the ordinate axis for the IEnKS. For readability, smoothing RMSE beyond 0.10 and filtering RMSE beyond 0.20 are in the same color and the scale is logarithmic.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** IEnKS (lower triangles) and IEnKS<sub>QS</sub> (upper triangles, $N_Q = S$) smoothing and filtering aRMSEs as a function of *L* and *S* with the L63 model. The *L* parameter is on the abscissa axis for the IEnKS and on the ordinate axis for the IEnKS<sub>QS</sub>. The *S* parameter is on the abscissa axis for the IEnKS<sub>QS</sub> and on the ordinate axis for the IEnKS. For readability, smoothing and filtering RMSE beyond 0.10 are in the same color and the scale is logarithmic.</figcaption></figure>

more and the IEnKS<sub>QS</sub> aRMSE degrades. Hence, the IEnKS-MDA *S* = 1 still has the best performance. With the L63 model, the IEnKS<sub>QS</sub> is always better than the IEnKS-MDA suggesting that *L* could still be increased.

Figure 11 compares the smoothing aRMSE (first column), the filtering aRMSE (second column) and the number of ensemble propagations (third column) of the IEnKS<sub>QS</sub> and IEnKS<sub>QC</sub> as a function of the $N_Q$ parameter. The IEnKS<sub>QS</sub> aRMSE quickly decreases after a point for both algorithms

and for both models. Before this point, the algorithms fail to find the global minimum and the RMSE is close to the climatological variance. After this point, the algorithms succeed in finding the global minimum and the RMSE is low. For the IEnKS<sub>QS</sub> with the L95 model, this point can be estimated using results of Sect. 3 by $S/S_{max} \approx 3.6$, in remarkable agreement with Fig. 11. This point comes later for the IEnKS<sub>QC</sub> but it demands less ensemble propagations making this algorithm numerically more efficient than the IEnKS<sub>QS</sub>. However
<figure><img src="image_13.png" /><figcaption>**Figure 10.** IEnKS-MDA ($S = 1$), IEnKS<sub>QS</sub> and IEnKS ($S = L, N_Q = S$) filtering and smoothing aRMSEs and number of ensemble propagations as a function of $L$ for the L95 and L63 models (logarithmic scale). The ensemble propagations are in units of $\Delta t$. For instance, with $L = 50$, a single ensemble propagation through the DAW counts for 50.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** IEnKS<sub>QS</sub> and IEnKS<sub>QC</sub> ($S = L = 50$) filtering and smoothing aRMSEs and number of ensemble propagations as a function of $N_Q$ for the L95 and L63 models (logarithmic scale). The ensemble propagations are in units of $\Delta t$.</figcaption></figure>
those ensemble propagations have different behavior for the L95 and L63 models when the minimizations fail to find the global minimum. For the L95 model, the number of ensemble propagations is high meaning that the minimization takes a lot of iterations and fails to converge. For the L63 model, the number of ensemble propagations is low indicating that minimizations converge but to a non-global extrema.

# 5 Conclusions

In this paper, we have extended the study of Pires et al. (1996) on quasi-static variational data assimilation, focused on 4D-Var technique, to cycled data assimilation schemes and specifically four-dimensional nonlinear ensemble variational techniques, an exemplar of this being the iterative ensemble Kalman smoother (IEnKS).

The long term impact of cycling was first theoretically investigated in a linear context for 4D-Var and the IEnKS, then numerically for the IEnKS in a nonlinear context. The way information is propagated between data assimilation cycles indeed explains the difference between 4D-Var and the IEnKS. Both reveal performance improvements with the DAW parameter *S*, which counts the number of observation vectors within the DAW, as well as the time shift between cycles. This is a consequence of the Gaussian background approximation; the larger *S* is, the less the assimilation relies on it.

However, it is observed that this improvement has a limit in the chaotic, perfect model case. The cost function global minimum basin of attraction appears to shrink with increasing *L*. This causes the Gauss-Newton procedure to miss the cost function global minimum, which deteriorates the assimilation performance.

Quasi-static minimizations lead slowly but surely to the global minimum by repeated cost function minimizations. As the DAW length *L* is gradually increased, the starting point of the minimization remains in the global minimum basin of attraction. For most *S*, *L* couples, the quasi-static IEnKS turns out to be a more accurate substitute for the multiple data assimilation IEnKS (IEnKS-MDA).

Unfortunately, this method (IEnKS<sub>QS</sub>) adds an outer loop which could significantly increase the numerical cost. Precision on intermediate minima being superfluous, one can limit the intermediate Gauss-Newton number of loops. The unavoidable space increments required to minimize the non-quadratic cost function are thus reported in time in the *quasi-convergent* IEnKS (IEnKS<sub>QC</sub>).

We did not focus on the applicability of the methods to high-dimensional and imperfect models. In particular, we considered very long DAWs, which, even if of high mathematical interest or for low-order reliable models, is less relevant for significantly noisy models. However, we know from Swanson et al. (1998), that the perfect model results are expected to extend to the imperfect model case provided that the growth rate of the model error is similar to that of the leading Lyapunov vectors of the model. This is likely to also apply to a (strong-constraint) IEnKS<sub>QS</sub>. Finally, an extension to this work would therefore consist in investigating the same ideas but using a weak constraint 4D-Var and IEnKS (Trémolet, 2006; Sakov and Bocquet, 2018; Sakov et al., 2018).

*Data availability.* No data sets were used in this article.
## Appendix A: Performance of 4D-Var in the linear, univariate case

The objective of this Appendix is to establish a recurrence relation between the 4D-Var eMSE of each cycle. From this relation we will get an expression for the 4D-Var asymptotic eMSE.

We assume $m = 1$ and $\mathcal{M}^l(x) = \alpha^l x$. At the $k$th cycle, the gradient and Hessian of the 4D-Var cost function Eq. (7) are

$$ \nabla J(x_{kS}; y_{:}, x_{kS}^b) = -\frac{1}{b}(x_{kS}^b - x_{kS}) - \frac{h}{r}\sum_{l=K}^{L}\alpha^l(y_{kS+l} - h\alpha^l x_{kS}), \quad (A1) $$

$$ \nabla^2 J(x_{kS}; y_{:}, x_{kS}^b) = \frac{1}{b} + \frac{h^2}{r} \sum_{l=K}^{L} \alpha^{2l}, \quad (A2) $$

where $y_{kS+L:kS+K}$ is temporarily denoted $y_{:}$. Because of the operators' linearity, $J$ is quadratic with respect to $x_{kS}$ and convex. Hence, its minimizer $x_{kS}^a$ exists and is characterized by the null gradient equation:

$$ \nabla J(x_{kS}^a; y_{:}, x_{kS}^b) = 0. \quad (A3) $$

With an exact Taylor expansion around the state $x_{kS}$ we obtain the following:

$$ 0 = \nabla J(x_{kS}; y_{:}, x_{kS}^b) + \nabla^2 J(x_{kS}; y_{:}, x_{kS}^b) \times (x_{kS}^a - x_{kS}). \quad (A4) $$

Note that $\nabla^2 J(x_{kS}; y_{:}, x_{kS}^b)$, given by Eq. (A2), is not random and does not depend on $x_{kS}$, $y_{:}$, $x_{kS}^b$. That is why it is simply noted $\nabla^2 J$. Using Eqs. (1) and (A4), we have

$$ \nabla^2 J \cdot (x_{kS}^a - x_{kS}) = \frac{1}{b} (x_{kS}^b - x_{kS}) + \frac{h}{r} \sum_{l=K}^{L} \alpha^l \varepsilon_{kS+l}. \quad (A5) $$

The random variable $x_{kS}^b - x_{kS} = \alpha^S (x_{(k-1)S}^a - x_{(k-1)S})$ is independent from the errors $\varepsilon_{kS+K}, \dots, \varepsilon_{kS+L}$. Thus, taking the expectation of the square of Eq. (A5) gives the following expression for the eMSE $P_{kS}^{4D-Var}$ of 4D-Var at time $t_{kS}$:

$$ P_{kS}^{4D-Var} = \mathbb{E} \left[ (x_{kS} - x_{kS}^a)^2 \right], \\ = \left( \nabla^2 J \right)^{-2} \left( \frac{\alpha^{2S}}{b^2} P_{(k-1)S}^{4D-Var} + \frac{h^2}{r} \sum_{l=K}^{L} \alpha^{2l} \right), \quad (A6a) $$

$$ P_{-S}^{4D-Var} = b. \quad (A6b) $$

Introducing the notation

$$ \Sigma_K^L = \frac{h^2}{r} \sum_{l=K}^{L} \alpha^{2l} = \frac{h^2}{r} \alpha^{2K} \frac{\alpha^{2S} - 1}{\alpha^2 - 1}, \quad (A7) $$

$$ \Delta = \frac{\alpha^{2S}}{(1 + b \Sigma_K^L)^2}, \quad (A8) $$

we obtain

$$ P_{kS}^{4D-Var} = \Delta P_{(k-1)S}^{4D-Var} + \frac{\Sigma_K^L}{\left(\frac{1}{b} + \Sigma_K^L\right)^2}. \quad (A9) $$

Thus, $(P_{kS}^{4D-Var})_k$ is an arithmetico-geometric sequence. Its limit $P_{\infty S}^{4D-Var}$ depends on the value of $\Delta$ in the following way:

$$ P_{\infty S}^{4D-Var} = \begin{cases} \infty & \text{if } \Delta \ge 1, \\ \frac{b^2 \sum_K^L \frac{\Delta}{1-\Delta}}{\alpha^{2S}} & \text{otherwise.} \end{cases} \quad (A10) $$

The generalization to any asymptotic eMSE with lag $L-l$ is straightforward:

$$ P_{\infty S+l}^{4D-Var} = \begin{cases} \infty & \text{if } \Delta \ge 1, \\ \frac{b^2 \sum_K^L \frac{\Delta}{1-\Delta}}{\alpha^{2(S-l)}} & \text{otherwise.} \end{cases} \quad (A11) $$

In the multivariate, diagonal case the algebra can be conducted on each direction independently. The eMSE in this case is the sum of the univariate eMSEs of each direction.

## Appendix B: Performance of the IEnKS in the linear, univariate case

The objective of this appendix is to establish a recurrence relation between the IEnKS eMSE of each cycle. From this relation we will get an expression for the IEnKS asymptotic eMSE.

First, it is proven by recurrence that for all $k \ge 0$, $G(x_{kS}|y_{kS+L:K})$ is Gaussian with moments

$$ \mathbb{E}[x_{kS} | y_{kS+L:K}] = \bar{x}_{kS}^a, \quad (B1) $$

$$ \mathbb{V}[x_{kS} | y_{kS+L:K}] = (X_{kS}^a)^2, \quad (B2) $$

where $\mathbb{V}$ is the variance of a random variable and $\bar{x}_{kS}^a$, $X_{kS}^a$ are defined by Eqs. (15) and (16). Because of the assumptions Eqs. (25c) and (25d) with Eq. (13) one gets

$$ G(x_0 | y_{L:K}) = J(w_0; y_{L:K}, \mathbf{E}_0^b) + c_0, \quad (B3) $$

where $c_0$ is a constant independent from $x_0$ and $w_0$. Hence, $G(x_0|y_{L:K})$ is Gaussian and its moments are given by Eqs. (B1) and (B2) with $k=0$. Now, assume $G(x_{kS}|y_{kS+L:K})$ is Gaussian with moments given by
Eqs. (B1) and (B2) for a $k \ge 0$. Because $x_{(k+1)S} = M^S(x_{kS})$ with $M$ affine, one gets

$$ \mathbb{E}[x_{(k+1)S}|y_{kS+L:K}] = M^S(\bar{x}_{kS}^a) = \bar{x}_{(k+1)S}^b, \quad (B4) $$

$$ \mathbb{V}[x_{(k+1)S}|y_{kS+L:K}] = \alpha^{2S}(X_{kS}^a)^2 = X_{(k+1)S}^b, \quad (B5) $$

using the conditional expectation properties and Eqs. (25a) and (25b). This result together with Eq. (6) and assumption Eq. (13) yields

$$ \begin{aligned} G(x_{(k+1)S}|y_{(k+1)S+L:K}) &= J(w_{(k+1)S}; y_{(k+1)S+L:(k+1)S+K}, \mathbf{E}_{(k+1)S}^{\mathbf{b}}) + c_{k+1}, \end{aligned} \quad (B6) $$

where $c_{k+1}$ is a constant independent from $x_{(k+1)S}$ and $w_{(k+1)S}$. Hence $G(x_{(k+1)S}|y_{(k+1)S+L:K})$ is Gaussian with moments given by Eqs. (B1) and (B2).

The conditional variance Eq. (B2) is therefore related to the IEnKS performance by the total law of expectation:

$$ \begin{aligned} P_{kS}^{\text{IEnKS}} &= \mathbb{E}[(\bar{x}_{kS}^a - x_{kS})^2], \\ &= \mathbb{E}[\mathbb{E}[(\bar{x}_{kS}^a - x_{kS})^2|y_{kS+L:K}]], \\ &= \mathbb{E}[\mathbb{V}[x_{kS}|y_{kS+L:K}]], \\ &= \mathbb{E}[(X_{kS}^a)^2]. \end{aligned} \quad (B7) $$

Then, from Eq. (16) we get the recurrence relation:

$$ (X_{(k+1)S}^a)^{-2} = \alpha^{-2S} (X_{kS}^a)^{-2} + \frac{h^2}{r} \sum_{l=K}^{L} \alpha^{2l}, \quad (B8a) $$

$$ (X_0^a)^{-2} = \frac{1}{b} + \frac{h^2}{r} \sum_{l=K}^{L} \alpha^{2l}. \quad (B8b) $$

Thus, $X_{kS}^a$ is not random and $P_{kS}^{\text{IEnKS}} = (X_{kS}^a)^2$. Equation (B8) shows that the sequence of inverse IEnKS eMSEs is arithmetico-geometric:

$$ (P_{kS}^{\text{IEnKS}})^{-1} = \alpha^{-2S} (P_{(k-1)S}^{\text{IEnKS}})^{-1} + \Sigma_K^L, \quad (B9a) $$

$$ (P_{-S}^{\text{IEnKS}})^{-1} = b^{-1}, \quad (B9b) $$

where the notation Eq. (A7) has been used. Properties of arithmetico-geometric sequences allow one to obtain the IEnKS asymptotic eMSE:

$$ P_{\infty S}^{\text{IEnKS}} = \begin{cases} 0 & \text{if } |\alpha| \le 1, \\ \frac{r}{h^{2\alpha^2 L}} \frac{\alpha^2 - 1}{\alpha^2} & \text{otherwise,} \end{cases} \quad (B10) $$

and the generalization to asymptotic eMSEs with lag $L - l$ is straightforward:

$$ P_{\infty S+l}^{\text{IEnKS}} = \begin{cases} 0 & \text{if } |\alpha| \le 1, \\ \frac{r}{h^{2\alpha^2(L-l)}} \frac{\alpha^2 - 1}{\alpha^2} & \text{otherwise.} \end{cases} \quad (B11) $$

Let us now show that the IEnKS eMSE is optimal. Let $x_{kS}^a$ ($y_{kS+L:K}$) be the 4D-Var analysis or any other function of $y_{kS+L:K}$. A bias-variance decomposition (e.g., Bishop, 2006) of this estimator yields

$$ \begin{aligned} \mathbb{E}_{x_{kS}, y_{kS+L:K}} [(x_{kS} - x_{kS}^a)^2] = \mathbb{E}_{y_{kS+L:K}} [\mathbb{V}_{x_{kS}} [x_{kS}|y_{kS+L:K}]] \\ + \mathbb{E}_{y_{kS+L:K}} [(\mathbb{E}_{x_{kS}} [x_{kS}|y_{kS+L:K}] - x_{kS}^a)^2]. \end{aligned} \quad (B12) $$

Replacing the moments with Eqs. (B1) and (B2) yields

$$ P_{kS}^{\text{4D-Var}} = P_{kS}^{\text{IEnKS}} + \mathbb{E}[(\bar{x}_{kS}^a - x_{kS}^a)^2]. \ge P_{kS}^{\text{IEnKS}}. \quad (B13) $$

In the multivariate, diagonal case the algebra can be conducted on each direction independently. Thus, the eMSE in this case is the sum of the univariate eMSEs of each direction.

## Appendix C: Expression of the averaged cost function

The IEnKS averaged cost function $J_{\infty S}$ is the $N$ goes to $\infty$ limit of

$$ \frac{1}{N} \sum_{k=0}^{N-1} J(\mathbf{w}; y_{kS+L:kS+K}, \mathbf{E}_{kS}^{\mathbf{b}}) = \frac{1}{2} \|\mathbf{w}\|^2 \quad (C1) \\ + \frac{1}{2N} \sum_{l=K}^{L} \sum_{k=0}^{N-1} \| y_{kS+l} - \mathcal{H} \circ \mathcal{M}^l (\bar{x}_{kS}^b + \mathbf{X}_{kS}^b \mathbf{w}) \|^2_{\mathbf{R}^{-1}}. $$

Expanding the squared norm around $\mathcal{H} \circ \mathcal{M}^l(\mathbf{x}_{kS})$ using Eq. (1) gives

$$ \begin{aligned} & \frac{1}{2} \| \mathbf{y}_{kS+l} - \mathcal{H} \circ \mathcal{M}^l (\bar{\mathbf{x}}_{kS}^b + \mathbf{X}_{kS}^b \mathbf{w}) \|^2_{\mathbf{R}^{-1}} \\ & = \frac{1}{2} \| \boldsymbol{\varepsilon}_{kS+l} \|^2_{\mathbf{R}^{-1}} + \frac{1}{2} \| \delta \mathbf{y}_{kS+l} \|^2_{\mathbf{R}^{-1}} + \boldsymbol{\varepsilon}_{kS+l}^{\mathrm{T}} \mathbf{R}^{-1} \delta \mathbf{y}_{kS+l}, \end{aligned} \quad (C2) $$

where $\delta \mathbf{y}_{kS+l} = \mathcal{H} \circ \mathcal{M}^l(\mathbf{x}_{kS}) - \mathcal{H} \circ \mathcal{M}^l(\bar{\mathbf{x}}_{kS}^b + \mathbf{X}_{kS}^b \mathbf{w})$. We assume that random variables exist ($\boldsymbol{\varepsilon}_{\infty S}, \mathbf{x}_{\infty S}, \mathbf{E}_{\infty S}^a$) whose distribution is invariant and ergodic with respect to the shift transformation:

$$ T : (\boldsymbol{\varepsilon}_{kS}, \mathbf{x}_{kS}, \mathbf{E}_{kS}^{\mathbf{a}}) \mapsto (\boldsymbol{\varepsilon}_{(k+1)S}, \mathbf{x}_{(k+1)S}, \mathbf{E}_{(k+1)S}^{\mathbf{a}}). \quad (C3) $$

Then because the $(\boldsymbol{\varepsilon}_{kS})_k$ are mutually independent, independent from the $(\mathbf{x}_{kS}, \mathbf{E}_{kS}^a)_k$ and identically distributed, $p(\boldsymbol{\varepsilon}_{\infty S}, \mathbf{x}_{\infty S}, \mathbf{E}_{\infty S}^a) = p(\boldsymbol{\varepsilon}_0) p(\mathbf{x}_{\infty S}, \mathbf{E}_{\infty S}^a)$. By Birkhoff's ergodic theorem (see Walters, 1982) we get

$$ \frac{1}{N} \sum_{k=0}^{N-1} \| \boldsymbol{\varepsilon}_{kS+l} \|^2_{\mathbf{R}^{-1}} \xrightarrow{N \to \infty} \mathbb{E} \left[ \| \boldsymbol{\varepsilon}_0 \|^2_{\mathbf{R}^{-1}} \right] = d, \quad (C4) $$

$$ \frac{1}{N} \sum_{k=0}^{N-1} \boldsymbol{\varepsilon}_{kS+l}^{\mathrm{T}} \delta \mathbf{y}_{kS+l} \xrightarrow{N \to \infty} \mathbb{E}[\boldsymbol{\varepsilon}_0^{\mathrm{T}}] \mathbb{E}[\delta \mathbf{y}_{\infty S+l}] = 0, \quad (C5) $$

$$ \frac{1}{N} \sum_{k=0}^{N-1} \| \delta \mathbf{y}_{kS+l} \|^2_{\mathbf{R}^{-1}} \xrightarrow{N \to \infty} \mathbb{E} \left[ \| \delta \mathbf{y}_{\infty S+l} \|^2_{\mathbf{R}^{-1}} \right], \quad (C6) $$
where $\delta \mathbf{y}_{\infty S+l} = \mathcal{H} \circ \mathcal{M}^l(\mathbf{x}_{\infty S}) - \mathcal{H} \circ \mathcal{M}^l(\bar{\mathbf{x}}_{\infty S}^b + \mathbf{X}_{\infty S}^b \mathbf{w})$
and $\bar{\mathbf{x}}_{\infty S}^b, \mathbf{X}_{\infty S}^b$ are respectively the mean and normalized
anomaly of $\mathcal{M}^S(\mathbf{E}_{\infty S}^a)$. Finally,

$$
J_{\infty S}(\mathbf{w}) = \frac{1}{2} \|\mathbf{w}\|^2 + \frac{\mathrm{d}S}{2} + \frac{1}{2} \sum_{l=K}^{L} \mathbb{E} \left[ \| \delta \mathbf{y}_{\infty S+l} \|_{\mathbf{R}^{-1}}^2 \right]. \quad (C7)
$$
*Competing interests.* The authors declare that they have no conflict of interest.

*Special issue statement.* This article is part of the special issue “Numerical modeling, predictability and data assimilation in weather, ocean and climate: A special issue honoring the legacy of Anna Trevisan (1946–2016)”. It is not associated with a conference.

*Acknowledgements.* The authors are grateful to Selime Gürol for her comments and suggestions on the manuscript. The authors are also grateful to two reviewers, Carlos Pires and an anonymous reviewer, for raising very interesting questions and for their suggestions, and to Alberto Carrassi, acting as editor, for his remarks and suggestions. CEREA is a member of Institut Pierre-Simon Laplace (IPSL).

Edited by: Alberto Carrassi

Reviewed by: Carlos Pires and one anonymous referee

## References

Asch, M., Bocquet, M., and Nodet, M.: Data assimilation: methods, algorithms, and applications, Fundamentals of Algorithms, Society for Industrial and Applied Mathematics, Philadelphia, USA, 306 pp., 2016.

Bishop, C.: Pattern Recognition and Machine Learning, Information Science and Statistics, Springer-Verlag, New York, USA, 738 pp., 2006.

Björck, Å.: Numerical methods for least squares problems, Society for Industrial and Applied Mathematics, Philadelphia, USA, 408 pp., https://doi.org/10.1137/1.9781611971484, 1996.

Bocquet, M.: Ensemble Kalman filtering without the intrinsic need for inflation, Nonlin. Processes Geophys., 18, 735–750, https://doi.org/10.5194/npg-18-735-2011, 2011.

Bocquet, M.: Localization and the iterative ensemble Kalman smoother, Q. J. Roy. Meteor. Soc., 142, 1075–1089, https://doi.org/10.1002/qj.2711, 2016.

Bocquet, M. and Carrassi, A.: Four-dimensional ensemble variational data assimilation and the unstable subspace, Tellus A, 69, 1304504, https://doi.org/10.1080/16000870.2017.1304504, 2017.

Bocquet, M. and Sakov, P.: Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems, Nonlin. Processes Geophys., 19, 383–399, https://doi.org/10.5194/npg-19-383-2012, 2012.

Bocquet, M. and Sakov, P.: Joint state and parameter estimation with an iterative ensemble Kalman smoother, Nonlin. Processes Geophys., 20, 803–818, https://doi.org/10.5194/npg-20-803-2013, 2013.

Bocquet, M. and Sakov, P.: An iterative ensemble Kalman smoother, Q. J. Roy. Meteor. Soc., 140, 1521–1535, https://doi.org/10.1002/qj.2236, 2014.

Bocquet, M., Raanes, P. N., and Hannart, A.: Expanding the validity of the ensemble Kalman filter without the intrinsic need for inflation, Nonlin. Processes Geophys., 22, 645–662, https://doi.org/10.5194/npg-22-645-2015, 2015.

Carrassi, A., Bocquet, M., Hannart, A., and Ghil, M.: Estimating model evidence using data assimilation, Q. J. Roy. Meteor. Soc., 143, 866–880, https://doi.org/10.1002/qj.2972, 2017.

Cosme, E., Verron, J., Brasseur, P., Blum, J., and Auroux, D.: Smoothing Problems in a Bayesian Framework and Their Linear Gaussian Solutions, Mon. Weather Rev., 140, 683–695, https://doi.org/10.1175/MWR-D-10-05025.1, 2012.

Goodliff, M., Amezcua, J., and Van Leeuwen, P. J.: Comparing hybrid data assimilation methods on the Lorenz 1963 model with increasing non-linearity, Tellus A, 67, 26928, https://doi.org/10.3402/tellusa.v67.26928, 2015.

Judd, K., Smith, L., and Weisheimer, A.: Gradient free descent: shadowing, and state estimation using limited derivative information, Physica D, 190, 153–166, https://doi.org/10.1016/j.physd.2003.10.011, 2004.

Lorenz, A.: Four-dimensional variational data assimilation, in: Advanced Data Assimilation for Geosciences, Lecture Notes of the Les Houches School of Physics: Special Issue, June 2012, Oxford University Press, Les Houches, France, 31–73, 2014.

Lorenz, E.: Deterministic Nonperiodic Flow, J. Atmos. Sci., 20, 130–141, https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2, 1963.

Lorenz, E. and Emanuel, K.: Optimal Sites for Supplementary Weather Observations: Simulation with a Small Model, J. Atmos. Sci., 55, 399–414, https://doi.org/10.1175/1520-0469(1998)055<0399:OSFSWO>2.0.CO;2, 1998.

Miller, R. N., Ghil, M., and Gauthiez, F.: Advanced Data Assimilation in Strongly Nonlinear Dynamical Systems, J. Atmos. Sci., 51, 1037–1056, https://doi.org/10.1175/1520-0469(1994)051<1037:ADAISN>2.0.CO;2, 1994.

Pires, C., Vautard, R., and Talagrand, O.: On extending the limits of variational assimilation in nonlinear chaotic systems, Tellus A, 48, 96–121, https://doi.org/10.3402/tellusa.v48i1.11634, 1996.

Sakov, P. and Bocquet, M.: Asynchronous data assimilation with the EnKF in presence of additive model error, Tellus A, 70, 1414545, https://doi.org/10.1080/16000870.2017.1414545, 2018.

Sakov, P., Haussaire, J.-M., and Bocquet, M.: An iterative ensemble Kalman filter in presence of additive model error, Q. J. Roy. Meteor. Soc., https://doi.org/10.1002/qj.3213, accepted for publication, 2018.

Swanson, K., Vautard, R., and Pires, C.: Four-dimensional variational assimilation and predictability in a quasi-geostrophic model, Tellus A, 50, 369–390, https://doi.org/10.1034/j.1600-0870.1998.t01-4-00001.x, 1998.

Trémolet, Y.: Accounting for an imperfect model in 4D-Var, Q. J. Roy. Meteor. Soc., 132, 2483–2504, https://doi.org/10.1256/qj.05.224, 2006.

Trevisan, A., D'Isidoro, M., and Talagrand, O.: Four-dimensional variational assimilation in the unstable subspace and the optimal subspace dimension, Q. J. Roy. Meteor. Soc., 136, 487–496, https://doi.org/10.1002/qj.571, 2010.

Walters, P.: An Introduction to Ergodic Theory, Graduate texts in mathematics, Springer-Verlag, New York, USA, 250 pp., 1982.

Ye, J., Rey, D., Kadakia, N., Eldridge, M., Morone, U., Rozdeba, P., Abarbanel, H., and Quinn, J.: Systematic variational method for statistical nonlinear state and parameter estimation, Phys. Rev. E, 92, 052901, https://doi.org/10.1103/PhysRevE.92.052901, 2015.