Hence, an alternative to physics-based models is to use
available observations of the system and statistical ap-
proaches to discover equations and then make predictions.
This has been introduced in several papers using combina-
tions and polynomials of observed variables as well as sparse
regressions or model selection strategies (Brunton et al.,
2016; Rudy et al., 2017; Mangiarotti and Huc, 2019). Those
methods have then been extended to the case of noisy and ir-
regular observation sampling, using a Bayesian framework as
in data assimilation (Bocquet et al., 2019; North et al., 2022).
Alternatively, some authors used data assimilation and lo-
cal linear regressions based on analogs (Tandeo et al., 2015;
Lguensat et al., 2017) or iterative data assimilation coupled
with neural networks (Brajard et al., 2020; Fablet et al., 2021;
Brajard et al., 2021) to make data-driven predictions without
discovering equations.

However, many approaches cited above assume that the
full state of the system is observed, which is a strong as-
sumption. Indeed, in a lot of applications in geophysics, im-
portant components of the system are never or only partially
observed, such as the deep ocean (see, e.g., Jayne et al.,
2017), and data-driven methods fail to make good predic-
tions. To deal with those strong constraints, i.e., when the
model is unknown and when some components of the system
are never observed, combination of data assimilation and ma-
chine learning shows potential (see, e.g., Wikner et al., 2021).
Additionally, an option is to use time-delay embedding of the
available components of the system (Takens, 1981; Brunton
et al., 2017), whereas another option is to find latent repre-
sentations of the dynamical system (see, e.g., Talmon et al.,
2015; Ouala et al., 2020). In this study, we will show that
there are strong relationships between those two approaches.

Here, we propose a simple algorithm using linear and
Gaussian assumptions based on a state-space formulation.
This classic Bayesian framework, used in data assimilation,
is able to deal with a dynamical model (physics- or data-
driven) and observations (partial and noisy). Three main
ideas are used: (i) augmented state formulation (Kitagawa,
1998), (ii) global linear approximation of the dynamical sys-
tem (Korda and Mezić, 2018) and (iii) estimation of the dy-
namical parameters using an iterative algorithm combined
with Kalman recursions (Shumway and Stoffer, 1982). The
current paper is thus an extension of Shumway and Stoffer
(1982) to never-observed components of a dynamical system,
using a state-augmentation strategy. The proposed frame-
work is probabilistic, where the state of the system is approx-
imated using a Gaussian distribution (with a mean vector and
a covariance matrix). The algorithm is iterative, where a cat-
alog is updated at each iteration and used to learn a linear
dynamical model. The final estimate of this catalog corre-
sponds to a new system of variables, including latent ones.

The proposed methodology is based on an important as-
sumption: the surrogate model is linear. Although it can be
considered a disadvantage compared to nonlinear models,
this linear assumption also has interesting properties. Indeed,

nonlinear models combined with state augmentation are a
very broad family of models and may lead to identifiabil-
ity issues. Using linear dynamics already leads to a very
flexible family of models since the latent variable may de-
scribe nonlinearities and include, for example, any transfor-
mation of the observed or non-observed components of a dy-
namical model. Furthermore, it allows rigorous estimation of
the parameters using well-established statistical algorithms
which can be run at a low computational cost. The proposed
methodology is evaluated on a low-dimensional and weakly
nonlinear chaotic model. As this paper is a proof of concept,
a linear surrogate model is certainly well suited for this situ-
ation.

The paper is organized as follows. Firstly, the methodol-
ogy is explained in Sect. 2. Secondly, Sect. 3 describes the
experiment using the Lorenz-63 system. Thirdly, the results
are reported in Sect. 4. The conclusions and perspectives are
given in Sect. 5.

## 2 Methods

The methodology proposed in this paper is borrowed from
data assimilation, machine learning and dynamical systems.
It is summarized in Fig. 1 and explained below.

In data assimilation, the goal is to estimate, from partial
and noisy observations **y**, the full state of a system **x**. When
the dynamical model used to propagate **x** in time is avail-
able (i.e., when model equations are given), classic data as-
similation techniques are used to retrieve unobserved com-
ponents of the system. For instance, in the Lorenz-63 system
(Lorenz, 1963), if only two variables ($x_2$ and $x_3$ in the exam-
ple defined below) are observed, knowing the Lorenz equa-
tions (system of three ordinary differential equations), it is
possible to retrieve the unobserved one ($x_1$ in our example
below). However, this estimation requires good estimates of
model and observation error statistics (see, e.g., Dreano et al.,
2017; Pulido et al., 2018).

Now, if the model equations are not known and observa-
tions of the system are available over a sufficient period of
time, it is possible to use data-driven methods to mathemati-
cally approximate the system dynamics. In this paper, a linear
approximation is used to model the relationship of the state
vector **x** between two time steps. It is parameterized with
the matrix **M**, whose dimension is equal to the square of the
state space. Moreover, a linear observation operator is intro-
duced to relate the partial observations **y** and the state **x**. It
is written using a matrix **H**, with its dimension equal to the
observation-space times the state-space dimensions. Nonlin-
ear and adaptive operators and noisy observations could be
taken into account but, for the sake of simplicity, only the
linear and non-noisy case is considered in this paper.

Mathematically, matrices (M, H) and vectors (x, y) are
linked using a Gaussian and linear state-space model such
<figure><img src="image_42.png" /><figcaption>**Figure 1.** Schematic of the proposed methodology, illustrated using the Lorenz-63 system. The algorithm is initialized with a Gaussian random noise for the hidden component (i.e., $z_1$) and with partial observations of the system (i.e., $y_2$ and $y_3$). Then, an iterative procedure is applied with a linear regression, a covariance computation, the Kalman recursions and a random sampling. This algorithm iteratively maximizes the likelihood of the observations denoted $\mathcal{L}$. After convergence of the algorithm, a hidden component $z_1$ is stabilized and represented by a Gaussian distribution represented by the mean $\mathbf{x}_1^s$ and variance $P_1^s$.</figcaption></figure>

that

$$ \mathbf{x}_t = \mathbf{M}\mathbf{x}_{t-1} + \eta_t, \quad (1a) $$

$$ \mathbf{y}_t = \mathbf{H}\mathbf{x}_t + \epsilon_t, \quad (1b) $$

where $t$ is the time index and $\eta_t$ and $\epsilon_t$ are unbiased Gaussian vectors, representing the model and observation errors, respectively. Their error covariance matrices are denoted $\mathbf{Q}$ and $\mathbf{R}$, respectively. Those matrices indirectly control the respective weight given to the model and to the observations. It constitutes an important tuning part of the state-space models (see Tandeo et al., 2020, for a more in-depth discussion).

In such a data-driven problem where only a part of the system is observed, a first natural step is to consider that the state $\mathbf{x}$ is directly related to the observations $\mathbf{y}$. For instance, in the example of the Lorenz-63 system introduced previously, observations correspond to the second and third components of the system (i.e., $\mathbf{x}_2$ and $\mathbf{x}_3$, formally defined later).

In this paper, we propose introducing a hidden vector denoted $\mathbf{z}$, corresponding to one or more hidden components that are not observed. For this purpose, the state is augmented using this hidden component $\mathbf{z}$, the observation vector $\mathbf{y}$ does not change, and the operator $\mathbf{H}$ is a truncated identity matrix. The use of augmented state space is classic in data assimilation and mostly refers to the estimation of unknown parameters of the dynamical model (see Ruiz et al., 2013, for further details).

The hidden vector $\mathbf{z}$ is now accounted in the linear model $\mathbf{M}$ given in Eq. (1a) whose dimension has increased. The hidden components are completely unknown and thus randomly initialized using Gaussian white noises and are parameterized by $\sigma^2$, their level of variance. The next step is to infer $\mathbf{z}$

using a statistical estimation method. Starting from the random initialization, an iterative procedure is proposed based on the maximization of the likelihood.

The proposed approach is based on a linear and Gaussian state-space model given in Eq. (1) and thus uses the classic Kalman filter and smoother equations. The Kalman filter (forward in time) is used to get the information of the likelihood, whereas the Kalman smoother (forward and backward in time) is used to get the best estimate of the state. The proposed approach is inspired by the expectation-maximization algorithm (denoted EM; see Shumway and Stoffer, 1982) and is able to iteratively estimate the matrices $\mathbf{M}$ and $\mathbf{Q}$. In this paper, $\mathbf{R}$ is assumed to be known and negligible. The criterion used to update those matrices is based on the innovations defined by the difference between the observations $\mathbf{y}$ and the forecast of the model $\mathbf{M}$, denoted $\mathbf{x}^f$. The likelihood of the innovations, denoted $\mathcal{L}$, is computed using $T$ time steps such that

$$ \mathcal{L} \triangleq p(\mathbf{y}_1, \dots, \mathbf{y}_T | \mathbf{x}_1^f, \dots, \mathbf{x}_T^f) \propto \prod_{t=1}^{T} \exp \left( - (\mathbf{y}_t - \mathbf{H}\mathbf{x}_t^f)^{\top} \mathbf{\Sigma}_t^{-1} (\mathbf{y}_t - \mathbf{H}\mathbf{x}_t^f) \right), \quad (2) $$

where $\mathbf{\Sigma}_t = \mathbf{H}\mathbf{P}_t^f\mathbf{H}^\top + \mathbf{R}$, with $\mathbf{P}_t^f = \mathbf{M}\mathbf{P}_{t-1}^a\mathbf{M}^\top + \mathbf{Q}$ and $\mathbf{P}_{t-1}^a$ corresponding to the state covariance estimated by the Kalman filter at time $t-1$. The innovation likelihood given in Eq. (2) is interesting because it corresponds to the squared distance between the observations and the forecast normalized by their uncertainties, represented by the covariance $\mathbf{\Sigma}_t$.

At each iteration of the augmented Kalman procedure, the estimate of the matrix $\mathbf{M}$ is given by the least-square estima-
tor, using a linear regression such that

$$M^{(i)} = \sum_{t=2}^{T} \frac{\left( \mathbf{x}_{t-1}^{(i-1)} (\mathbf{x}_{t-1}^{(i-1)})^{\top} \right)^{-1} \mathbf{x}_{t}^{(i-1)} (\mathbf{x}_{t-1}^{(i-1)})^{\top}}{T-1}, \quad (3)$$

where $\mathbf{x}^{(i-1)}$ corresponds to the output catalog of the previous iteration (the result of a Kalman smoothing and a Gaussian sampling, explained in more detail below). Following Eq. (1a), the covariance $Q$ is estimated empirically using the estimate of $M$ given in Eq. (3), such that

$$Q^{(i)} = \sum_{t=2}^{T} \frac{\left( \mathbf{x}_{t}^{(i-1)} - M^{(i)} \mathbf{x}_{t-1}^{(i-1)} \right) \left( \mathbf{x}_{t}^{(i-1)} - M^{(i)} \mathbf{x}_{t-1}^{(i-1)} \right)^{\top}}{T-1}. \quad (4)$$

Then, a Kalman smoother is applied using the $M^{(i)}$ and $Q^{(i)}$ matrices estimated in Eqs. (3) and (4). At each time $t$, it results in a Gaussian mean vector $\mathbf{x}_t^s$ and a covariance matrix $\mathbf{P}_t^s$. As input of the next iteration of the algorithm, the catalog $\mathbf{x}^{(i)}$ is updated using a Gaussian random sampling using $\mathbf{x}_t^s$ and $\mathbf{P}_t^s$ at each time $t$. This random sampling is used to exploit the linear correlations between the components of the state vector that appear in the nondiagonal terms of $\mathbf{P}^s$. The random sampling is also used to avoid being trapped in a local maximum, as in stochastic EM procedures (Delyon et al., 1999).

The likelihood calculated at each iteration of the procedure increases until convergence. The algorithm is stopped when the likelihood difference between two iterations becomes small. The solutions of the proposed method are the last Gaussian mean vectors $\mathbf{x}_t^s$ and covariance matrices $\mathbf{P}_t^s$ calculated at each time $t$. The component corresponding to the latent component $z$ is finally retrieved with information on its uncertainty.

### 3 Experiment and evaluation metrics

The methodology is tested on the Lorenz-63 system (Lorenz, 1963). This three-dimensional dynamical system models the evolution of the convection ($x_1$) as a function of horizontal ($x_2$) and vertical temperature gradients ($x_3$). The evolution of the system is governed by three ordinary differential equations, i.e.,

$$\dot{x}_1 = 10(x_2 - x_1), \quad (5a)$$

$$\dot{x}_2 = x_1(28 - x_3) - x_2, \quad (5b)$$

$$\dot{x}_3 = x_1 x_2 - \frac{8}{3} x_3. \quad (5c)$$

Runge-Kutta 4-5 is used to integrate the Lorenz-63 equations to generate $x_1, x_2$ and $x_3$. In this paper, it is assumed that $x_1$ is never observed: only $x_2$ and $x_3$ are observed on 10 model time units of the Lorenz-63 system every $dt = 0.001$ time steps (Fig. 2a). The observation vector is thus

$y = [y_2, y_3]$. In what follows, only those data are available, not the set of Eq. (5).

The methodology is applied to the Lorenz-63 system, adding sequentially a new hidden component in the state of the system as follows. At the beginning, the state is augmented such that $\mathbf{x} = [x_2, x_3, z_1]$, where $z_1$ is randomly initialized with a white noise, with variance $\sigma^2 = 5$. The observations are stored in the vector $\mathbf{y} = [y_2, y_3]$. The observation operator is thus the $2 \times 3$ matrix $\mathbf{H} = [1, 0, 0|0, 1, 0]$. After 30 iterations of the algorithm presented in Sect. 2, the hidden component $z_1$ has converged. After that, a new white noise $z_2$ is used to augment the state such that $\mathbf{x} = [x_2, x_3, z_1, z_2]$, the vector $\mathbf{y} = [y_2, y_3]$ remains the same, and the iterative algorithm is applied until stabilization of $z_2$. As long as the stabilized likelihood continues to increase with the addition of a hidden component, this state-augmentation procedure is repeated.

Note that several hidden components can be added all at once, with a similar performance to the sequential procedure described above (results not shown). In this all-at-once case, the interpretation of the retrieved components is not as informative, and thus we decided to retain the sequential case. Note also that the methodology has been tested with larger $dt$ (i.e., 0.01 and 0.1). The conclusion is that, by increasing the time delay between observations, it significantly increases the number of latent variables (results not shown). Finally, the assimilation window length corresponds to $10^4$ time steps. By reducing this length (e.g., to $10^3$, $10^2$ or $10^1$), the conclusions remain the same as for $dt = 0.001$.

### 4 Results

Using the experiment presented in Sect. 3, three hidden components $z_1, z_2$ and $z_3$ were sequentially added. They are reported in Fig. 2 with the true Lorenz components $x_1, x_2$ and $x_3$. Although they do not fit the hidden variable $x_1$ of the Lorenz system, the first two hidden components $z_1$ and $z_2$ show time variations. By contrast, $z_3$ remains close to 0, with a large confidence interval. This suggests that our method has identified that two hidden variables are enough to retrieve the dynamics of the two observed variables. This result is consistent with the effective dimension of the Lorenz-63 system, which is between two and three. Here, as the estimated dynamical model $M$ is a linear approximation, the dimension of the augmented state and the observed components is higher than the effective one.

This is confirmed by the evaluation of the likelihood of the observations $y_2$ and $y_3$ with different linear models, obtained with or without the use of hidden components $z$ (Fig. 3). This likelihood is useful for diagnosing the optimal number of dimensions needed to emulate the dynamics of the observed components. As the proposed method is stochastic, 50 independent realizations of the likelihood are shown for each experiment. The 50 realizations vary from the random values
<figure><img src="image_43.png" /><figcaption>Figure 2. True components of the Lorenz-63 model (a) and hidden components estimated using the iterative and augmented Kalman procedure (b). The shaded colors correspond to the 95 % Gaussian confidence intervals.</figcaption></figure>

given to the added hidden variable at the beginning of the iterative procedure. In the naive case where the state of the system is $[x_2, x_3]$ (black dashed line), the likelihood is small. Then, adding successively $z_1$ (green lines) and $z_2$ (red lines), after 30 iterations of the proposed algorithm, the likelihood significantly increases. Finally, due to a significant increase in the forecast covariance $\mathbf{P}^f$ in Eq. (2), the inclusion of $z_3$ reduces the likelihood (purple lines). This suggests that a third variable is not needed and is even detrimental to the skill of the reconstruction. Those results indicate that the best linear model for predicting the variations of the observations $y_2$ and $y_3$ is the one using two hidden components. Thus, for the rest of the paper, the focus is placed on the model with the following augmented state: $\mathbf{x} = [x_2, x_3, z_1, z_2]$.

The question is now the following: what is the significance of those hidden components $z_1$ and $z_2$ estimated using the proposed methodology? Are they correlated with the unobserved component $x_1$ or with the observed ones $x_2$ and $x_3$? Are they somehow proxies of the unobserved component? Using symbolic regression (i.e., using basic mathematical transformations of $x_2$ and $x_3$ as regressors to explain $z_1$ and $z_2$), it has been found that the hidden components $z$ correspond to linear combinations of the derivatives of the observations such that

$$z_1 = a_2 \dot{x}_2 + a_3 \dot{x}_3, \qquad (6a)$$

$$z_2 = b_1 \dot{z}_1 + b_2 \dot{x}_2 + b_3 \dot{x}_3. \qquad (6b)$$

When developing Eq. (6b) using Eq. (6a), the second hidden component is written as $z_2 = b_2 \dot{x}_2 + b_3 \dot{x}_3 + b_1 a_2 \dot{x}_2 + b_1 a_3 \dot{x}_3$. It shows that $z_1$ uses the first derivative of $x_2$ and $x_3$, whereas $z_2$ uses the second derivatives. This result makes the link with the Taylor and Takens theorem, which shows that an unobserved component (i.e., $x_1$) can be replaced by the observed components (i.e., $x_2$ and $x_3$) at different time lags.

Note that, due to the stochastic behavior of the algorithm, the $a$ and $b$ coefficients are not fixed, and several combinations of them can reach the same performance in terms of likelihood. This is illustrated in Fig. 3a, with 50 independent realizations of the proposed algorithm. When considering only $z_1$ (green lines), the algorithm converges to various solutions but is mainly restricted around two solutions (corresponding to a minimum and a maximum of likelihood). As shown in Fig. 3b, the minimum likelihood corresponds to $a_3 = 0$ and the maximum likelihood corresponds to $a_2 = 0$. Thus, the likelihood of $z_1 = a_3 \dot{x}_3$ is higher than $z_1 = a_2 \dot{x}_2$. This suggests that $\dot{x}_3$ is more important than $\dot{x}_2$ in explaining the variations of the Lorenz system (this is consistent with the investigation of Sévellec and Fedorov, 2014, in a modified version of the Lorenz-63 model). Interestingly, the scatter plot between $a_2$ and $a_3$ shows a circular relationship. This is also the case for $b_2$ and $b_3$ (results not shown). Then, in Fig. 3a, when considering $z_1$ and $z_2$ (red lines), the 50 independent realizations reach the same likelihood after 30 iterations. This means that if $a_3 = 0$ when considering only $z_1$, then $b_3 \neq 0$ when introducing $z_2$. In terms of forecast performance, this is similar to $a_2 = 0$ and $b_2 \neq 0$, because the likelihoods converge to the same value (red lines after 30 iterations).

To compare the performance of the naive linear model $\mathbf{M}$ with $[x_2, x_3]$ and the ones with $[x_2, x_3, z_1]$ or $[x_2, x_3, z_1, z_2]$, their forecasts are evaluated. After applying the proposed algorithm, the $\hat{\mathbf{M}}$ and $\hat{\mathbf{Q}}$ estimated matrices are used to derive probabilistic forecast, starting from the last available observation $\mathbf{y}_t$, using

$$E[\mathbf{x}_{t+1} | \mathbf{y}_1, \dots, \mathbf{y}_t] = \hat{\mathbf{M}} E[\mathbf{x}_t | \mathbf{y}_1, \dots, \mathbf{y}_t], \qquad (7a)$$

$$\text{Cov}[\mathbf{x}_{t+1} | \mathbf{y}_1, \dots, \mathbf{y}_t] = \hat{\mathbf{M}} \text{Cov}[\mathbf{x}_t | \mathbf{y}_1, \dots, \mathbf{y}_t] \hat{\mathbf{M}}^T + \hat{\mathbf{Q}}, \qquad (7b)$$
<figure><img src="image_44.png" /><figcaption>**Figure 3.** Likelihoods as a function of the iteration of the augmented Kalman procedure (a) and estimation of the $a_2$ and $a_3$ parameters (b). Different dynamical models are considered, from none to three hidden components in $z$, whereas only $x_2$ and $x_3$ are observed in the Lorenz-63 model. The likelihoods of 50 independent realizations of the iterative and augmented Kalman procedure are shown.</figcaption></figure>

with $E$ and Cov the expectation and the covariance, respectively. To test the predictability of the different linear models (i.e., with or without hidden components $z$), a test set has been created, starting from the end of the sequence of observations ($y_1, \dots, y_T$) used in the assimilation window. This test set also corresponds to $10^4$ time steps with $dt = 0.001$. It is used to compute two metrics, the root mean square error (RMSE) and the coverage probability at 50%. The RMSE is used to evaluate the precision of the forecasts, comparing the true $x_2$ and $x_3$ components to the estimated ones, whereas the coverage probability is used to evaluate the reliability of the prediction, evaluating the proportion of true trajectories falling within the 50% prediction interval of $x_2$ and $x_3$. Examples of predictions are given in Fig. 4. It shows bad linear predictions of the model with only [$x_2, x_3$] (dashed black lines). As the $M$ operator is not time-dependent, the predictions are quite similar, close to the persistence. Then, adding one (green) or two (red) hidden components in the $M$ operators creates some nonlinearities in the forecasts.

In Fig. 5, the predictions are evaluated over the whole test dataset for different lead times. By introducing hidden components, the RMSE decreases for both $x_2$ and $x_3$ components (panels a and b). For instance, for a lead time of 0.05, when considering two hidden components, the RMSE is halved when it is compared to the naive linear model without hidden components. The coverage probability metric is also largely improved (panels c and d). Indeed, the results with two hidden components are close to 50%, the optimal value.

To evaluate where the linear model with [$x_2, x_3, z_1, z_2$] performs better than the one with [$x_2, x_3$], the Euclidean distances between the forecasts (for a lead time of 0.1) and the truth are computed. Those errors are evaluated at each time step of the test dataset, in the ($x_2, x_3$) space. Based on those errors, Fig. 6 shows the relative improvement between the model without and the model with hidden components. When the two models have similar performance, values are close to 0 (white), and when the model including $z_1$ and $z_2$ is better, values are close to 1 (red). Figure 6 clearly shows that error reduction is not homogeneous in the attractor. The

improvement is moderate on the outside of the wings of the attractor but important in the wing transition. This suggests that the introduction of the hidden components $z_1$ and $z_2$ makes it possible to provide information on the position in the attractor and thus to make better predictions, especially in bifurcation regions.

# 5 Conclusions

In this article, the goal is to retrieve hidden components of a dynamical system that is partially observed. The proposed methodology is purely data-driven, not physics-driven (i.e., without the use of any equations of the dynamical model). It is based on the combination of data assimilation and machine learning techniques. Three main ideas are used in the methodology: an augmented state strategy, a linear approximation of a dynamical system and an iterative procedure. The methodology is easy to implement using simple strategies and well-established algorithms: Kalman filter and smoother, linear regression using least squares, an iterative procedure inspired by the EM recursions and Gaussian random sampling for the stochastic aspect.

The methodology is tested on the Lorenz-63 system, where only two components of the system are observed in a short period of time. Several hidden components are introduced sequentially in the system. Although the hidden components are initialized randomly, only a few iterations of the proposed algorithm are necessary to retrieve relevant information. The recovered components are expressed with Gaussian distributions. The new components correspond to linear combinations of successive derivatives of the observed variables. This result is consistent with the theorems of Taylor and Takens, which show that time-delay embedding is useful for improving the forecasts of the system. In our case, this is evaluated using the likelihood, a metric that evaluates the innovation (i.e., the difference between Gaussian forecasts and Gaussian observations).

Using our methodology, we do not retrieve the true missing Lorenz component and need two hidden variables to rep-
<figure><img src="image_45.png" /><figcaption>**Figure 4.** Example of three statistical forecasts of $x_2$ (a) and $x_3$ (b) with their 50% prediction interval using three different linear operators with no hidden component (dashed black), one hidden component (green) and two hidden components (red). These predictions are obtained using sequential statistical forecasts, as explained in Eq. (7), on an independent test dataset.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>**Figure 5.** Root mean square error (a, b) and 50% coverage probability (c, d) as a function of the lead time (x axis) for the reconstruction of the components $x_2$ (a, c) and $x_3$ (b, d). These metrics are evaluated on an independent test dataset.</figcaption></figure>

present a single missing one. The reason for this mismatch is two-fold and is mainly the linear approximation of the dynamical system, which implies that (1) the true missing component, which does not have to be linear combinations of the observed variables, is impossible to retrieve in our framework and (2) two variables, using combinations of the time derivatives of the observed variables, are needed to accurately represent the complexity of the dynamics. However, it is important to note that, even if two variables are needed to replace a single one, the dynamical evolution of the system is relatively well captured, for short lead times, with our methodology. This correct representation of the evolution

might ultimately be the most important (e.g., for accurate and reliable forecasting).

The proposed methodology uses a strong assumption: the linear approximation of the dynamical system is global (i.e., fixed for the whole observation period). A perspective is to use adaptive approximations of the model using local linear regressions. This strategy is computationally more expensive because a linear regression is adjusted at each time step but shows some improvements in chaotic systems (see Platzer et al., 2021a, b). In this context of an adaptive linear dynamical model, the proposed methodology could be easily plugged into an ensemble Kalman procedure based on
<figure><img src="image_47.png" /><figcaption>**Figure 6.** Relative forecast improvement measured as 1 minus the ratio between two Euclidean distances: the one calculated with model [$x_2, x_3, z_1, z_2$] (at the numerator) and the one calculated with model [$x_2, x_3$] (at the denominator). The Euclidean distances are calculated in the ($x_2, x_3$) space and correspond to the error between the forecasts (for a lead time of 0.1) and the truth, evaluated on an independent test dataset.</figcaption></figure>

analog forecasts (Lguensat et al., 2017). In future works, we plan to compare the global and local linear approaches (i.e., a fixed or adaptive linear surrogate model). We also plan to compare them to nonlinear surrogate models, based on neural network architectures with latent information encoded in an augmented space or in hidden layers (e.g., long short-term memory – LSTM).

In this paper, we have demonstrated the feasibility of the method on an idealized and comprehensive problem using the Lorenz-63 system. In the future, we plan to apply the methodology to more challenging problems, like the Lorenz-96 system or a quasi-geostrophic model. For application to real data, we plan to use a database of observed climate indices and try to find latent variables that help to make data-driven predictions.

**Code and data availability.** The Python code is available at https://github.com/ptandeo/Kalman under the GNU license and the data are generated using the Lorenz-63 system (Lorenz, 1963).

**Supplement.** The supplement related to this article is available online at: https://doi.org/10.5194/npg-30-129-2023-supplement.

**Author contributions.** PT wrote the article. PT and PA developed the algorithm. FS and PA helped with the redaction of the paper.

**Competing interests.** At least one of the (co-)authors is a member of the editorial board of *Nonlinear Processes in Geophysics*.

The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Acknowledgements.** This paper is the result of a project proposed in a course on "Data Assimilation" in the masters program "Ocean Data Science" at Univ Brest, ENSTA Bretagne, and IMT Atlantique, France. The authors would like to thank the students for their participation in the project: Nils Niebaum, Zackary Vanche, Benoit Presse, Dimitri Vlahopoulos, Yanis Grit and Joséphine Schmutz. The authors would like to thank Noémie Le Carrer for her proofreading of the paper and Paul Platzer, Said Ouala, Lucas Drumetz, Juan Ruiz, Manuel Pulido and Takemasa Miyoshi for their valuable comments.

**Financial support.** This work was supported by ISblue project, Interdisciplinary graduate school for the blue planet (ANR-17-EURE-0015) and co-funded by a grant from the French government under the program "Investissements d'Avenir" embedded in France 2030. This work was also supported by LEFE program (LEFE IMAGO projects ARVOR).

**Review statement.** This paper was edited by Natale Alberto Car-rassi and reviewed by two anonymous referees.
## References

Bocquet, M., Brajard, J., Carrassi, A., and Bertino, L.: Data assimilation as a learning tool to infer ordinary differential equation representations of dynamical models, Nonlin. Processes Geophys., 26, 143–162, https://doi.org/10.5194/npg-26-143-2019, 2019.

Brajard, J., Carrassi, A., Bocquet, M., and Bertino, L.: Combining data assimilation and machine learning to emulate a dynamical model from sparse and noisy observations: A case study with the Lorenz 96 model, Journal of Computational Science, 44, 101171, https://doi.org/10.1016/j.jocs.2020.101171, 2020.

Brajard, J., Carrassi, A., Bocquet, M., and Bertino, L.: Combining data assimilation and machine learning to infer unresolved scale parametrization, Philos. T. Roy. Soc. A, 379, 2194, https://doi.org/10.1098/rsta.2020.0086, 2021.

Brunton, S. L., Proctor, J. L., and Kutz, J. N.: Discovering governing equations from data by sparse identification of nonlinear dynamical systems, P. Natl. Acad. Sci. USA, 113, 3932–3937, 2016.

Brunton, S. L., Brunton, B. W., Proctor, J. L., Kaiser, E., and Kutz, J. N.: Chaos as an intermittently forced linear system, Nat. Commun., 8, 19, https://doi.org/10.1038/s41467-017-00030-8, 2017.

Delyon, B., Lavielle, M., and Moulines, E.: Convergence of a stochastic approximation version of the EM algorithm, Ann. Stat., 27, 94–128, 1999.

Dreano, D., Tandeo, P., Pulido, M., Ait-El-Fquih, B., Chonavel, T., and Hoteit, I.: Estimating model-error covariances in nonlinear state-space models using Kalman smoothing and the expectation–maximization algorithm, Q. J. Roy. Meteor. Soc., 143, 1877–1885, 2017.

Fablet, R., Chapron, B., Drumetz, L., Mémin, E., Pannekoucke, O., and Rousseau, F.: Learning variational data assimilation models and solvers, J. Adv. Model. Earth Sy., 13, e2021MS002572, https://doi.org/10.1029/2021MS002572, 2021.

Jayne, S. R., Roemmich, D., Zilberman, N., Riser, S. C., Johnson, K. S., Johnson, G. C., and Piotrowicz, S. R.: The Argo program: present and future, Oceanography, 30, 18–28, 2017.

Kitagawa, G.: A self-organizing state-space model, J. Am. Stat. Assoc., 93, 1203–1215, 1998.

Korda, M. and Mezić, I.: Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control, Automatica, 93, 149–160, 2018.

Lguensat, R., Tandeo, P., Ailliot, P., Pulido, M., and Fablet, R.: The analog data assimilation, Mon. Weather Rev., 145, 4093–4107, 2017.

Lorenz, E. N.: 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 (data available at: https://github.com/ptandeo/Kalman, last access: 26 May 2023).

Mangiarotti, S. and Huc, M.: Can the original equations of a dynamical system be retrieved from observational time series?, Chaos: An Interdisciplinary Journal of Nonlinear Science, 29, 023133, https://doi.org/10.1063/1.5081448, 2019.

North, J. S., Wikle, C. K., and Schliep, E. M.: A Bayesian Approach for Data-Driven Dynamic Equation Discovery, J. Agr. Biol. Envir. S., 27, 728–747, https://doi.org/10.1007/s13253-022-00514-1, 2022.

Ouala, S., Nguyen, D., Drumetz, L., Chapron, B., Pascual, A., Collard, F., Gaultier, L., and Fablet, R.: Learning latent dynamics for partially observed chaotic systems, Chaos: An Interdisciplinary Journal of Nonlinear Science, 30, 103121, https://doi.org/10.1063/5.0019309, 2020.

Platzer, P., Yiou, P., Naveau, P., Filipot, J.-F., Thiébaut, M., and Tandeo, P.: Probability distributions for analog-to-target distances, J. Atmos. Sci., 78, 3317–3335, 2021a.

Platzer, P., Yiou, P., Naveau, P., Tandeo, P., Filipot, J.-F., Ailliot, P., and Zhen, Y.: Using local dynamics to explain analog forecasting of chaotic systems, J. Atmos. Sci., 78, 2117–2133, 2021b.

Pulido, M., Tandeo, P., Bocquet, M., Carrassi, A., and Lucini, M.: Stochastic parameterization identification using ensemble Kalman filtering combined with maximum likelihood methods, Tellus A, 70, 1442099, https://doi.org/10.1080/16000870.2018.1442099, 2018.

Rudy, S. H., Brunton, S. L., Proctor, J. L., and Kutz, J. N.: Data-driven discovery of partial differential equations, Science Advances, 3, e1602614, https://doi.org/10.1126/sciadv.1602614, 2017.

Ruiz, J. J., Pulido, M., and Miyoshi, T.: Estimating model parameters with ensemble-based data assimilation: A review, J. Meteorol. Soc. Jpn., 91, 79–99, 2013.

Sévellec, F. and Fedorov, A. V.: Millennial variability in an idealized ocean model: predicting the AMOC regime shifts, J. Climate, 27, 3551–3564, 2014.

Shumway, R. H. and Stoffer, D. S.: An approach to time series smoothing and forecasting using the EM algorithm, J. Time Ser. Anal., 3, 253–264, 1982.

Takens, F.: Detecting strange attractors in turbulence, in: Dynamical systems and turbulence, Warwick 1980, Springer, 366–381, ISBN 978-3-540-38945-3, https://doi.org/10.1007/BFb0091924, 1981.

Talmon, R., Mallat, S., Zaveri, H., and Coifman, R. R.: Manifold learning for latent variable inference in dynamical systems, IEEE T. Signal Proces., 63, 3843–3856, 2015.

Tandeo, P., Ailliot, P., Ruiz, J., Hannart, A., Chapron, B., Cuzol, A., Monbet, V., Easton, R., and Fablet, R.: Combining analog method and ensemble data assimilation: application to the Lorenz-63 chaotic system, in: Machine learning and data mining approaches to climate science, Springer, 3–12, ISBN 978-3-319-17220-0, https://doi.org/10.1007/978-3-319-17220-0_1, 2015.

Tandeo, P., Ailliot, P., Bocquet, M., Carrassi, A., Miyoshi, T., Pulido, M., and Zhen, Y.: A review of innovation-based methods to jointly estimate model and observation error covariance matrices in ensemble data assimilation, Mon. Weather Rev., 148, 3973–3994, 2020.

Wikner, A., Pathak, J., Hunt, B. R., Szunyogh, I., Girvan, M., and Ott, E.: Using data assimilation to train a hybrid forecast system that combines machine-learning and knowledge-based components, Chaos: An Interdisciplinary Journal of Nonlinear Science, 31, 053114, https://doi.org/10.1063/5.0048050, 2021.
## 4.2 Tandeo, Ailliot, Bocquet, Carrassi, Miyoshi, Pulido et Zhen (2020) [MWR]

**Contexte** D'un point de vue mathématique, l'assimilation de données se base sur une modélisation espace-état, avec une équation de la dynamique et une équation d'observation. Ces deux équations ne sont pas déterministes et prennent en compte d'éventuelles erreurs du modèle et des observations. Depuis mes travaux de thèse, j'ai pris conscience de la nécessité d'une bonne quantification de ces incertitudes. Cet aspect, primordial pour toutes les méthodes d'assimilation de données, était jusque là relativement peu référencé dans la bibliographie. C'est pour cela qu'avec des collègues travaillant sur des aspects méthodologiques, nous avons décidé en 2017, suite à un workshop, de rédiger un papier d'état de l'art sur la quantification jointe des incertitudes du modèle (l'ébauche) et des observations. Au même moment, MWR incitait les chercheurs à publier ce type d'articles : ceci motiva le choix de ce journal. Après plusieurs mois de rédaction et de nombreuses révisions, l'article parut finalement en 2020.

**Résumé** L'assimilation de données combine les prévisions d'un modèle numérique avec des observations. La plupart des algorithmes d'assimilation considèrent les termes d'erreur de modèle et d'observation comme étant des bruits additifs Gaussiens, spécifiés par leurs matrices de covariance Q et R. Ces covariances d'erreurs, plus particulièrement leurs amplitudes, déterminent les pondérations accordées aux termes d'ébauche et d'observation. Par conséquent, les matrices Q et R influent considérablement sur le résultat de l'assimilation de données. Cet article de l'état de l'art vise à présenter et à discuter, dans un cadre mathématique unifié, différentes méthodes pour estimer conjointement les matrices Q et R à l'aide de méthodes d'assimilation de données ensemblistes. La plupart des méthodes développées à ce jour utilisent l'innovation, définie comme la différence entre les observations et les prévisions, projetées dans l'espace des observations. Ces méthodes reposent sur deux principaux critères statistiques : 1) la méthode des moments, dans laquelle les moments théoriques et empiriques des innovations sont supposées égaux et 2) les méthodes qui utilisent la vraisemblance des observations, contenue dans les innovations. Les méthodes examinées supposent que les innovations sont des variables aléatoires Gaussiennes, bien que l'extension à d'autres distributions soit possible pour les méthodes basées sur la vraisemblance. Les méthodes montrent aussi certaines différences en termes de niveaux de complexité et d'applicabilité aux systèmes de grande dimension. En conclusion, nous discutons des principaux défis pour améliorer les méthodes d'estimation de Q et R. Ces défis incluent la prise en compte des covariances d'erreur variant dans le temps, l'utilisation d'un nombre limité d'observation, l'estimation de termes d'erreur déterministes supplémentaires ou encore la prise en compte de bruits corrélés.
# REVIEW

## A Review of Innovation-Based Methods to Jointly Estimate Model and Observation Error Covariance Matrices in Ensemble Data Assimilation

PIERRE TANDEO,ª,ᵇ PIERRE AILLIOT,ᶜ MARC BOCQUET,ᵈ ALBERTO CARRASSI,ᵉ,,f,g

TAKEMASA MIYOSHI,ᵇ MANUEL PULIDO,ᵇ,i,e AND YICUN ZHENª

(Manuscript received 18 July 2019, in final form 12 June 2020)

### ABSTRACT

Data assimilation combines forecasts from a numerical model with observations. Most of the current data assimilation algorithms consider the model and observation error terms as additive Gaussian noise, specified by their covariance matrices **Q** and **R**, respectively. These error covariances, and specifically their respective amplitudes, determine the weights given to the background (i.e., the model forecasts) and to the observations in the solution of data assimilation algorithms (i.e., the analysis). Consequently, **Q** and **R** matrices significantly impact the accuracy of the analysis. This review aims to present and to discuss, with a unified framework, different methods to jointly estimate the **Q** and **R** matrices using ensemble-based data assimilation techniques. Most of the methods developed to date use the innovations, defined as differences between the observations and the projection of the forecasts onto the observation space. These methods are based on two main statistical criteria: 1) the method of moments, in which the theoretical and empirical moments of the innovations are assumed to be equal, and 2) methods that use the likelihood of the observations, themselves contained in the innovations. The reviewed methods assume that innovations are Gaussian random variables, although extension to other distributions is possible for likelihood-based methods. The methods also show some differences in terms of levels of complexity and applicability to high-dimensional systems. The conclusion of the review discusses the key challenges to further develop estimation methods for **Q** and **R**. These challenges include taking into account time-varying error covariances, using limited observational coverage, estimating additional deterministic error terms, or accounting for correlated noise.

KEYWORD: Data assimilation

## 1. Introduction

In meteorology and other environmental sciences, an important challenge is to estimate the state of the system as accurately as possible. In meteorology, this state includes pressure, humidity, temperature and wind at different locations and elevations in the atmosphere. Data assimilation (DA) refers to mathematical methods

that use both model predictions (also called background information) and partial observations to retrieve the current state vector with its associated error. An accurate estimate of the current state is crucial to get good forecasts, and it is particularly so whenever the system dynamics is chaotic, such as it is the case for the atmosphere.

The performance of a DA system to estimate the state depends on the accuracy of the model predictions, the observations, and their associated error terms. A simple, popular and mathematically justifiable way of modeling these errors is to assume them to be independent and

<footnote>ªIMT Atlantique, Lab-STICC, UMR CNRS 6285, Brest, France</footnote>

<footnote>ᵇRIKEN Center for Computational Science, Kobe, Japan</footnote>

<footnote>ᶜLMBA, UMR CNRS 6205, University of Brest, Brest, France</footnote>

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

<footnote>ᵉDepartment of Meteorology, University of Reading, Reading, United Kingdom</footnote>

<footnote>fNational Centre for Earth Observation, University of Reading, Reading, United Kingdom</footnote>

<footnote>gMathematical Institute, University of Utrecht, Utrecht, Netherlands</footnote>

<footnote>hUniversidad Nacional del Nordeste, Corrientes, Argentina</footnote>

<footnote>iCONICET, Corrientes, Argentina</footnote>

<footnote>Corresponding author: Pierre Tandeo, pierre.tandeo@imt-atlantique.fr</footnote>
<figure><img src="image_48.png" /><figcaption>FIG. 1. Sketch of sequential and ensemble DA algorithms in the observation space (i.e., in the space of the observations $\mathbf{y}$), where the observation operator $\mathcal{H}$ is omitted for simplicity. The ellipses represent the forecast $\mathbf{P}^f$ and analysis $\mathbf{P}^a$ error covariances, whereas the model $\mathbf{Q}$ and observation $\mathbf{R}$ error covariances are the unknown entries of the state-space model in Eqs. (1) and (2). The forecast error covariance matrix is written $\mathbf{P}^f$ and is the sum of $\mathbf{P}^m$, the forecast state $\mathbf{x}^f$ spread, and the model error $\mathbf{Q}$. This scheme is a modified version that is based on Fig. 1 from Carrassi et al. (2018).</figcaption></figure>

unbiased Gaussian white noise, with covariance matrices $\mathbf{Q}$ for the model and $\mathbf{R}$ for the observations. Given the aforementioned importance of $\mathbf{Q}$ and $\mathbf{R}$ in estimating the analysis state and error, a number of studies dealing with this problem has arisen in the last decades. This review presents and summarizes the different techniques used to estimate simultaneously the $\mathbf{Q}$ and $\mathbf{R}$ covariances. Before discussing the methods to achieve this goal, the mathematical formulation of DA is briefly introduced.

### a. Problem statement

Hereinafter, the unified DA notation proposed in Ide et al. (1997) is used.<fnref n="1" /> Data assimilation algorithms are used to estimate the state of a system, $\mathbf{x}$, conditionally on observations, $\mathbf{y}$. A classic strategy is to use sequential and ensemble DA frameworks, as illustrated in Fig. 1, and to combine two sources of information: model forecasts (in green) and observations (in blue). The ensemble framework uses different realizations, also called members, to track the state of the system at each assimilation time step.

The forecasts of the state are based on the usually incomplete and approximate knowledge of the system dynamics. The evolution of the state from time $k-1$ to $k$ is given by the model equation

$$ \mathbf{x}(k) = \mathcal{M}_k[\mathbf{x}(k-1)] + \mathbf{\eta}(k), \quad (1) $$

where the model error $\mathbf{\eta}$ implies that the dynamic model operator $\mathcal{M}_k$ is not perfectly known. Model error is usually assumed to follow a Gaussian distribution with zero mean (i.e., the model is unbiased) and covariance $\mathbf{Q}$. The dynamic model operator $\mathcal{M}_k$ in Eq. (1) has also an explicit dependence on $k$, because it may depend on time-dependent external forcing terms. At time $k$, the forecast state is characterized by the mean of the forecast states, $\mathbf{x}^f$, and its uncertainty matrix, namely $\mathbf{P}^f$, which is also called the background error covariance matrix and is noted as $\mathbf{B}$ in DA.

The forecast covariance $\mathbf{P}^f$ is determined by two processes. The first is the uncertainty propagated from $k-1$ to $k$ by the model $\mathcal{M}_k$ (the green shade within the dashed ellipse in Fig. 1 and denoted by $\mathbf{P}^m$). The second process is the model error covariance $\mathbf{Q}$ accounted by the noise term at time $k$ in Eq. (1). Given that model error is largely unknown and originated by various and diverse sources, the matrix $\mathbf{Q}$ is also poorly known. Model error sources encompass the model $\mathcal{M}$ deficiencies to represent the underlying physics, including deficiencies in the numerical schemes, the cumulative effects of errors in the parameters, and the lack of knowledge of the unresolved scales. Its estimation is a challenge in general, but it is

<footnote n="1">¹Other notations are also used in practice.</footnote>
particularly so in geosciences because we usually have far fewer observations than those needed to estimate the entries of **Q** (Daley 1992; Dee 1995). The sum of the two covariances **P**<sup>m</sup> and **Q** gives the forecast covariance matrix, **P**<sup>f</sup> (full green ellipse in Fig. 1). In the illustration given here, a large contribution of the forecast covariance **P**<sup>f</sup> is due to **Q**. This situation reflects what is common in ensemble DA, where **P**<sup>m</sup> can be too small, as a consequence of the ensemble undersampling of the initial condition error (i.e., the covariance estimated at the previous analysis). In that case, inflating **Q** could partially compensate for the bad specification of **P**<sup>m</sup>.

DA uses a second source of information, the observations **y**, which are assumed to be linked to the true state **x** through the time-dependent operator $\mathcal{H}_k$. This step in DA algorithms is formalized by the observation equation

$$ \mathbf{y}(k) = \mathcal{H}_k[\mathbf{x}(k)] + \epsilon(k), \quad (2) $$

where the observation error $\epsilon$ describes the discrepancy between what is observed and the truth. In practice, it is important to remove as much as possible the large-scale bias in the observation before DA. Then, it is common to state that the remaining error $\epsilon$ follows a Gaussian and unbiased distribution with a covariance **R** (the blue ellipse in Fig. 1). This covariance takes into account errors in the observation operator $\mathcal{H}$, the instrumental noise and the representation error associated with the observation, typically measuring a higher-resolution state than the model represents. Operationally, a correct estimation of **R** that takes into account all of these effects is often challenging (Janjić et al. 2018).

DA algorithms combine forecasts with observations, based on the model and observation equations, given in Eqs. (1) and (2), respectively. The corresponding system of equations is a nonlinear state-space model. As illustrated in Fig. 1, this Gaussian DA process produces a posterior Gaussian distribution with mean $\mathbf{x}^a$ and covariance $\mathbf{P}^a$ (red ellipse). The system given in Eqs. (1) and (2) is representative of a broad range of DA problems, as described in seminal papers such as Ghil and Malanotte-Rizzoli (1991), and still relevant today as referenced by Houtekamer and Zhang (2016) and Carrassi et al. (2018). The assumptions made in Eqs. (1) and (2) about model and observation errors (additive, Gaussian, unbiased, and mutually independent) are strong, yet convenient from the mathematical and computational point of view. Nevertheless, these assumptions are not always realistic in real DA problems. For instance, in operational applications, systematic biases in the model and in the observations

are recurring problems. Indeed, biases affect significantly the DA estimations and a specific treatment is required; see Dee (2005) for more details.

From Eqs. (1) and (2), noting that $M$, $H$, and $\mathbf{y}$ are given, the only parameters that influence the estimation of $\mathbf{x}$ are the covariance matrices **Q** and **R**. These covariances play an important role in DA algorithms. Their importance was early put forward in Hollingsworth and Lönnberg (1986), in section 4.1 of Ghil and Malanotte-Rizzoli (1991), and in section 4.9 of Daley (1991). The results of DA algorithms highly depend on the two error covariance matrices **Q** and **R**, which have to be specified by the users. But these covariances are not easy to tune. Indeed, their impact is hard to grasp in real DA problems with high-dimensionality and nonlinear dynamics. We thus illustrate the problem with a simple example first.

## b. Illustrative example

In either variational or ensemble-based DA methods, the quality of the reconstructed state (or hidden) vector $\mathbf{x}$ largely depends on the relative amplitudes between the assumed observation and model errors (Desroziers and Ivanov 2001). In Kalman filter-based methods, the signal-to-noise ratio $||\mathbf{P}^f||/||\mathbf{R}||$, where $\mathbf{P}^f$ depends on **Q**, impacts the Kalman gain, which gives the relative weights of the observations against the model forecasts. Here, the $|| ||$ operator represents a matrix norm. For instance, Berry and Sauer (2013) used the Frobenius norm to study the effect of this ratio in the reconstruction of the state in toy models.

The importance of **Q**, **R**, and $||\mathbf{P}^f||/||\mathbf{R}||$ is illustrated with the aid of a toy example, using a scalar state $x$ and simple linear dynamics. This simplified setup avoids several issues typical of realistic DA applications: the large dimension of the state, the strong nonlinearities and the chaotic behavior. In this example, the dynamic model in Eq. (1) is a first-order autoregressive model, denoted by AR(1) and defined by

$$ x(k) = 0.95x(k-1) + \eta(k), \quad (3) $$

with $\eta \sim N(0, Q)$, where the superscript $t$ means “true” and $Q^t = 1$. Furthermore, observations $\mathbf{y}$ of the state are contaminated with an independent additive zero-mean and unit-variance Gaussian noise, such that $R^t = 1$ in Eq. (2) with $\mathcal{H}(x) = x$. The goal is to reconstruct $x$ from the noisy observations $\mathbf{y}$ at each time step. The AR(1) dynamic model defined by Eq. (3) has an autoregressive coefficient close to one, representing a process that evolves slowly over time, and a stochastic noise term $\eta$ with variance $Q^t$. Although the knowledge of these two sources of noise is crucial for the estimation
<figure><img src="image_49.png" /><figcaption>FIG. 2. Example of a univariate AR(1) process generated using Eq. (3) with $Q' = 1$ (red line), noisy observations as in Eq. (2) with $R' = 1$ (black dots), and reconstructions with a Kalman smoother (black lines and gray 95% confidence interval) with different values of $Q$ and $R$, from 0.1 to 10. The optimal values of RMSE and coverage probabilities are 0.71% and 95%, respectively.</figcaption></figure>

problem, identifying them is not an easy task. Given that the dynamic model is linear and the error terms are additive and Gaussian in this simple example, the Kalman smoother provides the best estimation of the state (see section 2 for more details). To evaluate the effect of badly specified $Q$ and $R$ errors on the reconstructed state with the Kalman smoother, different experiments were conducted with values of {0.1, 1, 10} for the ratio $Q/R$ (in this toy example, we use $Q/R$ instead of $||P^t||/||R^t||$ for simplicity).

Figure 2 shows, as a function of time, the true state (red line) and the smoothing Gaussian distributions represented by the 95% confidence intervals (gray shaded) and their means (black lines). We also report the root-mean-square error (RMSE) of the reconstruction and the so-called coverage probability, or percentage of $x$ that falls in the 95% confidence intervals (defined as the mean $\pm$ 1.96 the standard deviation in the Gaussian case). In this synthetic experiment, the best RMSE and coverage probability obtained, applying the Kalman smoother with true $Q' = R' = 1$, are 0.71% and 95%, respectively. Using a small model error variance $Q = 0.1Q'$ in Fig. 2a, the filter gives a large weight to the forecasts given by the quasi-persistent autoregressive dynamic model. On the other hand, with a small

observation error variance $R = 0.1R'$ in Fig. 2b, excessive weight is given to the observation and the reconstructed state is close to the noisy measurements. These results show the negative impact of independently badly scaled $Q$ and $R$ error variances. In the case of overestimated model error variance as in Fig. 2c, the mean reconstructed state vector and thus its RMSE are identical to Fig. 2b. In the same way, overestimated observation error variance like in Fig. 2d gives similar mean reconstruction, as in Fig. 2a. These last two results are due to the fact that in both cases, the ratio $Q/R$ are equal, respectively, to 10 and 0.1. Now, we consider in Figs. 2e and 2f the case where the $Q/R$ ratio is equal to 1, but, respectively, using the simultaneous underestimation and overestimation of model and observation errors. In both cases, the mean reconstructed state is equal to that obtained with the true error variances (i.e., RMSE = 0.71). The main difference is the gray confidence interval, which is supposed to contain 95% of the true trajectory: the spread is clearly underestimated in Fig. 2e and overestimated in Fig. 2f, with coverage probability of 36% and 100%, respectively.

We used a simple synthetic example, but for large dimensional and highly nonlinear dynamics, such an
<figure><img src="image_50.png" /><figcaption>FIG. 3. Timeline of the main methods used in geophysical data assimilation for the joint estimation of **Q** and **R** over the last 15 years. Dee (1995) and Desroziers and Ivanov (2001) are not represented here but are certainly the seminal work of this research field in data assimilation.</figcaption></figure>

underestimation or overestimation of uncertainty may have a strong effect and may cause filters to collapse. The main issue in ensemble-based DA is an underdispersive spread, as in Fig. 2e. In that case, the initial condition spread is too narrow, and model forecasts (starting from these conditions) would be similar and potentially out of the range of the observations. In the case of an overdispersive spread, as in Fig. 2f, the risk is that only a small portion of model forecasts would be accurate enough to produce useful information on the true state of the system. This illustrative example shows how important is the joint tuning of model and observation errors in DA. Since the 1990s, a substantial number of studies have dealt with this topic.

### c. Seminal work in the data assimilation community

In a seminal paper, Dee (1995) proposed an estimation method for parametric versions of **Q** and **R** matrices. The method, based on maximizing the likelihood of the observations, yields an estimator that is a function of the innovation defined by $y - H(x^f)$. Maximization is performed at each assimilation step, with the current innovation computed from the available observations. This technique was later extended to estimate the mean of the innovation, which depends on the biases in the forecast and in the observations (Dee and da Silva 1999). The method was then applied to realistic cases in Dee et al. (1999), making the maximization of innovation likelihood a promising technique for the estimation of errors in operational forecasts.

Following a distinct path, Desroziers and Ivanov (2001) proposed using the observation-minus-analysis diagnostic. It is defined by $y - H(x^a)$ with $x^a$ being the analysis (i.e., the output of DA algorithms). The authors proposed

an iterative optimization technique to estimate a scaling factor for the background **B** = **P**^f and observation **R** matrices. The procedure was shown to converge to a proper fixed point. As in Dee's work, the fixed-point method presented in Desroziers and Ivanov (2001) is applied at each assimilation step, with the available observations at the current step.

Later, Chapnik et al. (2004) showed that the maximization of the innovation likelihood proposed by Dee (1995) makes the observation-minus-analysis diagnostic of Desroziers and Ivanov (2001) optimal. Moreover, the techniques of Dee (1995) and Desroziers and Ivanov (2001) have been further connected to the generalized cross-validation method previously developed by statisticians (Wahba and Wendelberger 1980).

These initial studies clearly nurtured the discussion of the estimation of observation **R**, model **Q**, or background **B** = **P**^f error covariance matrices in the modern DA literature. For demonstration purposes, the algorithms proposed in Dee (1995) and Desroziers and Ivanov (2001) were tested on realistic DA problems, using a shallow-water model on a plane with a simplified Kalman filter, and using the French ARPEGE three-dimensional variational framework, respectively. In both cases, although good performances have been obtained with a small number of iterations, the proposed algorithms have shown some limits, in particular with regard to the simultaneous estimation of the two sources of errors: observation and model (or background). In this context, Todling (2015) pointed out that using only the current innovation is not enough to distinguish the impact of **Q** and **R**, which still makes their simultaneous estimation challenging. Given that our preliminary focus here is to review methods for the joint estimate of **Q**
<table><caption>TABLE 1. Comparison of several methods to estimate error covariance **Q** and **R** in data assimilation.</caption><thead><tr><th>Estimation method</th><th>Criteria</th><th>Estimation of covariance <b>Q</b></th><th>Suitable for non-Gaussian errors</th><th>Application to the highest complexity model</th></tr></thead><tbody><tr><td>Method of moments</td><td>Innovation statistics in the observation space</td><td>No (inflation of <b>P</b>/<i>f</i> instead)</td><td>No</td><td>NWP</td></tr><tr><td>Method of moments</td><td>Lag innovation between consecutive times</td><td>Yes</td><td>No</td><td>Lorenz-96</td></tr><tr><td>Likelihood methods</td><td>Bayesian update of the posterior distribution</td><td>No (or joint parameter with <b>R</b>)</td><td>Yes (using particle filters and not EnKF)</td><td>Shallow water</td></tr><tr><td>Likelihood methods</td><td>Maximization of the total likelihood</td><td>Yes</td><td>Yes (using particle filters and not EnKF)</td><td>Two-scale Lorenz-96</td></tr></tbody></table>

and **R**, the work Dee (1995) and Desroziers and Ivanov (2001) are not further detailed hereinafter. After these two seminal studies, various alternatives were proposed. They are based on the use of several types of innovations and are discussed in this review.

### d. Methods presented in this review

The main topic of this review is the “joint estimation of **Q** and **R**.” Thus, only methods based on this specific goal are presented in detail. A history of what have been, in our opinion, the most relevant contributions and the key milestones for **Q** and **R** covariance estimation in DA is sketched in Fig. 3. The highlighted papers are discussed in this review, with a summary of the different methods, given in Table 1. We distinguish four methods and we can classify them into two categories: those that rely on moment-based methods, and those using likelihood-based methods. Both methods make use of the innovations. The main concepts of the techniques are briefly introduced below.

On the one hand, moment-based methods assume equality between theoretical and empirical statistical moments. A first approach is to study different type of innovations in the observation space (i.e., working in the space of the observations instead of the space of the state). It has been initiated in DA by Rutherford (1972) and Hollingsworth and Lönnberg (1986). A second approach extracts information from the correlation between lag innovations, namely innovations between consecutive times. On the other hand, likelihood-based methods aim to maximize likelihood functions with statistical algorithms. One option is to use a Bayesian framework, assuming prior distributions for the parameters of **Q** and **R** covariance matrices. Another option is to use the iterative expectation-maximization algorithm to maximize a likelihood function.

The four methods listed in Fig. 3 will be examined in this paper. Before doing that, it is worth mentioning existing review work that have attempted to summarize the methods in DA context and beyond.

### e. Other review papers

Other review papers on parameter estimation (including **Q** and **R** matrices) in state-space models have appeared in the statistical and signal processing communities. The first one (Mehra 1972) introduces moment- and likelihood-based methods in the linear and Gaussian case [i.e., when $η$ and $ε$ are Gaussians and $M$ is a linear operator in Eqs. (1) and (2)]. Many extensions to nonlinear state-space models have been proposed since the seminal work of Mehra, and these studies are summarized in the recent review by Duník et al. (2017), with a focus on moment-based methods and the extended Kalman filter (Jazwinski 1970). The book chapter by Buehner (2010) presents another review of moment-based methods, with a focus on the modeling and estimation of spatial covariance structures **Q** and **R** in DA with the ensemble Kalman filter algorithm (Evensen 2009).

In the statistical community, the recent development of powerful simulation techniques, known as sequential Monte Carlo algorithms or particle filters, has led to an extensive literature on the statistical inference in nonlinear state-space models relying on likelihood-based approaches. A recent and detailed presentation of this literature can be found in Kantas et al. (2015). However, these methods typically require a large number of particles, which make them impractical for geophysical DA applications.

The review presented here focuses on methods proposed in DA, especially the moment- and likelihood-based techniques that are suitable for geophysical systems (i.e., with high dimensionality and strong nonlinearities).

### f. Structure of this review

The paper is organized as follows. Section 2 briefly presents the filtering and smoothing DA algorithms used in this work. The main families of methods used in the literature to jointly estimate error covariance matrices **Q** and **R** are then described. First, moment-based methods are introduced in section 3. Then, we describe in section 4 the likelihood-based methods. We also mention other alternatives in section 5, along with
methods used in the past but not exactly matching the scope of this review, and diagnostic tools to check the accuracy of Q and R. In section 6, we provide a summary and discussion on what we consider to be the forthcoming challenges in this area.

## 2. Filtering and smoothing algorithms

This review paper focuses on the estimation of Q and R in the context of ensemble-based DA methods. For the overall discussion of the methods and to set the notation, a short description of the ensemble version of the Kalman recursions is presented in this section: the ensemble Kalman filter (EnKF) and ensemble Kalman smoother (EnKS).

The EnKF and EnKS estimate various state vectors $\mathbf{x}^f(k)$, $\mathbf{x}^a(k)$, $\mathbf{x}^s(k)$ and covariance matrices $\mathbf{P}^f(k)$, $\mathbf{P}^a(k)$, $\mathbf{P}^s(k)$, at each time step $1 \le k \le K$, where K represents the total number of assimilation steps. Kalman-based algorithms assume a Gaussian prior distribution

$$p[\mathbf{x}(k)|\mathbf{y}(1:k-1)] \sim \mathcal{N}[\mathbf{x}^f(k), \mathbf{P}^f(k)].$$

Then, filtering and smoothing estimates correspond to the Gaussian posterior distributions

$$p[\mathbf{x}(k)|\mathbf{y}(1:k)] \sim \mathcal{N}[\mathbf{x}^a(k), \mathbf{P}^a(k)] \quad \text{and}$$

$$p[\mathbf{x}(k)|\mathbf{y}(1:K)] \sim \mathcal{N}[\mathbf{x}^s(k), \mathbf{P}^s(k)]$$

of the state conditionally to past/present observations and past/present/future observations, respectively.

The basic idea of the EnKF and EnKS is to use an ensemble $\mathbf{x}_1, \dots, \mathbf{x}_{N_e}$ of size $N_e$ to track Gaussian distributions over time with the empirical mean vector $\bar{\mathbf{x}} = (1/N_e) \sum_{i=1}^{N_e} \mathbf{x}_i$ and the empirical error covariance matrix $[1/(N_e - 1)] \sum_{i=1}^{N_e} (\mathbf{x}_i - \bar{\mathbf{x}})(\mathbf{x}_i - \bar{\mathbf{x}})^T$.

The EnKF/EnKS equations are divided into three main steps, $\forall i = 1, \dots, N_e$ and $\forall k = 1, \dots, K$—the forecast step (forward in time):

$$\mathbf{x}_i^f(k) = \mathcal{N}_k[\mathbf{x}_i^a(k-1)] + \eta_i(k); \quad (4a)$$

the analysis step (forward in time):

$$\mathbf{d}_i(k) = \mathbf{y}(k) - \mathcal{H}_k[\mathbf{x}_i^f(k)] + \epsilon_i(k), \quad (4b)$$

$$\mathbf{K}^f(k) = \mathbf{P}^f(k) \mathcal{H}_k^T [\mathcal{H}_k \mathbf{P}^f(k) \mathcal{H}_k^T + \mathbf{R}(k)]^{-1}, \quad \text{and} \quad (4c)$$

$$\mathbf{x}_i^a(k) = \mathbf{x}_i^f(k) + \mathbf{K}^f(k) \mathbf{d}_i(k); \quad (4d)$$

and the reanalysis step (backward in time):

$$\mathbf{K}^s(k) = \mathbf{P}^a(k) \mathcal{H}_k^T [\mathbf{P}^f(k+1)]^{-1} \quad \text{and} \quad (4e)$$

$$\mathbf{x}_i^s(k) = \mathbf{x}_i^a(k) + \mathbf{K}^s(k)[\mathbf{x}_i^s(k+1) - \mathbf{x}_i^f(k+1)], \quad (4f)$$

with $\mathbf{K}^f(k)$ and $\mathbf{K}^s(k)$ being the filter and smoother Kalman gains, respectively. Here, $\mathbf{P}^f(k)$ and $\mathcal{H}_k \mathbf{P}^f(k) \mathcal{H}_k^T$ denote the empirical covariance matrices of $\mathbf{x}_i^f(k)$ and $\mathcal{H}_k[\mathbf{x}_i^f(k)]$, respectively. Then, $\mathbf{P}^f(k) \mathcal{H}_k^T$ and $\mathbf{P}^a(k) \mathcal{H}_k^T$ denote the empirical cross-covariance matrices between $\mathbf{x}_i^f(k)$ and $\mathcal{H}_k[\mathbf{x}_i^f(k)]$ and between $\mathbf{x}_i^a(k)$ and $\mathcal{H}_k[\mathbf{x}_i^a(k)]$, respectively. These quantities are estimated using $N_e$ ensemble members.

In some of the methods presented in this review, the ensembles are also used to approximate $\mathcal{H}_k$ and $\mathcal{H}_k$ by linear operators $\mathbf{M}_k$ and $\mathbf{H}_k$ such as

$$\mathbf{M}_k = \mathbf{E}_k^{\mathcal{N}(a)} (\mathbf{E}_{k-1}^a)^{\dagger} \quad \text{and} \quad (5a)$$

$$\mathbf{H}_k = \mathbf{E}_k^{\mathcal{N}(f)} (\mathbf{E}_k^f)^{\dagger}, \quad (5b)$$

with the dagger indicating the pseudoinverse and $\mathbf{E}_k^{\mathcal{N}(a)}$, $\mathbf{E}_{k-1}^a$, $\mathbf{E}_k^{\mathcal{N}(f)}$, and $\mathbf{E}_k^f$ being the matrices containing along their columns the ensemble perturbation vectors (the centered ensemble vectors) of $\mathcal{H}_k[\mathbf{x}_i^a(k-1)]$, $\mathbf{x}_i^a(k-1)$, $\mathcal{H}_k[\mathbf{x}_i^f(k)]$, and $\mathbf{x}_i^f(k)$, respectively.

In Eq. (4b), the innovation is denoted as $\mathbf{d}$ and is tracked by $\mathbf{d}_i(k)$, $\dots$, $\mathbf{d}_{N_e}(k)$. The innovation is the key ingredient of the methods presented in sections 3 and 4.

## 3. Moment-based methods

To constrain the model and observational errors in DA systems, initial efforts were focused on the statistics of relevant variables that could contain information on covariances. The innovation, given in Eq. (4b), corresponds to the difference between the observations and the forecast in the observation space. This variable implicitly takes into account the Q and R covariances. Unfortunately, as explained in Blanchet et al. (1997), by using only current observations, their individual contributions cannot be easily disentangled. Thus, the techniques with only the classic innovation $\mathbf{y}(k) - \mathcal{H}_k[\mathbf{x}^f(k)]$ are not discussed further in this review.

Two main approaches have been proposed in the literature to address this issue. They are based on the idea of producing multiple equations involving Q and R. The first approach uses different type of innovation statistics (i.e., not only the classic one). The second approach is based on lag innovations, or differences between consecutive innovations. From a statistical point of view, they refer to the “methods of moments,” where we construct a system of equations that links various moments of the innovations with the parameters and then replace the theoretical moments by the empirical ones in these equations.
a. Innovation statistics in the observation space

This first approach, based on the Desroziers di-
agnostic (Desroziers et al. 2005), is historical and
now popular in the DA community. It does not ex-
actly fit the topic of this review paper (i.e., estimat-
ing the model error **Q**), since it is based on the
inflation of the background covariance matrix **P**^f.
However, this forecast error covariance is defined by
**P**^f(k) = **M**_k**P**^a(k - 1)**M**_k^T + **Q** in the Kalman filter, con-
sidering a linear model operator **M**_k. Thus, even if DA
systems do not use an explicit model error perturba-
tion controlled by **Q**, the inflation of the background
covariance matrix **P**^f has similar effects, compensating
for the lack of an explicit model uncertainty.

Desroziers et al. (2005) proposed examining various innovation statistics in the observation space. It is based on a different type of innovation statistics between observations, forecasts, and analysis, with all of them defined in the observation space: namely, **d**^{o-f}(k) = **y**(k) - H_k[**x**^f(k)] as in Eq. (4b) and **d**^{o-a}(k) = **y**(k) - H_k[**x**^a(k)]. In theory, in the linear and Gaussian case, for unbiased forecast and observation, and when **P**^f(k) and **R**(k) are correctly specified, the Desroziers innovation statistics should verify the equalities:

$$
E[\mathbf{d}^{o-f}(k)\mathbf{d}^{o-f}(k)^{\mathrm{T}}] = \mathbf{H}_k \mathbf{P}^f(k) \mathbf{H}_k^{\mathrm{T}} + \mathbf{R}(k) \quad \text{and} \quad (6a)
$$

$$
E[\mathbf{d}^{o-a}(k)\mathbf{d}^{o-f}(k)^{\mathrm{T}}] = \mathbf{R}(k), \qquad (6b)
$$

with *E* being the expectation operator. Equation (6a) is
given by using Eq. (4b):

$$
\begin{equation}
\begin{split}
\mathbf{d}^{o-f}(k)\mathbf{d}^{o-f}(k)^{\mathrm{T}} ={}& -\mathbf{y}(k)\mathbf{x}^{f}(k)^{\mathrm{T}}\mathbf{H}_{k}^{\mathrm{T}} - \mathbf{H}_{k}\mathbf{x}^{f}(k)\mathbf{y}(k)^{\mathrm{T}} \\
& + \mathbf{H}_{k}\mathbf{x}^{f}(k)\mathbf{x}^{f}(k)^{\mathrm{T}}\mathbf{H}_{k}^{\mathrm{T}} + \mathbf{y}(k)\mathbf{y}(k)^{\mathrm{T}}
\end{split}
\tag{7}
\end{equation}
$$

and then applying the expectation operator and using
the definition of **P**^f and **R**. The observation-minus-
forecast innovation statistics in Eq. (6a) is not useful to
constrain model error **Q**. Indeed, **d**^o^f does not depend
explicitly on **Q** but rather on the forecast error covari-
ance matrix **P**^f. Thus, the combination of Eqs. (6a) and
(6b) can be used as a diagnosis of the forecast and ob-
servational error covariances in the system. A mismatch
between the Desroziers statistics and the actual covari-
ances, namely the left- and right-hand side terms in Eqs.
(6a) and (6b), indicates inappropriate estimated co-
variances **P**^f(k) and **R**(k).

The forecast covariance **P**^f is sometimes badly esti-
mated in ensemble-based assimilation systems. The
limitations may be attributed to a number of causes. The
limited number of ensemble members produces an over

or, most of the time, underestimation of the forecast
variance. Another limitation is the inaccuracies in
methods used to sample initial condition or model
error. The underestimation of the forecast covariance
produces negative feedback, and the estimated anal-
ysis covariance **P**^a is thus underestimated, which in
turn produces a further underestimation of the fore-
cast covariance in the next cycle. This feedback pro-
cess leads to filter divergence, as was pointed out by
Pham et al. (1998), Anderson and Anderson (1999) or
Anderson (2007). To avoid this filter divergence, in-
flating the forecast covariance **P**^f has been proposed.
This covariance inflation accounts for both sampling
errors and the lack of representation of model errors,
like a too-small amplitude for **Q** or the fact that a bias is
omitted in **η** and **ε** from Eqs. (1) and (2). In this context,
the diagnostics given by the Desroziers innovation sta-
tistics have been proposed as a tool to constrain the
required covariance inflation in the system.

We distinguish three inflation methods: multiplicative, additive and relaxation-to-prior. In the multiplicative case, the forecast error covariance matrix **P**^f is usually multiplied by a scalar coefficient greater than 1 (Anderson and Anderson 1999). Using innovation statistics in the observation space, adaptive procedures to estimate this coefficient have been proposed by Wang and Bishop (2003) and Anderson (2007, 2009) conditionally to the spatial location, Li et al. (2009), Miyoshi (2011), Bocquet (2011), Bocquet and Sakov (2012), Miyoshi et al. (2013), Bocquet et al. (2015), El Gharamti (2018) and Raanes et al. (2019). To prevent excessive inflation or deflation, some authors have proposed assuming a priori distribution for the multiplicative inflation factor. The most usual a priori distributions used by the authors are Gaussian in Anderson (2009), inversegamma in El Gharamti (2018) or inverse chi-square in Raanes et al. (2019).

In practice, multiplicative inflation tends to exces-
sively inflate in the data-sparse regions and inflate too
little in the densely observed regions. As a result, the
spread looks like exaggeration of data density (i.e., too
much spread in sparsely observed regions, and vice
versa). Additive inflation solves this problem but re-
quires many samples for additive noise; these drawbacks
and benefits are discussed in Miyoshi et al. (2010). In the
additive inflation case, the diagonal terms of the forecast
and analysis empirical covariance matrices is increased
(Mitchell and Houtekamer 2000; Corazza et al. 2003;
Whitaker et al. 2008; Houtekamer et al. 2009). This
regularization also avoids the problems corresponding
to the inversion of the covariance matrices.

The last alternative is the relaxation-to-prior method.
In application, this technique is more efficient than both
additive and multiplicative inflations because it maintains a reasonable spread structure. The idea is to relax the reduction of the spread at analysis. We distinguish the method proposed in Zhang et al. (2004), where the forecast and analysis ensemble perturbations are blended, from the one given in Whitaker and Hamill (2012), which multiplies the analysis ensemble without blending perturbations. This last method is thus a multiplicative inflation, but applied after the analysis, not the forecast. Ying and Zhang (2015) and Kotsuki et al. (2017b) proposed methods to adaptively estimate the relaxation parameters using innovation statistics. Their conclusions are that adaptive procedures for relaxation-to-prior methods are robust to sudden changes in the observing networks and observation error settings.

Closely connected to multiplicative inflation estimation is statistical modeling of the error variance terms proposed by Bishop and Satterfield (2013) and Bishop et al. (2013). From numerical evidence based on the 10-dimensional Lorenz-96 model, the authors assume an inverse-gamma prior distribution for these variances. This distribution allows for an analytic Bayesian update of the variances using the innovations. Building on Bocquet (2011), Bocquet et al. (2015), and Ménétrier and Auligné (2015), this technique was extended in Satterfield et al. (2018) to adaptively tune a mixing ratio between the true and sample variances.

Adaptive covariance inflations are estimation methods directly attached to a traditional filtering method (such as the EnKF used here), with almost negligible overhead computational cost. In practice, the use of this technique does not necessarily imply an additive error term $\eta$ in Eq. (1). Thus, it is not a direct estimation of $Q$ but rather an inflation applied to $P^f$ in order to compensate for model uncertainties and sampling errors in the EnKFs, as explained in Raanes et al. (2019, their section 4 and appendix C). Several DA systems work with an inflation method and use it for its simplicity, low cost, and efficiency. As an example of inflation techniques, the most straightforward inflation estimation is a multiplicative factor $\lambda$ of the incorrectly scaled $\tilde{P}^f(k)$ so that the corrected forecast covariance is given by $P^f(k) = \lambda(k)\tilde{P}^f(k)$. The estimate of the inflation factor is given by taking the trace of Eq. (6a):

$$ \tilde{\lambda}(k) = \frac{\mathbf{d}^{o-f}(k)^{\mathrm{T}}\mathbf{d}^{o-f}(k) - \mathrm{Tr}[\mathbf{R}(k)]}{\mathrm{Tr}[\mathbf{H}_k\tilde{\mathbf{P}}^f(k)\mathbf{H}_k^{\mathrm{T}}]} \quad (8) $$

The estimated inflation parameter $\tilde{\lambda}$ computed at each time $k$ can be noisy. The use of temporal smoothing of the form $\lambda(k+1) = \rho\tilde{\lambda}(k) + (1-\rho)\lambda(k)$ is crucial in

operational procedures. Alternatively, Miyoshi (2011) proposed calculating the estimated variance of $\lambda(k)$, denoted as $\sigma_{\lambda(k)}^2$, using the central limit theorem. Then, $\lambda(k+1)$ is updated using the previous estimate $\lambda(k)$ and the Gaussian distribution with mean $\tilde{\lambda}(k)$ and variance $\sigma_{\tilde{\lambda}(k)}^2$. From the Desroziers diagnostics, at each time step $k$ and when sufficient observations are available, an estimate of $\mathbf{R}(k)$ is possible using Eq. (6b). For instance, Li et al. (2009) proposed estimating each component of a diagonal and averaged $\mathbf{R}$ matrix. However, the diagonal terms cannot take into account spatial correlated error terms, and constant values for observation errors are not realistic. Then, Miyoshi et al. (2013) proposed additionally estimating the off-diagonal components of the time-dependent matrix $\mathbf{R}(k)$. The Miyoshi et al. (2013) implementation is summarized in the appendix as algorithm 1.

The Desroziers diagnostic method has been applied widely to estimate the real observation error covariance matrix $\mathbf{R}$ in numerical weather prediction (NWP). The observations are coming from different sources. In the case of satellite radiances, Bormann et al. (2010) applied three methods, including the Desroziers diagnostic and the method detailed in Hollingsworth and Lönnberg (1986) to estimate a constant diagonal term of $\mathbf{R}$ using the innovation $\mathbf{d}^{o-f}$ and its correlations in space, assuming that horizontal correlations in $\mathbf{d}^{o-f}$ samples are purely due to $\mathbf{P}^f$. Weston et al. (2014) and Campbell et al. (2017) then included the interchannel observation error correlations of satellite radiances in DA and obtained improved results when compared with the case using a diagonal $\mathbf{R}$. For spatial error correlations in $\mathbf{R}$, Kotsuki et al. (2017a) estimated the horizontal observation error correlations of satellite-derived precipitation data. Including horizontal observation error correlations in DA for densely observed data from satellites and radars is more challenging than including interchannel error correlations in DA. Indeed, the number of horizontally error-correlated observations is much larger, and some recent studies have been tackling this issue (e.g., Guillet et al. 2019).

To conclude, the Desroziers diagnostic is a consistency check and makes it possible to detect if the error covariances $\mathbf{P}^f$ and $\mathbf{R}$ are incorrect. When and how this method can result in accurate or inaccurate estimates, and convergence properties, have been studied in depth by Waller et al. (2016) and Ménard (2016). The Desroziers diagnostic is also useful to estimate off-diagonal terms of $\mathbf{R}$, for instance taking into account the spatial error correlations. However, covariance localization used in the ensemble Kalman filter might induce erroneous estimates of spatial correlations (Waller et al. 2017).
b. Lag innovation between consecutive times

Another way to estimate error covariances is to use multiple equations involving **Q** and **R**, exploiting cross correlations between lag innovations. More precisely, it involves the current innovation **d**(k) = **d**<sup>o-f</sup>(k) defined in Eq. (4b) and past innovations **d**(k - 1), ..., **d**(k - l). Lag innovations were introduced by Mehra (1970) to recover **Q** and **R** simultaneously for Gaussian, linear and stationary dynamic systems. In such a case, {**d**(k)}<sub>k≥1</sub> is completely characterized by the lagged covariance matrix **C**<sub>l</sub> = Cov[**d**(k), **d**(k - l)], which is independent of k. In other words, the information encoded in {**d**(k)}<sub>k≥1</sub> is completely equivalent to the information provided by {**C**<sub>l</sub>}<sub>l≥0</sub>. Moreover, for linear systems in a steady state, analytic relations exist between **Q**, **R** and E[**d**(k) **d**(k - l)<sup>T</sup>]. However, these linear relations can be dependent and redundant for different lags l. Therefore, as stated in Mehra (1970), only a limited number of **Q** components can be recovered.

Bélanger (1974) extended these results to the case of time-varying linear stochastic processes, taking **d**(k) **d**(k - l)<sup>T</sup> as “observations” of **Q** and **R** and using a secondary Kalman filter to update them iteratively. On the one hand, considering the time-varying case may increase the number of components in **Q** that can be estimated. On the other hand, as pointed out in Bélanger (1974), this method would no longer be analytically ex- act if **Q** and **R** were updated adaptively at each time step. One numerical difficulty of Bélanger’s method is that it needs to invert a matrix of size m<sup>2</sup> × m<sup>2</sup>, where m refers to the dimension of the observation vector. However, this difficulty has been largely overcome by Dee et al. (1985) in which the matrix inversion is reduced to O(m<sup>3</sup>) by taking the advantage of the fact that the big matrix comes from some tensor product.

More recent work has focused on high-dimensional and nonlinear systems using the extended or ensemble Kalman filters. Berry and Sauer (2013) proposed a fast and adaptive algorithm inspired by the use of lag innovations proposed by Mehra. Harlim et al. (2014) applied the original Bélanger algorithm empirically to a nonlinear system with sparse observations. Zhen and Harlim (2015) proposed a modified version of Bélanger’s method, by removing the secondary filter and alternatively solving **Q** and **R** in a least squares sense based on the averaged linear relation over a long term.

Here, we briefly describe the algorithm of Berry and Sauer (2013), considering the lag-zero and lag-one innovations. The following equations are satisfied in the linear and Gaussian case, for unbiased forecast and observation when **P**<sup>f</sup>(k) and **R**(k) are correctly specified:

$$
E[\mathbf{d}(k)\mathbf{d}(k)^{\mathrm{T}}] = \mathbf{H}_k \mathbf{P}^f(k) \mathbf{H}_k^{\mathrm{T}} + \mathbf{R}(k) = \mathbf{\Sigma}(k) \quad \text{and} \quad (9a)
$$

$$
\begin{equation}
\begin{split}
E[\mathbf{d}(k)\mathbf{d}(k-1)^{\mathrm{T}}] ={}& \mathbf{H}_k \mathbf{M}_k \mathbf{P}^f(k-1) \mathbf{H}_{k-1}^{\mathrm{T}} \\
& - \mathbf{H}_k \mathbf{M}_k \mathbf{K}^f(k-1) \mathbf{\Sigma}(k-1).
\end{split}
\tag{9b}
\end{equation}
$$

Equation (9a) is equivalent to Eq. (6a). Moreover, Eq. (9b) results from the fact that, developing the expression of **d**(k) using consecutively Eqs. (2), (1), (4a), and (4d), the innovation can be written as

$$
\begin{align}
\mathbf{d}(k) &= \mathbf{y}(k) - \mathbf{H}_k \mathbf{x}^f(k) \nonumber \\
&= \mathbf{H}_k [\mathbf{x}(k) - \mathbf{x}^f(k)] + \boldsymbol{\epsilon}(k) \nonumber \\
&= \mathbf{H}_k [\mathbf{M}_k \mathbf{x}(k-1) - \mathbf{x}^f(k) + \boldsymbol{\eta}(k)] + \boldsymbol{\epsilon}(k) \nonumber \\
&= \mathbf{H}_k \{\mathbf{M}_k [\mathbf{x}(k-1) - \mathbf{x}^a(k-1)] + \boldsymbol{\eta}(k)\} + \boldsymbol{\epsilon}(k) \nonumber \\
&= \mathbf{H}_k \mathbf{M}_k [\mathbf{x}(k-1) - \mathbf{x}^f(k-1) - \mathbf{K}^f(k-1) \mathbf{d}(k-1)] \nonumber \\
&\qquad + \mathbf{H}_k \boldsymbol{\eta}(k) + \boldsymbol{\epsilon}(k). \tag{10}
\end{align}
$$

Hence, the innovation product **d**(k) **d**(k - 1)<sup>T</sup> between
two consecutive times is given by

$$
\begin{equation}
\begin{split}
& \mathbf{H}_k \mathbf{M}_k [\mathbf{x}(k-1) - \mathbf{x}^f(k-1)] \mathbf{d}(k-1)^\mathrm{T} \\
& \quad - \mathbf{H}_k \mathbf{M}_k [\mathbf{K}^f(k-1) \mathbf{d}(k-1)] \mathbf{d}(k-1)^\mathrm{T} \quad (11) \\
& \quad + \mathbf{H}_k \boldsymbol{\eta}(k) \mathbf{d}(k-1)^\mathrm{T} + \boldsymbol{\epsilon}(k) \mathbf{d}(k-1)^\mathrm{T}
\end{split}
\end{equation}
$$

and, assuming that the model **η** and observation **ε** error
noises are white and mutually uncorrelated, then
E[**η**(*k*)**d**(*k* − 1)<sup>T</sup>] = 0 and E[**ε**(*k*)**d**(*k* − 1)<sup>T</sup>] = 0. Last,
developing E[**d**(*k*)**d**(*k* − 1)<sup>T</sup>], Eq. (9b) is satisfied.

The algorithm in Berry and Sauer (2013) is sum-
marized in the appendix as algorithm 2. It is based
on an adaptive estimation of **Q**(k) and **R**(k), which
satisfies the following relations in the linear and
Gaussian case:

$$
\begin{equation}
\begin{split}
\tilde{\mathbf{P}}(k) ={}& (\mathbf{H}_k \mathbf{M}_k)^{-1} \mathbf{d}(k) \mathbf{d}(k-1)^{\mathrm{T}} \mathbf{H}_{k-1}^{\mathrm{T}}, \\
& + \mathbf{K}^f(k-1) \mathbf{d}(k-1) \mathbf{d}(k-1)^{\mathrm{T}} \mathbf{H}_{k-1}^{\mathrm{T}},
\end{split}
\tag{12a}
\end{equation}
$$

$$
\tilde{\mathbf{Q}}(k) = \tilde{\mathbf{P}}(k) - \mathbf{M}_{k-1} \mathbf{P}^a(k-2) \mathbf{M}_{k-1}^{\mathrm{T}}, \quad \text{and} \quad (12\text{b})
$$

$$
\tilde{\mathbf{R}}(k) = \mathbf{d}(k)\mathbf{d}(k)^{\mathrm{T}} - \mathbf{H}_k \mathbf{P}^f(k)\mathbf{H}_k^{\mathrm{T}}. \quad (12c)
$$

In operational applications, when the number of ob-
servations is not equal to the number of components in
state **x**, **H** is not a square matrix and Eq. (12a) is ill-
defined. To avoid the inversion of **H**, Berry and Sauer
(2013) proposed considering parametric models for **Q**
and then solving a linear system associated with Eqs.
(12a) and (12b). It is written as a least squares problem
such that
$$ \tilde{\mathbf{Q}}(k) = \arg \min_{\mathbf{Q}} \|\mathbf{d}(k)\mathbf{d}(k-1)^{\mathrm{T}} + \mathbf{H}_k \mathbf{M}_k \mathbf{K}^f(k-1)\mathbf{d}(k-1)\mathbf{d}(k-1)^{\mathrm{T}} \\ - \mathbf{H}_k \mathbf{M}_k \mathbf{M}_{k-1} \mathbf{P}^a(k-2) \mathbf{M}_{k-1}^{\mathrm{T}} \mathbf{H}_{k-1}^{\mathrm{T}} - \mathbf{H}_k \mathbf{M}_k \mathbf{Q} \mathbf{H}_{k-1}^{\mathrm{T}} \| \quad (13) $$

In this adaptive procedure, joint estimations of $\tilde{\mathbf{Q}}(k)$ and $\tilde{\mathbf{R}}(k)$ can abruptly vary over time. Thus, the temporal smoothing of the covariances being estimated becomes crucial. As suggested by Berry and Sauer (2013), such temporal smoothing between current and past estimates is a reasonable choice:

$$ \mathbf{Q}(k+1) = \rho \tilde{\mathbf{Q}}(k) + (1-\rho) \mathbf{Q}(k) \quad \text{and} \quad (14a) $$

$$ \mathbf{R}(k+1) = \rho \tilde{\mathbf{R}}(k) + (1-\rho) \mathbf{R}(k), \quad (14b) $$

with $\mathbf{Q}(1)$ and $\mathbf{R}(1)$ being the initial conditions and $\rho$ being the smoothing parameter. When $\rho$ is large (close to 1), weight is given to the current estimates $\mathbf{Q}$ and $\tilde{\mathbf{R}}$, and when $\rho$ is small (close to 0) it gives smoother $\mathbf{Q}$ and $\mathbf{R}$ sequences. The value of $\rho$ is arbitrary and may depend on the system and how it is observed. For instance, in the case where the number of observations equals the size of the system, Berry and Sauer (2013) uses $\rho = 5 \times 10^{-5}$ in order to estimate the full matrix $\mathbf{Q}$ for the Lorenz-96 model.

The algorithm in Berry and Sauer (2013) only considers lag-zero and lag-one innovations. By incorporating more lags, Zhen and Harlim (2015) and Harlim (2018) showed that it makes it possible to deal with the case in which some components of $\mathbf{Q}$ are not identifiable from the method in Berry and Sauer (2013). For instance, let us consider the two-dimensional system with any stationary operator $\mathbf{M}$ and $\mathbf{H} = [1, 0]$, meaning that only the first component of the system is observed. This is a linear, Gaussian, stationary system, and Mehra's theory implies that two parameters of $\mathbf{Q}$ are identifiable. However, using only lag-one innovations as in Berry and Sauer (2013), Eq. (13) becomes a scalar equation and only one parameter of $\mathbf{Q}$ can be determined. The idea of considering more lag innovations to estimate more components of $\mathbf{Q}$ was tested in Zhen and Harlim (2015). Numerical results show that considering more than one lag can improve the estimates of $\mathbf{Q}$ and $\mathbf{R}$. For instance, Zhen and Harlim (2015) focused on the Lorenz-96 model. Results show that when $\mathbf{Q}$ is stationary, the trace of $\mathbf{Q}$ and $\mathbf{R}$ are equal, and when observations are taken at twenty fixed equally spaced grid points for every five integration time steps, the optimal RMSE of the estimates of $\mathbf{Q}$ and $\mathbf{R}$ is achieved when four time lags are considered. But with more lags, the performance is degraded.

To summarize, methods based on lag innovation between consecutive times have been studied for a long

time in the signal processing community. The original methods (Mehra 1970; Bélanger 1974) were analytically established for linear systems with Gaussian noises. Inspired by these foundational ideas, empirical methods have been established for nonlinear systems in DA (Berry and Sauer 2013; Harlim et al. 2014; Zhen and Harlim 2015). Although these methods have not been tested in any operational experiment, the idea of using lagged innovations seems to have significant potential.

## 4. Likelihood-based methods

This section focuses on methods based on the likelihood of the observations, given a set of statistical parameters. The conceptual idea behind what we refer to as likelihood-based methods is to determine the optimal statistical parameters (i.e., $\mathbf{Q}$ and $\mathbf{R}$) that maximize the likelihood function for a given set of observations that may be distributed over time. In this way, the aim is to derive estimation methods that use the observations to find the most suitable, or most likely parameters.

Early studies in Dee (1995), Blanchet et al. (1997), Mitchell and Houtekamer (2000) and Liang et al. (2012) proposed finding the optimal $\mathbf{Q}$ and $\mathbf{R}$ that maximize the current innovation likelihood at time $k$. Unfortunately, if only the current observations are used, the joint estimation of $\mathbf{Q}$ and $\mathbf{R}$ is not well constrained (Todling 2015). To tackle this issue, several solutions have been recently proposed where the likelihood function considers observations distributed in time over several assimilation cycles.

The likelihood-based methods are broadly divided into two categories. One approach uses a Bayesian framework. It assumes a priori knowledge about the parameters and estimate jointly the posterior distribution of $\mathbf{Q}$ and $\mathbf{R}$ together with the state of the system, or alternatively to estimate them in a two-stage process.<fnref n="2" /> The second one is based on the frequentist viewpoint and attempts a point estimate of the parameters by maximizing a total likelihood function.

<footnote n="2">² Some of the methods presented in section 3 also use the Bayesian philosophy; for instance, they assume a priori distribution for the multiplicative inflation parameter $\lambda$ (Anderson 2009; El Gharamti 2018).</footnote>
a. Bayesian inference

In the Bayesian framework, the elements of the co-variance matrices **Q** and **R** are assumed to have a priori distributions that are controlled by hyperparameters. In practice, it is difficult to have prior distributions for each element of **Q** and **R**, especially for large DA systems. Instead, parametric forms are used for the matrices, typically describing the shape and level noise. We denote the corresponding parameters as **θ**.

The inference in the Bayesian framework aims to determine the posterior density $p[\theta|\mathbf{y}(1:k)]$. Two techniques have appeared, the first based on a state augmentation and the second based on a rigorous Bayesian update of the posterior distribution.

1) STATE AUGMENTATION

In the Bayesian framework, **θ** is a random variable such that the state is augmented with these parameters by defining **z**(k) = [$\mathbf{x}(k), \theta$]. To define an augmented state-space model, one has to define an evolution equation for the parameters. This leads to a new state-space model of the form of Eqs. (1) and (2) with **x** replaced by **z**. Therefore, the state and the parameters are estimated jointly using the DA algorithms.

State augmentation was first proposed in Schmidt (1966) and is known as the Schmidt-Kalman filter. This technique was mainly used to estimate both the state of the system and additional parameters, including bias, forcing terms and physical parameters. These kinds of parameters are strongly related to the state of the system (Ruiz et al. 2013a). Therefore, they are identifiable and suitable for an augmented state approach. However, Stroud and Bengtsson (2007) and later DeISole and Yang (2010) formally demonstrated that augmentation methods fail for variance parameters like **Q** and **R**. The explanation is that in the EnKF, the empirical forecast covariance **P**^f is computed using all the ensemble members, each one with a different realization of the random variable **θ**. Thus, **P**^f and consequently the Kalman gain **K**^f, are mixing the effects of **Q** and **R** parameters contained in **θ**. Therefore, after applying Eq. (4d), the update of **z** corresponding to the **θ** parameters is the same for all of the parameters. To capture the impact of a single variance parameter on the prediction covariance and circumvent the limitation of the state augmentation, Scheffler et al. (2019) proposed to use an ensemble of states integrated with the same variance parameter. The choice of an ensemble of states for each variance parameter leads to two nested ensemble Kalman filters. The technique performs successfully under different model error covariance structures but has an important computational cost.

Another critical aspect of state augmentation is that one needs to define an evolution model for the augmented state **z**(k) = [$\mathbf{x}(k), \theta(k)$]. If persistence is assumed in the parameters such that they are constant in time, this leads to filter degeneracy, since the estimated variance of the error in **θ** is bound to decrease in time. To prevent or at least mitigate this issue, it was suggested to use an independent inflation factor on the parameters (Ruiz et al. 2013b) or to impose artificial stochastic dynamics for **θ**, typically a random walk or AR(1) model, as introduced in Eq. (3) and proposed in Liu and West (2001). The tuning of the parameters introduced in these artificial dynamics may be difficult, and this introduces bias into the procedure, which is hard to quantify.

2) BAYESIAN UPDATE OF THE POSTERIOR DISTRIBUTION

Instead of the inference of the joint posterior density using a state augmentation strategy, the state **x**(k) and parameters **θ** can be divided into a two-step inference procedure using the following formula:

$$p[\mathbf{x}(k), \theta|\mathbf{y}(1:k)] = p[\mathbf{x}(k)|\mathbf{y}(1:k), \theta]p[\theta|\mathbf{y}(1:k)], \quad (15)$$

which is a direct consequence of the conditional density definition. In Eq. (15), $p[\mathbf{x}(k)|\mathbf{y}(1:k), \theta]$ represents the posterior distribution of the state, given the observations and the parameter **θ**. It can be computed using a filtering DA algorithm. The second term on the right-hand side of Eq. (15) corresponds to the posterior distribution of the parameters, given the observations up to time k. The latter can be updated sequentially using the following Bayesian hierarchy:

$$p[\theta|\mathbf{y}(1:k)] \propto p[\mathbf{y}(k)|\mathbf{y}(1:k-1), \theta]p[\theta|\mathbf{y}(1:k-1)], \quad (16)$$

where $p[\mathbf{y}(k)|\mathbf{y}(1:k-1), \theta]$ is the likelihood of the innovations.

Different approximations have been used for $p[\theta|\mathbf{y}(1:k)]$ in Eq. (16); these include parametric models based on Gaussian (Stroud et al. 2018), inverse-gamma (Stroud and Bengtsson 2007) or Wishart distributions (Ueno and Nakamura 2016), particle-based approximations (Frei and Künsch 2012; Stroud et al. 2018) and grid-based approximation (Stroud et al. 2018).

The methods proposed in the literature also differ by the approximation used for the likelihood of the innovations. We emphasize that $p[\mathbf{y}(k)|\mathbf{y}(1:k-1), \theta]$ needs to be evaluated for different values of **θ** at each time step, and that this requires applying the filter from the initial time with a single value of **θ**, which is computationally impossible for applications in high dimensions.
To reduce computational time, it is generally assumed that $\mathbf{x}^f$ and $\mathbf{P}^f$ are independent of $\theta$, and only observations $\mathbf{y}(k-l:k-1)$ in a small time window from the current observation are used when computing the likelihood of the innovations [see Ueno and Nakamura (2016) and Stroud et al. (2018) for a more detailed discussion]. A summary of the Bayesian method from Stroud et al. (2018) is given in the appendix as algorithm 3. It was implemented within the EnKF framework and is one of the most recent studies based on the Bayesian approach.

Applications of the Bayesian method in the DA context are now discussed. It has mainly been used to estimate shape and noise parameters of **Q** and **R** error covariance matrices. For instance, Purser and Parrish (2003) and Solonen et al. (2014) estimated statistical parameters controlling the magnitude of the variance and the spatial dependencies in the model error **Q**, assuming that **R** is known. There are also applications aimed at estimating parameters governing the shape of the observation error covariance matrix **R** only: Frei and Künsch (2012) and Stroud et al. (2018) in the Lorenz-96 system, Winiarek et al. (2012, 2014) for the inversion of the source term of airborne radionuclides using a regional atmospheric model, and Ueno and Nakamura (2016) using a shallow-water model to assimilate satellite altimetry.

As pointed out in Stroud and Bengtsson (2007), Bayesian update algorithms work best when the number of unknown parameters in $\theta$ is small. This limitation may explain why the joint estimation of parameters controlling both model and observation error covariances is not systematically addressed. For instance, Stroud and Bengtsson (2007) used the EnKF with the Lorenz-96 model for the estimation of a common multiplicative scalar parameter for pre-defined matrices **Q** and **R**. Alternatively, Stroud et al. (2018) tested the Bayesian method on different spatiotemporal systems to estimate the signal-to-noise ratio between **Q** and **R**. Nevertheless, based on the experiments about the importance of the signal-to-noise ratio $||\mathbf{P}^f||/||\mathbf{R}||$ presented in Fig. 2, we know that this estimation of the ratio is not optimal.

Widely used in the statistical community, the Bayesian framework is useful incorporating physical knowledge about error covariance matrices and constraining their estimation process. In the DA literature, authors have used a priori distributions for the shape and noise parameters of **Q** and **R**, but rarely both. Operationally, only a limited number of parameters can be estimated. To address this issue, Stroud and Bengtsson (2007) suggested combining Bayesian algorithms with other techniques.

b. Maximization of the total likelihood

The innovation likelihood at time $k$, $p[\mathbf{y}(k)|\mathbf{y}(1:k-1), \theta]$ in Eq. (16), can be maximized to find the optimal $\theta$ (i.e., **Q** and **R** matrices or parameterizations of them). In practice, when this maximization is done at each time step, two issues arise. First, the innovation covariance matrix $\mathbf{\Sigma}(k) = \mathbf{H}_k \mathbf{P}^f(k) \mathbf{H}_k^T + \mathbf{R}(k)$ combines the information about **R** and **Q**, the latter being contained in $\mathbf{P}^f$. When using only time $k$, it is difficult to disentangle the model and observation error covariances; in application, the aforementioned studies only estimated one of them. Second, the number of observations at each time step is in general limited and, as pointed out by Dee (1995), available observations should exceed “the number of tunable parameters by two or three orders of magnitude.” To overcome these limitations, a reasonable alternative is to use a batch of observations within a time window and to assume $\theta$ to be constant in time. The resulting total likelihood expressed sequentially through conditioning is given by

$$p[\mathbf{y}(1:K)|\theta] = \prod_{k=1}^{K} p[\mathbf{y}(k)|\mathbf{y}(1:k-1), \theta]. \quad (17)$$

Because it is an integration of innovation likelihoods over a long period of time from $k=1$ to $k=K$, Eq. (17) provides more observational information to estimate **Q** and **R**. The maximization of this total likelihood has been applied for the estimation of deterministic and stochastic parameters (related to **Q**) using a direct sequential optimization procedure (DelSole and Yang 2010). Ueno et al. (2010) used a grid-based procedure to estimate noise levels and spatial correlation lengths of **Q** and a noise level for **R**. This grid-based method uses predefined sets of covariance parameters and evaluates the different combinations to find the one that maximizes the likelihood criterion. Brankart et al. (2010) also proposed a method using the same criterion but adding (at the initial time) information on scale and correlation length parameters of **Q** and **R**. This information is only given the first time, and is progressively forgotten over time, using a decreasing exponential factor. The marginalization of the hidden state in Eq. (17) considers all the previous observations, and it requires the use of a filter. The maximization of the total likelihood $p[\mathbf{y}(1:K)|\theta]$ to estimate model error covariance **Q** was conducted in Pulido et al. (2018), where they used a gradient-based optimization technique and the EnKF.

The likelihood function given in Eq. (17) only depends on the observations $\mathbf{y}$. This likelihood can be
written in a different way, taking into account both the observations and the hidden state **x**. Indeed, the marginalization of the hidden state to obtain the total likelihood can be produced using the whole trajectory of the state from *k* = 0 to the last time step *K* all at once. It is given by

$$p[\mathbf{y}(1:K)|\mathbf{\theta}] = \int p[\mathbf{x}(0:K), \mathbf{y}(1:K)|\mathbf{\theta}] d\mathbf{x}(0:K). \quad (18)$$

The maximization of the total likelihood as a function of statistical parameters **θ** is not possible, since the total likelihood cannot be evaluated directly, nor its gradient with regard to the parameters (Pulido et al. 2018). Shumway and Stoffer (1982) proposed using an iterative procedure based on the expectation-maximization algorithm (hereinafter denoted as EM). They applied it to estimate the parameters of a linear state-space model, with linear dynamics, and a linear observational operator and Gaussian errors. The EM algorithm was introduced by Dempster et al. (1977).

Each iteration of the EM algorithm consists of two steps. In the expectation step (E-step), the posterior density $p[\mathbf{x}(0:K)|\mathbf{y}(1:K), \mathbf{\theta}_{(n)}]$ is determined conditioned on the batch of observations $\mathbf{y}(1:K)$ and given the parameters $\mathbf{\theta}_{(n)} = [\mathbf{Q}_{(n)}, \mathbf{R}_{(n)}]$ from the previous iteration or initial guess. This is obtained through the application of a smoother like the EnKS. Then, the M-step relies on the maximization of an intermediate function, depending on the posterior density obtained in the E-step. The intermediate function is defined by the conditional expectation:

$$E(\log\{p[\mathbf{x}(0:K), \mathbf{y}(1:K)|\mathbf{\theta}]\}|\mathbf{y}(1:K), \mathbf{\theta}_{(n)}). \quad (19)$$

If, as in Eqs. (1) and (2), the observational and model errors are assumed to be additive, unbiased, and Gaussian, the expression for the logarithm of the joint density in Eq. (19) is given by

$$-\frac{1}{2}\left\{\sum_{k=1}^{K} \|\mathbf{x}(k) - \mathcal{M}[\mathbf{x}(k-1)]\|_{\mathbf{Q}}^{2} + \log|\mathbf{Q}| + \|\mathbf{y}(k) - \mathcal{H}[\mathbf{x}(k)]\|_{\mathbf{R}}^{2} + \log|\mathbf{R}|\right\} + c, \quad (20)$$

where $||\mathbf{v}||_A^2$ is defined to be equal to $\mathbf{v}^T\mathbf{A}^{-1}\mathbf{v}$ and $c$ is a constant independent of **Q** and **R**. In this case, an analytic expression for the optimal error covariances at each iteration of the EM algorithm can be obtained. The estimators of the parameters that maximize Eq. (19) using Eq. (20) are

$$\mathbf{Q}_{(n+1)} = \frac{1}{K} \sum_{k=1}^{K} E\left(\{\mathbf{x}(k) - \mathcal{M}[\mathbf{x}(k-1)]\} \{\mathbf{x}(k) - \mathcal{M}[\mathbf{x}(k-1)]\}^{\mathrm{T}} | \mathbf{y}(1:K), \mathbf{\theta}_{(n)}\right) \quad \text{and} \quad (21a)$$

$$\mathbf{R}_{(n+1)} = \frac{1}{K} \sum_{k=1}^{K} E\left(\{\mathbf{y}(k) - \mathcal{H}[\mathbf{x}(k)]\}\{\mathbf{y}(k) - \mathcal{H}[\mathbf{x}(k)]\}^{\mathrm{T}}|\mathbf{y}(1:K), \mathbf{\theta}_{(n)}\right). \quad (21b)$$

The application of the EM algorithm for the estimation of **Q** and **R** is straightforward. Starting from $\mathbf{Q}_{(1)}$ and $\mathbf{R}_{(1)}$, an ensemble Kalman smoother is applied with this first guess and the batch of observations $\mathbf{y}(1:K)$ to obtain the posterior density $p[(\mathbf{x}(0:K)|\mathbf{y}(1:K), \mathbf{\theta}_{(1)}]$. Then Eqs. (21a) and (21b) are used to update and obtain $\mathbf{Q}_{(2)}$ and $\mathbf{R}_{(2)}$. Next, a new application of the smoother is conducted using the parameters $\mathbf{Q}_{(2)}$ and $\mathbf{R}_{(2)}$ and the observations $\mathbf{y}(1:K)$, the new resulting states are used in Eqs. (21a) and (21b) to estimate $\mathbf{Q}_{(3)}$ and $\mathbf{R}_{(3)}$, and so on. As a diagnostic of convergence or as a stop criterion, the product of innovation likelihood functions given in Eq. (17) is evaluated using a filter. The EM algorithm guarantees that the total likelihood increases in each iteration and that the sequence $\mathbf{\theta}_{(n)}$ converges to a local maximum (Wu 1983). A summary of the EM method (using EnKF and EnKS) from Dreano et al. (2017) is given in the appendix as algorithm 4.

EM is a well-known algorithm used in the statistical community. This procedure is parameter-free and robust, due to the large number of observations used to approximate the likelihood when using a long batch period (Shumway and Stoffer 1982). Although the use of the EM algorithm is still limited in DA, it is becoming more and more popular. Some studies have implemented the EM algorithm for estimating only the observation error matrix **R**. For instance, Ueno and Nakamura (2014) used the model proposed in Zebiak and Cane (1987) and satellite altimetry observations, whereas Liu et al. (2017) used an air quality model for accidental pollutant source retrieval. But the estimation of only the observation error covariance is limited, and other studies have tried to jointly estimate model error **Q** and **R** matrices, for instance as in Tandeo et al. (2015) for an orographic subgrid-scale nonlinear observation operator. Then, Dreano et al. (2017) and Pulido et al. (2018) used the EM procedure to produce joint estimation of **Q** and **R** matrices in the Lorenz-63 and stochastic parameters of the Lorenz-96 systems, respectively. Recently, Yang and Mémin (2019) extended the EM procedure for the estimation of physical
parameters in a one-dimensional shallow-water model, more specifically for the identification of stochastic subgrid terms. Last, an online adaptation of the EM algorithm for the estimation of **Q** and **R** at each time step, after the filtering procedure, has been proposed in Cocucci et al. (2020). In this adaptive case, the likelihood is averaged locally over time, see Cappé (2011) for more details.

To our knowledge, EM has not been tested yet on operational systems with large observation and state space. In that case, the use of parametric forms for the matrices **Q** and **R** is essential to reduce the number of statistical parameters **θ** to estimate. For instance, Dreano et al. (2017) and Liu et al. (2017) showed that in the particular cases where covariances are diagonal or of the form **αA**, with **A** being a positive definite matrix, expressions in Eqs. (21a) and (21b) are simplified, and a sub-optimal **θ** in the space of the parametric covariance form can be obtained.

## 5. Other methods

In this section, we describe other methods that have been used to estimate **Q** and **R**, and that cannot be included in the categories presented in the previous sections. In particular, we report here about methods that are applied either a posteriori, after DA cycles, or without applying any DA algorithms.

### a. Analysis (or reanalysis) increment approach

This first method is based on previous DA outputs. The key idea here is to use the analysis (or reanalysis) increments to provide a realistic sample of model errors from which statistical moments, such as the covariance matrix **Q**, can be empirically estimated. This assumes that the sequence of reanalysis **x**<sup>s</sup> (or analysis **x**<sup>a</sup>) is the best available representation of the true process **x**. In that case, the following approximation in Eq. (1) is made:

$$ \eta(k) = \mathcal{N}[\mathbf{x}(k-1)] - \mathbf{x}(k) \\ \approx \mathcal{N}[\mathbf{x}^s(k-1)] - \mathbf{x}^s(k). \qquad (22) $$

In this approximation, it is implicitly assumed that the estimated state is the truth so that the initial condition at time $k-1$ is neglected. A similar approximation of the true process by **x**<sup>a</sup> or **x**<sup>s</sup> in Eq. (2) can be used to estimate the observation error covariance matrix **R**.

Operationally, the analysis (or reanalysis) increment method is applied after a DA filter (or smoother) to estimate the **Q** matrix. This method was originally introduced by Leith (1978), and later used

to account for model error in the context of ensemble Kalman filters, using analysis and reanalysis increments by Mitchell and Carrassi (2015), and in the context of weak-constraint variational assimilation by Bowler (2017). Along this line, Rodwell and Palmer (2007) also proposed evaluating the average of instantaneous analysis increments to represent the systematic forecast tendencies of a model.

### b. Covariance matching

The covariance matching method was introduced by Fu et al. (1993). It involves matching sample covariance matrices to their theoretical expectations. Thus, it is a method of moments, similar to the work in Desroziers et al. (2005), except that covariance matching is performed on a set of historical observations and numerical simulations (noted **x**<sup>sim</sup>), without applying any DA algorithms. It has been extended by Menemenlis and Chechelnitsky (2000) to time-lagged innovations, as first considered in Bélanger (1974).

In the case of a constant and linear observation operator **H**, the basic idea in Fu et al. (1993) is to assume the following system:

$$ \mathbf{x}^{\text{sim}}(k) = \mathbf{x}(k) + \eta^{\text{sim}}(k), \qquad (23a) $$

$$ \eta^{\text{sim}}(k) = \mathbf{A}\eta^{\text{sim}}(k-1) + \eta(k), \quad \text{and} \qquad (23b) $$

$$ \mathbf{H}\mathbf{x}^{\text{sim}}(k) - \mathbf{y}(k) = \mathbf{H}\eta^{\text{sim}}(k) + \varepsilon(k), \qquad (23c) $$

with **A** being a transition matrix close to the identity matrix, assuming slow variations of the numerical simulation errors $\eta^{\text{sim}}$. In Eqs. (23b) and (23c), the definitions of $\eta$ and $\varepsilon$ errors remain similar, as in the general Eqs. (1) and (2).

Assuming that **Q** and **R** are constant over time, $\varepsilon$ is uncorrelated from **x** and from $\eta^{\text{sim}}$, then Eqs. (23c) and (23a) yield to the following estimates of **R** and $\mathbf{P}^{\text{sim}}$ (the latter represents the error covariance of the numerical simulations):

$$ \hat{\mathbf{R}} = \frac{1}{2} \left\{ \begin{aligned} & E\left[(\mathbf{y} - \mathbf{H}\mathbf{x}^{\text{sim}})(\mathbf{y} - \mathbf{H}\mathbf{x}^{\text{sim}})^{\text{T}}\right] \\ & - E\left[(\mathbf{H}\mathbf{x}^{\text{sim}})(\mathbf{H}\mathbf{x}^{\text{sim}})^{\text{T}}\right] + E(\mathbf{y}\mathbf{y}^{\text{T}}) \end{aligned} \right\} \quad \text{and} \qquad (24a) $$

$$ \mathbf{H}\hat{\mathbf{P}}^{\text{sim}}\mathbf{H}^{\text{T}} = \frac{1}{2} \left\{ \begin{aligned} & E\left[(\mathbf{y} - \mathbf{H}\mathbf{x}^{\text{sim}})(\mathbf{y} - \mathbf{H}\mathbf{x}^{\text{sim}})^{\text{T}}\right] \\ & + E\left[(\mathbf{H}\mathbf{x}^{\text{sim}})(\mathbf{H}\mathbf{x}^{\text{sim}})^{\text{T}}\right] - E(\mathbf{y}\mathbf{y}^{\text{T}}) \end{aligned} \right\}, \qquad (24b) $$

where E is the expectation operator over time. Then, an estimate of **Q** is obtained using Eqs. (23b) and
(24b) and assuming that $\mathbf{P}^{\text{sim}}$ has a unique time-invariant limit.

### c. Forecast sensitivity

In operational meteorology, it is critical to learn the sensitivity of the forecast accuracy to various parameters of a DA system, in particular the error statistics of both the model and the observations. This is why a significant portion of literature considers the tuning problem of **R** and **Q** through the lens of the sensitivity of the forecast to these parameters. The computation of those sensitivities can be seen as a first-order correction or diagnostic for such an estimation. The forecast sensitivities are computed either using the adjoint model (Daescu and Todling 2010; Daescu and Langland 2013) in the context of variational methods, or a forecast ensemble (Hotta et al. 2017) in the context of the EnKF.

The basic idea is to compute at each assimilation cycle an innovation between forecast and analysis, noted $\mathbf{d}^{l-^a}(k) = \mathbf{x}^l(k) - \mathbf{x}^a(k)$. Then, the forecast sensitivity is given by $\mathbf{d}^{l-^a}(k)^\mathrm{T}\mathbf{S}\mathbf{d}^{l-^a}(k)$ with **S** being a diagonal scaling matrix, to normalize the components of $\mathbf{d}^{l-^a}$. The **Q** and **R** estimates are the matrices that minimize $\mathbf{d}^{l-^a}(k)$. The adjoint or the ensemble are thus used to compute the partial derivatives of this forecast sensitivity. w.r.t. **Q** and **R**.

## 6. Conclusions and perspectives

As often considered in data assimilation, this review paper also deals with model and observation errors that are assumed additive and Gaussian with covariance matrices **Q** and **R**. The model error corresponds to the dynamic model deficiencies to represent the underlying physics, whereas the observation error corresponds to the instrumental noise and the representativity error. Model and observation errors are assumed to be uncorrelated and white in time. The model and observations are also assumed unbiased, a strong assumption for real data assimilation applications.

The discussion starts with the aid of an illustration of the individual and joint impacts of improperly calibrated covariances using a linear toy model. The experiments clearly showed that to achieve reasonable filter accuracy (i.e., in terms of root-mean-squared error), it is crucial to carefully define both **Q** and **R**. The effect on the coverage probability of a misspecification of **Q** and **R** is also highlighted. This coverage probability is related to the estimated covariance of the reconstructed state, and thus to the uncertainty quantification in data assimilation. After

the one-dimensional illustration, the core of the paper gives an overview of various methods to jointly estimate the **Q** and **R** error covariance matrices: they are summarized and compared below.

### a. Comparison of existing methods for estimating **Q** and **R**

We mainly focused in this review on four methods for the joint estimation of the error covariances **Q** and **R**. The methods are summarized in Table 1. They correspond to classic estimation methods, based on statistical moments or likelihoods. The main difference between the four methods comes from the innovations taken into account: the total innovation, as in the EM algorithm proposed by Shumway and Stoffer (1982); lag innovations, following the idea given in Mehra (1970); or different type of innovations in the observation space, as in Desroziers et al. (2005). Additionally, to constrain the estimation, hierarchical Bayesian approaches use prior distributions for the shape parameters of **Q** and **R**.

Most of the methods estimate the model error **Q**. The exception is the one using the Desroziers diagnostic, dealing with different type of innovations in the observation space, which instead estimates an inflation factor for $\mathbf{P}^l$. Moreover, the methods are mainly defined online, meaning that they aim to estimate **Q** and **R** adaptively, together with the current state of the system. Consequently, these methods require additional tunable parameters to smooth the estimated covariances over time. However, most of the methods presented in this review also have an offline variant. In that case, a batch of observations is used to estimate **Q** and **R**. In some methods, such as the EM algorithm, the parameters are determined iteratively. These offline approaches avoid the use of additional smoothing parameters.

Throughout this review paper, as usually stated in DA, it is assumed that model error $\eta$ and observation error $\epsilon$, defined in Eqs. (1) and (2), are Gaussian. Consequently, the distribution of the innovations is also Gaussian. The four presented methods use this property to build estimates of **Q** and **R** adequately. But, if $\eta$ and $\epsilon$ are non-Gaussian, Desroziers diagnostic and lag-innovation methods are not suitable anymore. However, the EM procedures and Bayesian methods are still relevant, although they must be used with an appropriate filter (e.g., particle filters), not Kalman-based algorithms (i.e., assuming a Gaussian distribution of the state). Recently, the treatment of non-Gaussian error distributions in DA has been explored in Katzfuss et al. (2020), using hierarchical state-space models. This Bayesian framework allows
to handle unknown variables that cannot be easily included in the state vector (e.g., parameters of **Q** and **R**) and to model non-Gaussian observations.

The four methods have been applied at different levels of complexity. For instance, Bayesian inference methods (due to their algorithm complexity) and the EM algorithm (due to its computational cost) have so far only been applied to small dynamic models. However, the online version of the EM algorithm is less consuming and opens new perspectives of applications on larger models. On the other hand, methods using innovation statistics in the observation space have already been applied to NWP models.

The four methods summarized in Table 1 show differences in maturity in terms of applications and methodological aspects. This review also shows that there are still remaining challenges and possible improvements for the four methods.

## b. Remaining challenges for each method

The first challenge concerns the improvements of adaptive techniques regarding additional parameters that control the variations of **Q** and **R** estimates over time. Instead of using fixed values for these parameters, for instance fixed $\rho$ in the lag innovations or $\sigma_{\lambda}^2$ in the inflation methods, we suggest using time-dependent adaptations. This adaptive solution could avoid the problems of instabilities close to the solution. Another option could be to adapt these procedures, working with stable parameter values (small $\rho$, low $\sigma_{\lambda}^2$) and iterating the procedures on a batch of observations, as in the EM algorithm. This offline variant was suggested and tested in Desroziers et al. (2005) with encouraging results. To the best of our knowledge, it has not yet been tested with lag-innovation methods.

The second challenge concerns considering time-varying error covariance matrices. The adaptive procedures, based on online estimations with temporal smoothing of **Q** and **R**, are supposed to capture slowly evolving covariances. On the contrary, offline methods like the EM algorithm are working on a batch of observations, assuming that **Q** and **R** are constant over the batch period. Online solutions for the EM algorithm, with the likelihood averaged locally over time (Cocucci et al. 2020), could also capture slow evolution of the covariances. Another simple solution could be to work on small sets of observations, named as minibatches, and to apply the EM algorithm in each set using the previous estimates as an initial guess. These intermediate schemes are of common use in machine learning.

A third challenge has to do with the assumption, used by all of the methods described herein, that observation and model errors are mutually independent. Nevertheless,

as pointed out in Berry and Sauer (2018), observation and model error are often correlated in real data assimilation problems (e.g., for satellite retrieval of Earth observations that uses model outputs in the inversion process). Methods based on Bayesian inference can, in principle, exploit existing model-to-observation correlations by using a prior joint distribution (i.e., not two individual ones). The explicit taking into account of this correlation can then constrain the optimization procedure. This is not possible in the other approaches described in this review, at least not in their standard known formulations, and the presence of model-observation correlation can deteriorate their accuracy.

A fourth challenge is common to all the methods presented in this review. Iterative versions of the presented algorithms need initial values or distributions for **R** and **Q** (or $\mathbf{B} = \mathbf{P}^f$ in the case of Desroziers). However, as mentioned in Waller et al. (2016) for the Desroziers diagnostics, there is no guarantee that the algorithms will converge to the optimal solution. Indeed, in such an optimization problem, there are possibly several local and nonoptimal solutions. Suboptimal specifications of **R**, **Q**, or **B** in the initial DA cycle will affect the final estimation results. There are several solutions to avoid this convergence problem: initialize the covariance matrices using physical expertise, execute the iterative algorithms several times with different initial covariance matrices, or use stochastic perturbations in the optimization algorithms to avoid to be trapped in local solutions. These aspects of convergence and sensitivity to initial conditions have so far been poorly addressed. It is therefore necessary to check which method is robust operationally.

The last remaining challenge concerns the estimation of other statistical parameters of the state-space model given in Eqs. (1) and (2) and associated filters. Indeed, the initial conditions $\mathbf{x}(0)$ and $\mathbf{P}(0)$ are crucial for certain satellite retrieval problems and have to be estimated. This is the case, for instance, when the time sequence of observations is short (i.e., shorter than the spinup time of the filter with an uninformative prior) or when filtering and smoothing are repeated on various iterations, as in the EM algorithm. Estimation methods should also consider the estimation of systematic or time-varying biases, the deterministic part of $\eta$ and $\epsilon$. This was initially proposed by Dee and da Silva (1999) and tested in Dee et al. (1999) in the case of maximizing the innovation likelihood, in Dee (2005) in a state augmentation formulation, and was adapted to a Bayesian update formulation in Liu et al. (2017) and in Berry and Harlim (2017). Recently, the joint estimation of bias and covariance error terms, for
the treatment of brightness temperatures from the European geostationary satellite, has been successfully applied in Merchant et al. (2020).

### c. Perspectives for geophysical DA

Beyond the aforementioned potential improvements in the existing techniques, specific research directions need to be taken by the data assimilation community. The main one concerns the realization of a comprehensive numerical evaluation of the different methods for the estimation of **Q** and **R**, built on an agreed experimental framework and a consensus model. Such an effort would help to evaluate (i) the pros and cons of the different methods (including their capability to deal with high dimensionality, localization in ensemble methods, and their practical feasibility), (ii) their effects on different error statistics (RMSE, coverage probabilities, and other diagnostics), (iii) the potential combination of the various methods (especially those considering constant or adaptive covariances), and (iv) the capability to take into account other sources of error (due for instance to improper parameterizations, multiplicative errors, or forcing terms).

The use of a realistic DA problem, with a high-dimensional state-space and a limited and heterogeneous observational coverage should be addressed in the future. In that realistic case, the observational information per degree of freedom will be significantly lower, and the estimates of **Q** and **R** will deteriorate. Parametric versions of these error covariance matrices will therefore be necessary. Among the parameters, some of them will control the variances, and will be different depending on the variable. Other parameters will control the spatial correlation lengths, that could be isotropic or anisotropic, depending on the region of interest and the considered variable. Cross correlations between variables will also have to be considered. Consequently, **Q** and **R** will be block matrices with as few parameters as possible.

A further challenge for future work is the evaluation of the feasibility of estimating nonadditive, non-Gaussian, and time-correlated noises under the current estimation frameworks. In this way, the need for observational constraints for the stochastic perturbation methods in the NWP community could be considered within the estimation framework discussed in this review.

*Acknowledgments.* This work has been carried out as part of the Copernicus Marine Environment Monitoring Service (CMEMS) 3DA project. CMEMS is implemented by Mercator Ocean in the framework of a delegation agreement with the European Union. This work

was also partially supported by FOCUS Establishing Supercomputing Center of Excellence. CEREA is a member of Institut Pierre Simon Laplace (IPSL). Author Carrassi has been funded by the project REDDA (250711) of the Norwegian Research Council. Carrassi was also supported by the Natural Environment Research Council (Agreement PR140015 between NERC and the National Centre for Earth Observation). We thank Paul Platzer, a second-year Ph.D. student, who helped to popularize the summary and the introduction, and John C. Wells, Gilles-Olivier Guégan, and Aimée Johansen for their English grammar corrections. We also thank the five anonymous reviewers for their precious comments and ideas to improve this review paper. We are immensely grateful to Prof. David M. Schultz, Chief Editor of *Monthly Weather Review*, for his detailed advice and careful reading of the paper.

# APPENDIX

## Four Main Algorithms to Jointly Estimate **Q** and **R** in Data Assimilation

Algorithm 1 is an adaptive algorithm for the EnKF as implemented by Miyoshi et al. (2013). The steps of the algorithm are the following:

- initialize inflation factor [for instance $\lambda(1) = 1$]);

**for** *k* in 1:K **do**

**for** *i* in 1:N<sub>e</sub> **do**

- compute forecast **x**<sub>*i*</sub><sup>*f*</sup>(*k*) using Eq. (4a);

- compute innovation **d**<sub>*i*</sub>(*k*) using Eq. (4b);

**end**

- compute empirical covariance $\tilde{\mathbf{P}}^f(k)$ of the $\mathbf{x}_i^f(k)$;

- compute $\mathbf{K}^f(k)$ using Eq. (4c) where $\mathbf{P}^f(k) \mathcal{H}_k^T$ and $\mathcal{H}_k \mathbf{P}^f(k) \mathcal{H}_k^T$ are inflated by $\lambda(k)$;

**for** *i* in 1:N<sub>e</sub> **do**

- compute analysis $\mathbf{x}_i^a(k)$ using Eq. (4d);

**end**

- compute mean innovations $\mathbf{d}_i^{o-f}(k)$ and $\mathbf{d}_i^{o-a}(k)$ with $\mathbf{d}_i^{o-f}(k) = \mathbf{y}(k) - \mathcal{H}_k[\mathbf{x}_i^f(k)]$ and $\mathbf{d}_i^{o-a}(k) = \mathbf{y}(k) - \mathcal{H}_k[\mathbf{x}_i^a(k)]$;

- update $\mathbf{R}(k)$ from Eq. (6b) using the cross-covariance between $\mathbf{d}_i^{o-f}(k)$ and $\mathbf{d}_i^{o-a}(k)$;

- estimate $\lambda(k)$ using Eq. (8) where $\mathcal{H}_k \tilde{\mathbf{P}}^f(k) \mathcal{H}_k^T$ is inflated by $\lambda(k)$;

- update $\lambda(k + 1)$ using temporal smoother;

**end**

Algorithm 2 is an adaptive algorithm for the EnKF by Berry and Sauer (2013). The steps of the algorithm are the following:
- initialize $\mathbf{Q}(1)$ and $\mathbf{R}(1)$;

for $k$ in 1:$K$ do

for $i$ in 1:$N_e$ do

- compute forecast $\mathbf{x}_i^f(k)$ using Eq. (4a);

- compute innovation $\mathbf{d}_i(k)$ using Eq. (4b);

end

- compute $\mathbf{K}^f(k)$ using Eq. (4c);

for $i$ in 1:$N_e$ do

- compute analysis $\mathbf{x}_i^a(k)$ using Eq. (4d);

end

- apply Eq. (12a) to get $\tilde{\mathbf{P}}(k)$ using linearizations of $\mathbf{M}_k$ and $\mathbf{H}_k$ given in Eqs. (5a) and (5b);

- estimate $\mathbf{Q}(k)$ using Eq. (12b);

- estimate $\tilde{\mathbf{R}}(k)$ using Eq. (12c);

- update $\mathbf{Q}(k+1)$ and $\mathbf{R}(k+1)$ using temporal smoothers;

end

Algorithm 3 is an adaptive algorithm for the EnKF from Stroud et al. (2018). The steps of the algorithm are the following:

- define a priori distributions for $\theta$ (shape parameters of $\mathbf{Q}$ and $\mathbf{R}$);

for $k$ in 1:$K$ do

do $i$ in 1:$N_e$ do

- draw samples $\theta_i(k)$ from $p[\theta|\mathbf{y}(1:k-1)];$

- compute forecast $\mathbf{x}_i^f(k)$ using Eq. (4a) with $\theta_i(k)$;

- compute innovation $\mathbf{d}_i(k)$ using Eq. (4b) with $\theta_i(k)$;

end

- compute $\mathbf{K}^f(k)$ using Eq. (4c);

for $i$ in 1:$N_e$ do

- compute analysis $\mathbf{x}_i^a(k)$ using Eq. (4d);

end

- approximate Gaussian likelihood of innovations $p[\mathbf{y}(k)|\mathbf{y}(1:k-1), \theta(k)]$ using empirical mean $\bar{\mathbf{d}}(k) = (1/N_e)\sum_{i=1}^{N_e}\mathbf{d}_i(k)$ and empirical covariance $\mathbf{\Sigma}(k) = [1/(N_e-1)]\sum_{i=1}^{N_e}[\mathbf{d}_i(k) - \bar{\mathbf{d}}(k)][\mathbf{d}_i(k) - \bar{\mathbf{d}}(k)]^T$ with $\mathbf{d}_i(k) = \mathbf{y}(k) - \mathcal{H}_k[\mathbf{x}_i^f(k)];$

- update $p[\theta|\mathbf{y}(1:k)]$ using Eq. (16);

end

Algorithm 4 is an EM algorithm for the EnKF/EnKS from Dreano et al. (2017). The steps of the algorithm are the following:

- define $\theta$ (shape parameters of $\mathbf{Q}$ and $\mathbf{R}$);

- set $p[\mathbf{y}(1:\mathbf{K})|\theta_{(0)}] = +\infty$;

- initialize $n = 1$, $\theta_{(1)}$ and $\epsilon$ (stop condition);

while $p[\mathbf{y}(1:\mathbf{K})|\theta_{(n)}] - p[\mathbf{y}(1:\mathbf{K})|\theta_{(n-1)}] > \epsilon$ do

for $k$ in 1:$K$ do

for $i$ in 1:$N_e$ do

- compute forecast $\mathbf{x}_i^f(k)$ using Eq. (4a);

- compute innovation $\mathbf{d}_i(k)$ using Eq. (4b);

end

- compute $\mathbf{K}^f(k)$ using Eq. (4c);

for $i$ in 1:$N_e$ do

- compute analysis $\mathbf{x}_i^a(k)$ using Eq. (4d);

end

end

for $k$ in $K:1$ do

- compute $\mathbf{K}^s(k)$ using Eq. (4e);

for $i$ in 1:$N_e$ do

- compute reanalysis $\mathbf{x}_i^s(k)$ using Eq. (4f);

end

end

- increment $n \leftarrow n + 1$;

- estimate $\mathbf{Q}_{(n)}$ using Eq. (21a);

- estimate $\mathbf{R}_{(n)}$ using Eq. (21b);

end

REFERENCES

Anderson, J. L., 2007: An adaptive covariance inflation error correction algorithm for ensemble filters. *Tellus*, **59A**, 210–224, https://doi.org/10.1111/j.1600-0870.2006.00216.x.

—, 2009: Spatially and temporally varying adaptive covariance inflation for ensemble filters. *Tellus*, **61A**, 72–83, https://doi.org/10.1111/j.1600-0870.2008.00361.x.

—, and S. L. Anderson, 1999: A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts. *Mon. Wea. Rev.*, **127**, 2741–2758, https://doi.org/10.1175/1520-0493(1999)127<2741:AMCIOT>2.0.CO;2.

Bélanger, P. R., 1974: Estimation of noise covariance matrices for a linear time-varying stochastic process. *Automatica*, **10**, 267–275, https://doi.org/10.1016/0005-1098(74)90037-5.

Berry, T., and T. Sauer, 2013: Adaptive ensemble Kalman filtering of non-linear systems. *Tellus*, **65A**, 20331, https://doi.org/10.3402/tellusa.v65i0.20331.

—, and J. Harlim, 2017: Correcting biased observation model error in data assimilation. *Mon. Wea. Rev.*, **145**, 2833–2853, https://doi.org/10.1175/MWR-D-16-0428.1.

—, and T. Sauer, 2018: Correlation between system and observation errors in data assimilation. *Mon. Wea. Rev.*, **146**, 2913–2931, https://doi.org/10.1175/MWR-D-17-0331.1.

Bishop, C. H., and E. A. Satterfield, 2013: Hidden error variance theory. Part I: Exposition and analytic model. *Mon. Wea. Rev.*, **141**, 1454–1468, https://doi.org/10.1175/MWR-D-12-00118.1.

—, —, and K. T. Shanley, 2013: Hidden error variance theory. Part II: An instrument that reveals hidden error variance distributions from ensemble forecasts and observations. *Mon. Wea. Rev.*, **141**, 1469–1483, https://doi.org/10.1175/MWR-D-12-00119.1.

Blanchet, I., C. Frankignoul, and M. A. Cane, 1997: A comparison of adaptive Kalman filters for a tropical Pacific Ocean model. *Mon. Wea. Rev.*, **125**, 40–58, https://doi.org/10.1175/1520-0493(1997)125<0040:ACOAKF>2.0.CO;2.

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

—, and P. Sakov, 2012: Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems. *Nonlinear Processes Geophys.*, **19**, 383–399, https://doi.org/10.5194/npg-19-383-2012.
—, P. N. Raanes, and A. Hannart, 2015: Expanding the validity of the ensemble Kalman filter without the intrinsic need for inflation. *Nonlinear Processes Geophys.*, **22**, 645–662, https://doi.org/10.5194/npg-22-645-2015.

Bormann, N., A. Collard, and P. Bauer, 2010: Estimates of spatial and interchannel observation-error characteristics for current sounder radiances for numerical weather prediction. II: Application to AIRS and IASI data. *Quart. J. Roy. Meteor. Soc.*, **136**, 1051–1063, https://doi.org/10.1002/qj.615.

Bowler, N. E., 2017: On the diagnosis of model error statistics using weak-constraint data assimilation. *Quart. J. Roy. Meteor. Soc.*, **143**, 1916–1928, https://doi.org/10.1002/qj.3051.

Brankart, J.-M., E. Cosme, C.-E. Testut, P. Brasseur, and J. Verron, 2010: Efficient adaptive error parameterizations for square root or ensemble Kalman filters: Application to the control of ocean mesoscale signals. *Mon. Wea. Rev.*, **138**, 932–950, https://doi.org/10.1175/2009MWR3085.1.

Buehner, M., 2010: Error statistics in data assimilation: Estimation and modelling. *Data Assimilation: Making Sense of Observations*, W. Lahoz, B. Khattatov, and R. Menard, Eds., Springer, 93–112.

Campbell, W. F., E. A. Satterfield, B. Ruston, and N. L. Baker, 2017: Accounting for correlated observation error in a dual-formulation 4D variational data assimilation system. *Mon. Wea. Rev.*, **145**, 1019–1032, https://doi.org/10.1175/MWR-D-16-0240.1.

Cappé, O., 2011: Online expectation-maximisation. *Mixtures: Estimation and Applications*, K. L. Mengersen, C. Robert, and M. Titterington, Eds.,Wiley Series in Probability and Statistics, John Wiley & Sons, 1–53.

Carrassi, A., M. Bocquet, L. Bertino, and G. Evensen, 2018: Data assimilation in the geosciences: An overview on methods, issues and perspectives. *Wiley Interdiscip. Rev.: Climate Change*, **9**, e535, https://doi.org/10.1002/wcc.535.

Chapnik, B., G. Desroziers, F. Rabier, and O. Talagrand, 2004: Properties and first application of an error-statistics tuning method in variational assimilation. *Quart. J. Roy. Meteor. Soc.*, **130**, 2253–2275, https://doi.org/10.1256/qj.03.26.

Cocucci, T. J., M. Pulido, M. Lucini, and P. Tandeo, 2020: Model error covariance estimation in particle and ensemble Kalman filters using an online expectation-maximization algorithm. arXiv:2003.02109, https://arxiv.org/pdf/2003.02109.pdf.

Corazza, M., E. Kalnay, D. J. Patil, R. Morss, M. Cai, I. Szunyogh, B. R. Hunt, and J. A. Yorke, 2003: Use of the breeding technique to estimate the structure of the analysis “errors of the day.” *Nonlinear Processes Geophys.*, **10**, 233–243, https://doi.org/10.5194/npg-10-233-2003.

Daescu, D. N., and R. Todling, 2010: Adjoint sensitivity of the model forecast to data assimilation system error covariance parameters. *Quart. J. Roy. Meteor. Soc.*, **136**, 2000–2012, https://doi.org/10.1002/qj.693.

—, and R. H. Langland, 2013: Error covariance sensitivity and impact estimation with adjoint 4D-Var: Theoretical aspects and first applications to NAVDAS-AR. *Quart. J. Roy. Meteor. Soc.*, **139**, 226–241, https://doi.org/10.1002/qj.1943.

Daley, R., 1991: *Atmospheric Data Analysis*. Cambridge University Press, 457 pp.

—, 1992: Estimating model-error covariances for application to atmospheric data assimilation. *Mon. Wea. Rev.*, **120**, 1735–1746, https://doi.org/10.1175/1520-0493(1992)120<1735:EMECFA>2.0.CO;2.

Dee, D. P., 1995: On-line estimation of error covariance parameters for atmospheric data assimilation. *Mon. Wea. Rev.*, **123**,

1128–1145, https://doi.org/10.1175/1520-0493(1995)123<1128:OLEOEC>2.0.CO;2.

—, 2005: Bias and data assimilation. *Quart. J. Roy. Meteor. Soc.*, **131**, 3323–3343, https://doi.org/10.1256/qj.05.137.

—, and A. M. da Silva, 1999: Maximum-likelihood estimation of forecast and observation error covariance parameters. Part I: Methodology. *Mon. Wea. Rev.*, **127**, 1822–1834, https://doi.org/10.1175/1520-0493(1999)127<1822:MLEOFA>2.0.CO;2.

—, S. E. Cohn, A. Dalcher, and M. Ghil, 1985: An efficient algorithm for estimating noise covariances in distributed systems. *IEEE Trans. Autom. Control*, **30**, 1057–1065, https://doi.org/10.1109/TAC.1985.1103837.

—, G. Gaspari, C. Redder, L. Rukhovets, and A. M. da Silva, 1999: Maximum-likelihood estimation of forecast and observation error covariance parameters. Part II: Applications. *Mon. Wea. Rev.*, **127**, 1835–1849, https://doi.org/10.1175/1520-0493(1999)127<1835:MLEOFA>2.0.CO;2.

DelSole, T., and X. Yang, 2010: State and parameter estimation in stochastic dynamical models. *Physica D*, **239**, 1781–1788, https://doi.org/10.1016/j.physd.2010.06.001.

Dempster, A. P., N. M. Laird, and D. B. Rubin, 1977: Maximum likelihood from incomplete data via the EM algorithm. *J. Roy. Stat. Soc.*, **39B**, 1–22, https://doi.org/10.1111/j.2517-6161.1977.tb01600.x.

Desroziers, G., and S. Ivanov, 2001: Diagnosis and adaptive tuning of observation-error parameters in a variational assimilation. *Quart. J. Roy. Meteor. Soc.*, **127**, 1433–1452, https://doi.org/10.1002/qj.49712757417.

—, L. Berre, B. Chapnik, and P. Poli, 2005: Diagnosis of observation, background and analysis-error statistics in observation space. *Quart. J. Roy. Meteor. Soc.*, **131**, 3385–3396, https://doi.org/10.1256/qj.05.108.

Dreano, D., P. Tandeo, M. Pulido, T. Chonavel, B. A. It-El-Fquih, and I. Hoteit, 2017: Estimating model error covariances in nonlinear state-space models using Kalman smoothing and the expectation-maximisation algorithm. *Quart. J. Roy. Meteor. Soc.*, **143**, 1877–1885, https://doi.org/10.1002/qj.3048.

Dunfjyk, J., O. Straka, O. Kost, and J. Havlík, 2017: Noise covariance matrices in state-space models: A survey and comparison of estimation methods-Part I. *Int. J. Adapt. Control Signal Process.*, **31**, 1505–1543, https://doi.org/10.1002/acs.2783.

El Gharamti, M., 2018: Enhanced adaptive inflation algorithm for ensemble filters. *Mon. Wea. Rev.*, **146**, 623–640, https://doi.org/10.1175/MWR-D-17-0187.1.

Evensen, G., 2009: *Data Assimilation: The Ensemble Kalman Filter*. Springer Science and Business Media, 332 pp.

Frei, M., and H. R. Künsch, 2012: Sequential state and observation noise covariance estimation using combined ensemble Kalman and particle filters. *Mon. Wea. Rev.*, **140**, 1476–1495, https://doi.org/10.1175/MWR-D-10-05088.1.

Fu, L.-L., I. Fukumori, and R. N. Miller, 1993: Fitting dynamic models to the Geosat sea level observations in the tropical Pacific Ocean. Part II: A linear, wind-driven model. *J. Phys. Oceanogr.*, **23**, 2162–2181, https://doi.org/10.1175/1520-0485(1993)023<2162:FDMTTG>2.0.CO;2.

Ghil, M., and P. Malanotte-Rizzoli, 1991: Data assimilation in meteorology and oceanography. *Advances in Geophysics*, Vol. 33, Academic Press, 141–266, https://doi.org/10.1016/S0065-2687(08)60442-2.

Guillet, O., A. T. Weaver, X. Vasseur, Y. Michel, S. Gratton, and S. Gürol, 2019: Modelling spatially correlated observation errors in variational data assimilation using a diffusion operator
on an unstructured mesh. *Quart. J. Roy. Meteor. Soc.*, **145**, 1947–1967, https://doi.org/10.1002/qj.3537.

Harlim, J., 2018: Ensemble Kalman filters. *Data-Driven Computational Methods: Parameter and Operator Estimations*, J. Harlim, Ed., Cambridge University Press, 31–59.

—, A. Mahdi, and A. J. Majda, 2014: An ensemble Kalman filter for statistical estimation of physics constrained nonlinear regression models. *J. Comput. Phys.*, **257**, 872–812, https://doi.org/10.1016/j.jcp.2013.10.025.

Hollingsworth, A., and P. Lönnberg, 1986: The statistical structure of short-range forecast errors as determined from radiosonde data. Part I: The wind field. *Tellus*, **38A**, 111–136, https://doi.org/10.3402/tellusa.v38i2.11707.

Hotta, D., E. Kalnay, Y. Ota, and T. Miyoshi, 2017: EFSR: Ensemble forecast sensitivity to observation error covariance. *Mon. Wea. Rev.*, **145**, 5015–5031, https://doi.org/10.1175/MWR-D-17-0122.1.

Houtekamer, P. L., and F. Zhang, 2016: Review of the ensemble Kalman filter for atmospheric data assimilation. *Mon. Wea. Rev.*, **144**, 4489–4532, https://doi.org/10.1175/MWR-D-15-0440.1.

—, H. L. Mitchell, and X. Deng, 2009: Model error representation in an operational ensemble Kalman filter. *Mon. Wea. Rev.*, **137**, 2126–2143, https://doi.org/10.1175/2008MWR2737.1.

Ide, K., P. Courtier, M. Ghil, and A. C. Lorenc, 1997: Unified notation for data assimilation: Operational, sequential and variational. *J. Meteor. Soc. Japan*, **75**, 181–189, https://doi.org/10.2151/jmsj1965.75.1B_181.

Janjić, T., and Coauthors, 2018: On the representation error in data assimilation. *Quart. J. Roy. Meteor. Soc.*, **144**, 1257–1278, https://doi.org/10.1002/qj.3130.

Jazwinski, A. H., 1970: *Stochastic Processes and Filtering Theory*. Academic Press, 376 pp.

Kantas, N., A. Doucet, S. S. Singh, J. Maciejowski, and N. Chopin, 2015: On particle methods for parameter estimation in state-space models. *Stat. Sci.*, **30**, 328–351, https://doi.org/10.1214/14-STS511.

Katzfuss, M., J. R. Stroud, and C. K. Wikle, 2020: Ensemble Kalman methods for high-dimensional hierarchical dynamic space-time models. *J. Amer. Stat. Assoc.*, **115**, 866–885, https://doi.org/10.1080/01621459.2019.1592753.

Kotsuki, S., T. Miyoshi, K. Terasaki, G.-Y. Lien, and E. Kalnay, 2017a: Assimilating the global satellite mapping of precipitation data with the Nonhydrostatic Icosahedral Atmospheric Model (NICAM). *J. Geophys. Res. Atmos.*, **122**, 631–650, https://doi.org/10.1002/2016jd025355.

—, Y. Ota, and T. Miyoshi, 2017b: Adaptive covariance relaxation methods for ensemble data assimilation: Experiments in the real atmosphere. *Quart. J. Roy. Meteor. Soc.*, **143**, 2001–2015, https://doi.org/10.1002/qj.3060.

Leith, C. E., 1978: Objective methods for weather prediction. *Annu. Rev. Fluid Mech.*, **10**, 107–128, https://doi.org/10.1146/annurev.fl.10.010178.000543.

Li, H., E. Kalnay, and T. Miyoshi, 2009: Simultaneous estimation of covariance inflation and observation errors within an ensemble Kalman filter. *Quart. J. Roy. Meteor. Soc.*, **135**, 523–533, https://doi.org/10.1002/qj.371.

Liang, X., X. Zheng, S. Zhang, G. Wu, Y. Dai, and Y. Li, 2012: Maximum likelihood estimation of inflation factors on error covariance matrices for ensemble Kalman filter assimilation. *Quart. J. Roy. Meteor. Soc.*, **138**, 263–273, https://doi.org/10.1002/qj.912.

Liu, J., and M. West, 2001: Combined parameter and state estimation in simulation-based filtering. *Sequential Monte Carlo Methods in Practice*, A. Doucet, N. Freitas, and N. Gordon, Eds., Springer, 197–223.

Liu, Y., J.-M. Haussaire, M. Bocquet, Y. Roustan, O. Saunier, and A. Mathieu, 2017: Uncertainty quantification of pollutant source retrieval: Comparison of Bayesian methods with application to the Chernobyl and Fukushima Daiichi accidental releases of radionuclides. *Quart. J. Roy. Meteor. Soc.*, **143**, 2886–2901, https://doi.org/10.1002/qj.3138.

Mehra, R. K., 1970: On the identification of variances and adaptive Kalman filtering. *IEEE Trans. Autom. Control*, **15**, 175–184, https://doi.org/10.1109/TAC.1970.1099422.

—, 1972: Approaches to adaptive filtering. *IEEE Trans. Autom. Control*, **17**, 693–698, https://doi.org/10.1109/TAC.1972.1100100.

Ménard, R., 2016: Error covariance estimation methods based on analysis residuals: Theoretical foundation and convergence properties derived from simplified observation networks. *Quart. J. Roy. Meteor. Soc.*, **142**, 257–273, https://doi.org/10.1002/qj.2650.

Menemenlis, D., and M. Chechelnitsky, 2000: Error estimates for an ocean general circulation model from altimeter and acoustic tomography data. *Mon. Wea. Rev.*, **128**, 763–778, https://doi.org/10.1175/1520-0493(2000)128<0763:EEFAOG>2.0.CO;2.

Ménétrier, B., and T. Auligné, 2015: Optimized localization and hybridization to filter ensemble-based covariances. *Mon. Wea. Rev.*, **143**, 3931–3947, https://doi.org/10.1175/MWR-D-15-0057.1.

Merchant, C. J., S. Saux-Picart, and J. Waller, 2020: Bias correction and covariance parameters for optimal estimation by exploiting matched in-situ references. *Remote Sens. Environ.*, **237**, 111590, https://doi.org/10.1016/j.rse.2019.111590.

Mitchell, H. L., and P. L. Houtekamer, 2000: An adaptive ensemble Kalman filter. *Mon. Wea. Rev.*, **128**, 416–433, https://doi.org/10.1175/1520-0493(2000)128<0416:AAEKF>2.0.CO;2.

Mitchell, L., and A. Carrassi, 2015: Accounting for model error due to unresolved scales within ensemble Kalman filtering. *Quart. J. Roy. Meteor. Soc.*, **141**, 1417–1428, https://doi.org/10.1002/qj.2451.

Miyoshi, T., 2011: The Gaussian approach to adaptive covariance inflation and its implementation with the local ensemble transform Kalman filter. *Mon. Wea. Rev.*, **139**, 1519–1535, https://doi.org/10.1175/2010MWR3570.1.

—, Y. Sato, and T. Kadowaki, 2010: Ensemble Kalman filter and 4D-Var intercomparison with the Japanese operational global analysis and prediction system. *Mon. Wea. Rev.*, **138**, 2846–2866, https://doi.org/10.1175/2010MWR3209.1.

—, E. Kalnay, and H. Li, 2013: Estimating and including observation-error correlations in data assimilation. *Inverse Probl. Sci. Eng.*, **21**, 387–398, https://doi.org/10.1080/17415977.2012.712527.

Pham, D. T., J. Verron, and M. C. Roubaud, 1998: A singular evolutive extended Kalman filter for data assimilation in oceanography. *J. Mar. Syst.*, **16**, 323–340, https://doi.org/10.1016/S0924-7963(97)00109-7.

Pulido, M., P. Tandeo, M. Bocquet, A. Carrassi, and M. Lucini, 2018: Stochastic parameterization identification using ensemble Kalman filtering combined with maximum likelihood methods. *Tellus*, **70A**, 1–17, https://doi.org/10.1080/16000870.2018.1442099.

Purser, R. J., and D. F. Parrish, 2003: A Bayesian technique for estimating continuously varying statistical parameters of a variational assimilation. *Meteor. Atmos. Phys.*, **82**, 209–226, https://doi.org/10.1007/s00703-001-0583-x.

Raanes, P. N., M. Bocquet, and A. Carrassi, 2019: Adaptive covariance inflation in the ensemble Kalman filter by Gaussian
scale mixtures. Quart. J. Roy. Meteor. Soc., 145, 53–75, https://doi.org/10.1002/qj.3386.

Rodwell, M. J., and T. N. Palmer, 2007: Using numerical weather prediction to assess climate models. Quart. J. Roy. Meteor. Soc., 133, 129–146, https://doi.org/10.1002/qj.23.

Ruiz, J. J., M. Pulido, and T. Miyoshi, 2013a: Estimating model parameters with ensemble-based data assimilation: A review. J. Meteor. Soc. Japan, 91, 79–99, https://doi.org/10.2151/jmsj.2013-201.

—, —, and —, 2013b: Estimating model parameters with ensemble-based data assimilation: Parameter covariance treatment. J. Meteor. Soc. Japan, 91, 453–469, https://doi.org/10.2151/jmsj.2013-403.

Rutherford, I. D., 1972: Data assimilation by statistical interpolation of forecast error fields. J. Atmos. Sci., 29, 809–815, https://doi.org/10.1175/1520-0469(1972)029<0809:DABSIO>2.0.CO;2.

Satterfield, E. A., D. Hodyss, D. D. Kuhl, and C. H. Bishop, 2018: Observation-informed generalized hybrid error covariance models. Mon. Wea. Rev., 146, 3605–3622, https://doi.org/10.1175/MWR-D-18-0016.1.

Scheffler, G., J. Ruiz, and M. Pulido, 2019: Inference of stochastic parametrizations for model error treatment using nested ensemble Kalman filters. Quart. J. Roy. Meteor. Soc., 145, 2028–2045, https://doi.org/10.1002/qj.3542.

Schmidt, S. F., 1966: Applications of state space methods to navigation problems. Adv. Control Syst., 3, 293–340, https://doi.org/10.1016/B978-1-4831-6716-9.50011-4.

Shumway, R. H., and D. S. Stoffer, 1982: An approach to time series smoothing and forecasting using the EM algorithm. J. Time Ser. Anal., 3, 253–264, https://doi.org/10.1111/j.1467-9892.1982.tb00349.x.

Solonen, A., J. Hakkarainen, A. Ilin, M. Abbas, and A. Bibov, 2014: Estimating model error covariance matrix parameters in extended Kalman filtering. *Nonlinear Processes Geophys.*, 21, 919–927, https://doi.org/10.5194/npg-21-919-2014.

Stroud, J. R., and T. Bengtsson, 2007: Sequential state and variance estimation within the ensemble Kalman filter. Mon. Wea. Rev., 135, 3194–3208, https://doi.org/10.1175/MWR3460.1.

—, M. Katzfuss, and C. K. Wikle, 2018: A Bayesian adaptive ensemble Kalman filter for sequential state and parameter estimation. Mon. Wea. Rev., 146, 373–386, https://doi.org/10.1175/MWR-D-16-0427.1.

Tandeo, P., M. Pulido, and F. Lott, 2015: Offline parameter estimation using EnKF and maximum likelihood error covariance estimates: Application to a subgrid-scale orography parametrization. Quart. J. Roy. Meteor. Soc., 141, 383–395, https://doi.org/10.1002/qj.2357.

Todling, R., 2015: A lag-1 smoother approach to system-error estimation: Sequential method. Quart. J. Roy. Meteor. Soc., 141, 1502–1513, https://doi.org/10.1002/qj.2460.

Ueno, G., and N. Nakamura, 2014: Iterative algorithm for maximum-likelihood estimation of the observation-error covariance matrix for ensemble-based filters. Quart. J. Roy. Meteor. Soc., 140, 295–315, https://doi.org/10.1002/qj.2134.

—, and —, 2016: Bayesian estimation of the observation-error covariance matrix in ensemble-based filters. Quart. J. Roy. Meteor. Soc., 142, 2055–2080, https://doi.org/10.1002/qj.2803.

—, T. Higuchi, T. Kagimoto, and N. Hirose, 2010: Maximum likelihood estimation of error covariances in ensemble-based filters and its application to a coupled atmosphere-ocean

model. Quart. J. Roy. Meteor. Soc., 136, 1316–1343, https://doi.org/10.1002/qj.654.

Wahba, G., and J. Wendelberger, 1980: Some new mathematical methods for variational objective analysis using splines and cross validation. Mon. Wea. Rev., 108, 1122–1143, https://doi.org/10.1175/1520-0493(1980)108<1122:SNMMFV>2.0.CO;2.

Waller, J. A., S. L. Dance, and N. K. Nichols, 2016: Theoretical insight into diagnosing observation error correlations using observation-minus-background and observation-minus-analysis statistics. Quart. J. Roy. Meteor. Soc., 142, 418–431, https://doi.org/10.1002/qj.2661.

—, —, and —, 2017: On diagnosing observation-error statistics with local ensemble data assimilation. Quart. J. Roy. Meteor. Soc., 143, 2677–2686, https://doi.org/10.1002/qj.3117.

Wang, X., and C. H. Bishop, 2003: A comparison of breeding and ensemble transform Kalman filter ensemble forecast schemes. J. Atmos. Sci., 60, 1140–1158, https://doi.org/10.1175/1520-0469(2003)060<1140:ACOBAE>2.0.CO;2.

Weston, P. P., W. Bell, and J. R. Eyre, 2014: Accounting for correlated error in the assimilation of high-resolution sounder data. Quart. J. Roy. Meteor. Soc., 140, 2420–2429, https://doi.org/10.1002/qj.2306.

Whitaker, J. S., and T. M. Hamill, 2012: Evaluating methods to account for system errors in ensemble data assimilation. Mon. Wea. Rev., 140, 3078–3089, https://doi.org/10.1175/MWR-D-11-00276.1.

—, —, X. Wei, Y. Song, and Z. Toth, 2008: Ensemble data assimilation with the NCEP global forecast system. Mon. Wea. Rev., 136, 463–482, https://doi.org/10.1175/2007MWR2018.1.

Winiarek, V., M. Bocquet, O. Saunier, and A. Mathieu, 2012: Estimation of errors in the inverse modeling of accidental release of atmospheric pollutant: Application to the reconstruction of the cesium-137 and iodine-131 source terms from the Fukushima Daiichi power plant. J. Geophys. Res., 117, D05122, https://doi.org/10.1029/2011JD016932.

—, —, N. Duhanyan, Y. Roustan, O. Saunier, and A. Mathieu, 2014: Estimation of the caesium-137 source term from the Fukushima Daiichi nuclear power plant using a consistent joint assimilation of air concentration and deposition observations. Atmos. Environ., 82, 268–279, https://doi.org/10.1016/j.atmosenv.2013.10.017.

Wu, C. F. J., 1983: On the convergence properties of the EM algorithm. Ann. Stat., 11, 95–103, https://doi.org/10.1214/ aos/1176346060.

Yang, Y., and E. Mémin, 2019: Estimation of physical parameters under location uncertainty using an ensemble-expectation-maximization algorithm. Quart. J. Roy. Meteor. Soc., 145, 418–433, https://doi.org/10.1002/qj.3438.

Ying, Y., and F. Zhang, 2015: An adaptive covariance relaxation method for ensemble data assimilation. Quart. J. Roy. Meteor. Soc., 141, 2898–2906, https://doi.org/10.1002/qj.2576.

Zebiak, S. E., and M. A. Cane, 1987: A model El Niño–Southern Oscillation. Mon. Wea. Rev., 115, 2262–2278, https://doi.org/10.1175/1520-0493(1987)115<2262:AMENO>2.0.CO;2.

Zhang, F., C. Snyder, and J. Sun, 2004: Impacts of initial estimate and observation availability on convective-scale data assimilation with an ensemble Kalman filter. Mon. Wea. Rev., 132, 1238–1253, https://doi.org/10.1175/1520-0493(2004)132<1238:IOIEAO>2.0.CO;2.

Zhen, Y., and J. Harlim, 2015: Adaptive error covariances estimation methods for ensemble Kalman filters. J. Comput. Phys., 294, 619–638, https://doi.org/10.1016/j.jcp.2015.03.061.
## 4.3 Lguensat, Tandeo, Ailliot, Pulido et Fablet (2017) [MWR]

**Contexte** En 2014, je suis allé à Boulder (Colorado) présenter mes travaux à la conférence Climate Informatics. J'avais choisi cette conférence car elle cherche à fédérer les communautés climat et apprentissage automatique. Le travail que j'y ai présenté discutait de méthodes d'assimilation de données sans utiliser de modèle dynamique basé sur les équations de la physique. L'alternative que je proposais était d'utiliser des approches statistiques pour émuler le modèle dynamique. Ce travail reçut un très bon accueil et fut publié comme premier chapitre d'un livre de recueil d'articles liés à cette conférence (TANDEO, AILLIOT, RUIZ et al., 2015). Suite à ce travail préliminaire, j'ai encadré une partie du travail du doctorant Redouane Lguensat et nous avons publié un article plus abouti dans MWR.

**Résumé** Compte tenu de l'intérêt croissant pour les méthodes basées sur les données dans les sciences océaniques, atmosphériques et climatiques, ce travail présente l'assimilation de données par analogues (AnDA). Le cadre proposé fournit une reconstruction de la dynamique n'ayant aucune connaissance explicite du modèle dynamique. Au lieu de cela, un catalogue représentatif des trajectoires du système est supposé être disponible. Sur la base de ce catalogue, l'assimilation de données par analogues combine l'échantillonnage non paramétrique de la dynamique à l'aide de méthodes de prévision par analogues et de méthodes ensemblistes d'assimilation. Cette étude explore différentes stratégies de prévision par analogues, par filtrage de Kalman ou par filtrage particulaire. Des expériences numériques sont réalisées pour deux systèmes dynamiques chaotiques : les modèles de Lorenz-63 et Lorenz-96. L'assimilation de données par analogues est comparée à l'assimilation de données classique, dans laquelle les équations des modèles sont disponibles. Une boîte à outils Matlab et une bibliothèque Python d'AnDA sont fournies pour faciliter les recherches futures en s'appuyant sur les présents résultats.
# The Analog Data Assimilation™

REDOUANE LGUENSAT AND PIERRE TANDEO

IMT Atlantique, Lab-STICC, Université Bretagne Loire, Brest, France

PIERRE AILLIOT

Laboratoire de Mathématiques de Bretagne Atlantique, University of Western Brittany, Brest, France

MANUEL PULIDO

Department of Physics, Universidad Nacional del Nordeste, and CONICET, Corrientes, Argentina

RONAN FABLET

IMT Atlantique, Lab-STICC, Université Bretagne Loire, Brest, France

(Manuscript received 23 November 2016, in final form 31 July 2017)

## ABSTRACT

In light of growing interest in data-driven methods for oceanic, atmospheric, and climate sciences, this work focuses on the field of data assimilation and presents the analog data assimilation (AnDA). The proposed framework produces a reconstruction of the system dynamics in a fully data-driven manner where no explicit knowledge of the dynamical model is required. Instead, a representative catalog of trajectories of the system is assumed to be available. Based on this catalog, the analog data assimilation combines the nonparametric sampling of the dynamics using analog forecasting methods with ensemble-based assimilation techniques. This study explores different analog forecasting strategies and derives both ensemble Kalman and particle filtering versions of the proposed analog data assimilation approach. Numerical experiments are examined for two chaotic dynamical systems: the Lorenz-63 and Lorenz-96 systems. The performance of the analog data assimilation is discussed with respect to classical model-driven assimilation. A Matlab toolbox and Python library of the AnDA are provided to help further research building upon the present findings.

## 1. Introduction

The reconstruction of the spatiotemporal dynamics of geophysical systems from noisy and/or partial observations is a major issue in geosciences. Variational and stochastic data assimilation schemes are the two main categories of methods considered to address this issue [see Evensen (2007) for more details]. A key feature of these data assimilation schemes is that they rely on repeated forward integrations of an explicitly known dynamical model. This may greatly limit their application

range as well as their computational efficiency. First, thorough and time-consuming simulations may be required to identify explicit representations of the dynamics, especially regarding finescale effects and subgrid-scale processes as for instance in regional geophysical models (Hong and Dudhia 2012). Such processes typically involve highly nonlinear and local effects (Wilby and Wigley 1997). The resulting numerical models may be computationally intensive and even prohibitive for assimilation problems, for instance regarding the time integration of members with different initial conditions at each time step. Second, as explained in Van Leeuwen (2010), “with ever-increasing resolution and complexity, the numerical models tend to be highly nonlinear and also observations become more complicated and their relation to the models more nonlinear” (p. 1991). In such situations, standard data assimilation techniques may find difficulties, including

Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/MWR-D-16-0441.s1.

Corresponding author: Redouane Lguensat, redouane.lguensat@imt-atlantique.fr

DOI: 10.1175/MWR-D-16-0441.1

© 2017 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).
nonlinear particle filters which are prone to the “curse of dimensionality.” Third, difficulties may occur when geophysical dynamics involve uncertain model param- eterizations or space-time switching between different dynamical modes that need to be estimated online (Ruiz et al. 2013) or offline (Tandeo et al. 2015b). Dealing with such situations may not be straightforward using classi- cal model-driven assimilation schemes.

Meanwhile, recent years have witnessed a prolifera- tion of satellite data, in situ monitoring, as well as numerical simulations. Large databases of valuable information have been collected and offer a major opportunity for oceanic, atmospheric, and climate sci- ences. As pioneered by Lorenz (1969), the availability of such datasets advocates for the development of analog forecasting strategies, which make use of “similar” states of the dynamical system of interest to generate realistic forecasts. Analog forecasting strategies have become more and more popular in oceanic and atmo- spheric sciences (Nagarajan et al. 2015; McDermott and Wikle 2016), and have benefited from recent advances in machine learning (Zhao and Giannakis 2014). They have been applied to a variety of systems and applica- tion domains, including among others, rainfall now- casting (Atencia and Zawadzki 2015), air quality analysis (Delle Monache et al. 2014), wind field down- scaling (He-Guelton et al. 2015), climate reconstruction (Schenk and Zorita 2012), and stochastic weather gen- erators (Yiou 2014).

In this work, we examine the extension of the analog forecasting paradigm for data assimilation issues. Given a representative dataset of the dynamics of the system, this extension that we call analog data assimili- ation (AnDA) consists of a combination of the implicit analog forecasting of the dynamics with stochastic fil- tering schemes, namely, ensemble Kalman and particle filtering schemes (Evensen and Van Leeuwen 2000). This idea was first introduced in Tandeo et al. (2015a) where the relevance of the proposed analog data as- simulation is shown for the reconstruction of complex dynamics from partial and noisy observations. Tandeo et al. derived filtering and smoothing algorithms called the analog ensemble Kalman filter and smoother , which combine analog forecasting and the ensemble Kalman filter and smoother. A similar philosophy was followed independently in Hamilton et al. (2016) where the au- thors combine ideas from Takens’s embedding theorem and ensemble Kalman filtering to infer the hidden dy- namics from noisy observations. Hamilton et al. called their algorithm the Kalman–Takens filter .

Whereas these two previous works provide proofs of concept, our study further investigates and evaluates different analog assimilation strategies and their detailed

implementation. Our contributions are threefold. First, we present and examine various analog forecasting strategies, including locally linear ones that were not considered in previous works, and evaluate their performance for analog data assimilation. Second, in addition to the ensemble Kalman algorithms, we propose and examine a novel implementation of the analog forecasting combined with a particle filter. Finally, in the online supplemental material, we provide a unified computational framework, through both a Matlab Toolbox and a Python Library, to pave the way for practical use and future research (https://github. com/ptandeo/AnDA).

The work is organized as follows. In section 2, we briefly present the general concepts of data assimilation and introduce the key ideas of analog data assimilation. Different analog forecasting strategies are introduced in section 3. Section 4 describes the different components of the proposed analog data assimilation framework and the associated algorithms. Numerical experiments for two classical chaotic dynamical systems are reported in section 5. Section 6 further discusses our work, high- lights our key contributions, and proposes possible di- rections for future work.

2. General context

a. Model-driven data assimilation

Classically, data assimilation is based on the following discrete state space (Bocquet et al. 2010):

$$
\mathbf{x}(t) = \mathcal{M}[\mathbf{x}(t-1), \eta(t)], \quad (1)
$$

$$
\mathbf{y}(t) = \mathcal{H}[\mathbf{x}(t)] + \varepsilon(t), \quad (2)
$$

where time t ∈ { 0 , . . . , T } refers to the times in which observations are available. For the sake of simplicity we assume observations are at regular time steps.

In (1), M characterizes the dynamical model of the true state x ( t ) , while η ( t ) is a random perturbation added to represent model uncertainty. The observation equation (2) describes the relationship between the observation y ( t ) and x ( t ) . Observation error is consid- ered through the random noise ε ( t ) . Here, for the sake of simplicity, we consider an additive Gaussian noise ε with covariance R in (2) and the observation operator . H = H is assumed linear.

Data assimilation aims to reconstruct the state se- quence { x ( t ) } from a series of observations { y ( t ) } . Two types of data assimilation schemes are extensively studied in the literature: variational and stochastic. Variational data assimilation proceeds by minimizing a cost function based on a continuous formulation of (1) and (2) (see Lorenc et al. 2000), while stochastic data assimilation
<figure><img src="image_51.png" /><figcaption>FIG. 1. The evolution in time of one particle or member. The catalog implicitly represents the dynamics of the system from exemplars of historical datasets. The observations are shown by black asterisks, and their variance is shown by the corresponding error bar.</figcaption></figure>

schemes rely on the sampling and/or maximization of the posterior likelihood of the state sequence given the observation series (see Kalnay 2003). These classical data assimilation schemes are regarded as “model driven,” in the sense that they combine observations with forecasts provided by a numerical model $\mathcal{M}$.

### b. Data-driven data assimilation

The proposed assimilation framework relies on a similar state-space formulation. The key feature is to substitute the explicit dynamical model $\mathcal{M}$ in (1) by a “data driven” dynamical model involving an analog forecasting operator, denoted by $\mathcal{A}$, namely,

$$ \mathbf{x}(t) = \mathcal{A}[\mathbf{x}(t-1), \eta(t)]. \quad (3) $$

Henceforth, this state-space model will be referred to as AnDA. A sequential and stochastic data assimilation scheme including filtering and smoothing, is used involving different Monte Carlo realizations of the state at each assimilation time. We sketch the proposed AnDA methodology for one realization in Fig. 1.

The analog forecasting operator $\mathcal{A}$ requires the existence of a representative dataset of exemplars of the considered dynamics. This dataset is referred to as the catalog and denoted by $\mathcal{C}$. The reference catalog is formed by pairs of consecutive state vectors, separated by the same time lag. The second component of each pair is referred to as the successor of the first component hereafter. The catalog may be issued from observational data as well as from numerical simulations. In the last case, one can have a catalog issued from numerical simulations (based on physical equations), and wants to perform data assimilation without running the model

again. This is for instance useful for operational prediction centers that do not have the computational resources to integrate a forecast model, but do have access to a large database of numerical simulations or analysis data of a large prediction center. In this respect, we discuss also the situation where the catalog comprises noisy versions of the true states (section 5d).

Given a catalog $\mathcal{C}$, the analog forecasting operator $\mathcal{A}$ is stated as an exemplar-based statistical emulator of the state $\mathbf{x}$ from time $t$ to time $t+dt$. For any state $\mathbf{x}(t)$, we emulate the following state at time $t+dt$ based on its nearest neighbors in catalog $\mathcal{C}$. Given the analog forecasting operator, we present associated stochastic assimilation schemes, namely the analog ensemble Kalman filter/smoothen (Tandeo et al. 2015a) and the analog particle filter.

## 3. Analog forecasting strategies

### a. Analog forecasting operator

Let us consider a kernel function, denoted by $g$, in the state space (Schölkopf and Smola 2001). Among the classical choices for kernels, we consider here a radial basis function (also referred to as a Gaussian kernel):

$$ g(u, v) = \exp(-\lambda \|u - v\|^2), \quad (4) $$

where $\lambda$ is a scale parameter, $(u, v)$ are variables in the state space $\mathcal{B}$, and $\|\cdot\|$ is the Euclidean distance or another appropriate distance function. Note that the proposed analog forecasting operator may be applied to other kernels or subspace reduction methods to efficiently retrieve relevant analog situations. This is discussed in section 6.

Given the considered kernel, the analog forecasting operator $\mathcal{A}$ is defined as follows: for a given state $\mathbf{x}(t)$, we denote by $a_k[\mathbf{x}(t)]$ its $k$th nearest neighbor (or analog situation) in the reference catalog of exemplars $\mathcal{C}$, and by $s_k[\mathbf{x}(t)]$ the known successor of state $a_k[\mathbf{x}(t)]$. Hereinafter, we refer by $K$ to the number of nearest neighbors (analogs), and by $\text{cov}_w$ the weighted covariance. The normalized kernel weight for every pair $\{a_k[\mathbf{x}(t)], s_k[\mathbf{x}(t)]\}$ is given by

$$ \omega_k[\mathbf{x}(t)] = \frac{g\{\mathbf{x}(t), a_k[\mathbf{x}(t)]\}}{\sum_{K=1}^{K} g\{\mathbf{x}(t), a_k[\mathbf{x}(t)]\}}. \quad (5) $$

Several ideas can be explored to define the analog forecasting operator $\mathcal{A}$. The natural first option consists in deriving the forecast using the weighted mean of the $K$ successors. This approach, that we call here the locally constant operator, was considered in many analog forecasting related works (McDermott and Wikle 2016;
<figure><img src="image_52.png" /><figcaption>FIG. 2. A simplified illustration of the considered analog forecasting strategies in the case of two analogs (nearest neighbors). Two situations for the state $x(t)$ are shown: (top) a situation where $x(t)$ lies in the convex hull spanned by catalog exemplars and (bottom) a situation where $x(t)$ lies farther from its analogs. The second situation is expected to occur more often for high-dimensional space as well as for states, which are less likely. The latter may model extreme events or outliers.</figcaption></figure>

Zhao and Giannakis 2014; Hamilton et al. 2016), and is also known in statistics as Nadaraya-Watson kernel regression. One can also use as analog forecasting operator the weighted mean of the anomalies between the K analogs and their successors and adding it to the state to derive the forecast. The operator, referred to as locally incremental, is seen as more physically sound and relates more closely to a finite-difference approximation of the underlying differential equations. Finally, we introduce in this work a new analog forecasting operator that makes use of local linear regression techniques based on weighted least squares estimates. This operator that we call the locally linear operator is known to make an efficient use of small datasets and to reduce biases (Cleveland 1979). Note that the locally constant and locally incremental operators are two special cases of the locally linear operator.

Figure 2 shows an illustration of the three analog forecasting operators used in this work. Hereafter, we denote the forecasted state as $\mathbf{x}^f(t+dt)$. The three analog forecasting operators are defined as follows for two

sampling schemes: a Gaussian sampling and a multinomial one. Hereinafter, $\delta_Z(\cdot)$ denotes a delta function centered on Z.

* Locally constant analog operator: for the Gaussian case, the forecasted state is sampled from a Gaussian distribution whose mean $m_{LC}$ and covariance $\Sigma_{LC}$ are the weighted mean and the weighted covariance estimated from the K successors and their weights:

$$ \mathbf{x}^f(t+dt) \sim \mathcal{N}(m_{LC}, \Sigma_{LC}), \quad (6) $$

where $m_{LC} = \sum_{k=1}^K \omega_k[\mathbf{x}(t)]s_k[\mathbf{x}(t)]$ and $\Sigma_{LC} = \text{cov}_\omega(s_k[\mathbf{x}(t)])_{k \in \{1,K\}}$. While in the multinomial case, the forecasted state is drawn from the multinomial discrete distribution that samples the successor $s_k[\mathbf{x}(t)]$ with a probability of $\omega_k$:

$$ \mathbf{x}^f(t+dt) \sim \sum_{k=1}^{K} \omega_k[\mathbf{x}(t)] \delta_{s_k[\mathbf{x}(t)]}(\cdot). \quad (7) $$

* Locally incremental analog operator: instead of considering a weighted mean of the K successors as
in the locally constant operator, we consider the value of the current state plus a weighted mean of the K increments $\tau_k$, that is, the differences between analogs and successors $\tau_k[\mathbf{x}(t)] = s_k[\mathbf{x}(t)] - a_k[\mathbf{x}(t)]$. The Gaussian sampling is given by

$$ \mathbf{x}^f(t+dt) \sim \mathcal{N}(m_{\text{LI}}, \Sigma_{\text{LI}}), \quad (8) $$

where $m_{\text{LI}} = \mathbf{x}(t) + \sum_{k=1}^K \omega_k[\mathbf{x}(t)]\tau_k[\mathbf{x}(t)] = \sum_{k=1}^K \omega_k[\mathbf{x}(t)]\{\mathbf{x}(t) + \tau_k[\mathbf{x}(t)]\}$ and $\Sigma_{\text{LI}} = \text{cov}_\omega(\{\mathbf{x}(t) + \tau_k[\mathbf{x}(t)]\}_{k \in \{1, S\}})$ and the multinomial sampling resorts to

$$ \mathbf{x}^f(t+dt) \sim \sum_{k=1}^{K} \omega_k[\mathbf{x}(t)] \delta_{\mathbf{x}(t)+\tau_k[\mathbf{x}(t)]}(\cdot). \quad (9) $$

* Locally linear analog operator: we fit a multivariate linear regression between the K analogs of the current state and their corresponding successors using weighted least squares estimates (see Cleveland 1979). The regression gives slope $\alpha[\mathbf{x}(t)]$ and intercept $\beta[\mathbf{x}(t)]$ parameters, and residuals $\xi_k[\mathbf{x}(t)] = s_k[\mathbf{x}(t)] - (\alpha[\mathbf{x}(t)]a_k[\mathbf{x}(t)] + \beta[\mathbf{x}(t)])$. The Gaussian sampling comes to

$$ \mathbf{x}^f[(t+dt)] \sim \mathcal{N}(m_{\text{LL}}, \Sigma_{\text{LL}}), \quad (10) $$

with $m_{\text{LL}} = \alpha[\mathbf{x}(t)]\mathbf{x}(t) + \beta[\mathbf{x}(t)]$ and $\Sigma_{\text{LL}} = \text{cov}(\xi_k[\mathbf{x}(t)]_{k \in \{1, K\}})$, while the multinomial sampling is given by

$$ \mathbf{x}^f(t+dt) \sim \sum_{k=1}^{K} \omega_k[\mathbf{x}(t)] \delta_{m_{\text{LL}}+\xi_k[\mathbf{x}(t)]}(\cdot). \quad (11) $$

The choice of one operator over another depends mostly on the available computational resource and the complexity of the application. Locally constant and locally increment operators are less time and memory consuming than the locally linear operator, and while they can be of comparable performance in case of a flat regression function, the locally linear is expected to better deal with curvier regression functions at the expense, however, of the requirement of a larger number of analogs to fit the regression (Hansen 2000). The locally linear and the locally incremental are more suitable for samples near or outside the boundary of the select analogs (as depicted in Fig. 2), this may be particularly relevant in geoscience applications where chaos and extreme events are of high interest.

## b. Global and local analogs

The global analog strategy is the direct application of the introduced analog forecasting strategies to the entire state vector. We also introduce a local analog forecasting operator. For a given state $\mathbf{x}(t)$, the analogs $a_k[\mathbf{x}(t)]$ in the reference catalog, and their associated

successors $s_k[\mathbf{x}_l(t)]$ for each component $l$ of the state $\mathbf{x}(t)$ are defined according to a component-wise local neighborhood, typically $\{\mathbf{x}_{l-\nu}(t), \dots, \mathbf{x}_l(t), \dots, \mathbf{x}_{l+\nu}(t)\}$ with $\nu$ being the width of the considered component-wise neighborhood, such that the evaluation of the kernel function and the computation of the associated normalized weights $\omega_k[\mathbf{x}_l(t)]$ only involve this local neighborhood.

The idea of using local analogs is motivated by the fact that points tends to scatter far away from each other in high dimensions, which make the search for skillful analogs nearly impossible for high-dimensional state space. For instance, Van den Dool (1994) has shown that finding a relevant analog at synoptic scale over the Northern Hemisphere for atmospheric data would require 10³⁰ years of data to match the observational errors at that time. Conversely, analog forecasting schemes may only apply to systems or subsystems associated with low-dimensional embedding. Following this analysis, the analog forecasting of the global state is split as a series of local and low-dimensional analog forecasting operations. Note that such local analogs also reduce possible spurious correlations.

## 4. Analog data assimilation

The analog data assimilation is stated as a sequential and stochastic assimilation scheme, using Monte Carlo methods. It amounts to estimating the so-called filtering and smoothing posterior likelihoods, respectively, $p[\mathbf{x}(t)|\mathbf{y}(1), \dots, \mathbf{y}(T)]$ the distribution of the current state knowing past and current observations and $p[\mathbf{x}(t)|\mathbf{y}(1), \dots, \mathbf{y}(T)]$ the distribution of the current state knowing past, current, and future observations. We investigate both ensemble Kalman filter/smoothen and particle filter.

### a. Analog ensemble Kalman filter and smoother (AnEnKF/AnEnKS)

Ensemble Kalman filters (EnKF) and smoothers (EnKS) (Burgers et al. 1998; Evensen 2007) are particularly popular in geoscience as they provide flexible assimilation strategies for high-dimensional states. They rely on the assumption that the filtering and smoothing posteriors are multivariate Gaussian distributions, such that the following forward and backward recursions are derived. The next two paragraphs present the AnEnKF and AnEnKS equations, which are equivalent to those of the EnKF and EnKS described in Tandeo et al. (2015b), except for the update step where we use the analog forecasting operator.

The forward recursions of the AnEnKF correspond to the stochastic EnKF algorithm proposed by Burgers
et al. (1998) in which observations are treated as random variables. The AnEnKF algorithm starts at time $t = 1$ by generating the vectors $\mathbf{x}_i^f(1) \forall i \in \{1, ..., N\}$ using a multivariate Gaussian random generator with mean vector $\mathbf{x}^b$ and covariance matrix $\mathbf{B}$. The index $i$ of the state vector corresponds to the $i$th realization of the Monte Carlo procedure (called member or particle). Then the update step proceeds from $t = 2$ to $t = T$ by applying the analog forecasting operator to each member of the ensemble following (3) to generate $\mathbf{x}_i^f(t)$. The forecast state is represented by the sample mean $\mathbf{x}^f(t)$ and the sample covariance $\mathbf{P}^f(t)$. In the analysis step, following (2), $N$ samples of $\mathbf{y}_i^f(t)$ are generated from a multivariate Gaussian random generator with mean $\mathbf{H}\mathbf{x}_i^f(t)$ and covariance $\mathbf{R}$. The observations are then used to update the $N$ members of the ensemble as $\mathbf{x}_i^a(t) = \mathbf{x}_i^f(t) + \mathbf{K}^a(t)[\mathbf{y}(t) - \mathbf{y}_i^f(t)]$, where $\mathbf{K}^a(t) = \mathbf{P}^f(t)\mathbf{H}'[\mathbf{H}\mathbf{P}^f(t)\mathbf{H}' + \mathbf{R}]^{-1}$ is the Kalman filter gain. The filtering posterior distribution is then represented by the sample mean $\mathbf{x}^a(t)$ and the sample covariance $\mathbf{P}^a(t)$.

The analog ensemble Kalman smoother combines the analog forecasting operator and the classical Kalman smoother, here, Rauch-Tung-Striebel smoother [see Cosme et al. (2012) for more details]. Given the forward recursion, the backward recursion starts from time $t = T$ with filtered state, $\forall i \in \{1, ..., N\}$, such as $\mathbf{x}_i^s(T) = \mathbf{x}_i^a(T)$ and $\mathbf{P}^s(T) = \mathbf{P}^a(T)$. Then, we proceed backward from $t = T - 1$ to $t = 1$. At each time $t$, we compute $\mathbf{x}_i^s(t) = \mathbf{x}_i^a(t) + \mathbf{K}^s(t)[\mathbf{x}_i^s(t + 1) - \mathbf{x}_i^f(t + 1)]$, where $\mathbf{K}^s(t) = \mathbf{P}^s(t)\mathcal{M}'[\mathbf{P}^f(t + 1)]^{-1}$ is the Kalman smoother gain. Note that we empirically estimate $\mathbf{P}^a(t)\mathcal{M}'$ as the sample covariance matrix of the ensemble members as in Pham (2001) or Tandeo et al. (2015b) in the case of a nonlinear operator $\mathcal{M}$. The smoothing posterior distribution is represented by the sample mean $\mathbf{x}^s(t)$ and the sample covariance $\mathbf{P}^s(t)$.

We note that the following way of extending EnKF and EnKS to become analog-based algorithms can be applied in the same way to other flavors of EnKF such as the square root ensemble Kalman filter (EnSRF). We chose stochastic ensemble-based Kalman filters and smoothers as an illustration in this work, even if they are not the first choice in practice for atmospheric and oceanic applications because of issues related to perturbing observations with noise (Bowler et al. 2013). Besides, the work of Hoteit et al. (2015), where the authors address this issue, suggests that the stochastic EnKF is worth a reevaluation for oceanic and atmospheric applications.

## b. Analog particle filter (AnPF)

We also implement particle filtering techniques for the proposed analog data assimilation strategy. Contrary to

the Kalman filters, particle filters do not assume a Gaussian distribution of the state. The key principle is to estimate the posteriors of the state from a set of particles (equivalent to members in the terminology used for ensemble Kalman filters).

Given an analog forecasting operator, we consider an application of the classical particle filter (Van Leeuwen 2009). From an initialization similar to the EnKF, the particle filter applies a forward recursion from time $t = 1$ to $t = T$ as follows. At time step $t$, we first apply the considered analog forecasting operator $\mathcal{A}$ to forecast $\mathbf{x}_i^f(t) \forall i \in \{1, ..., N\}$ from previous filtered particles $\mathbf{x}_i^a(t - 1)$. Then, following (2), we compute particle weights $\pi_i(t)$ as

$$ \pi_i(t) \propto \phi[\mathbf{y}(t) - \mathbf{H}\mathbf{x}_i^f(t); \mathbf{R}], \quad (12) $$

where $\phi(\cdot; \mathbf{R})$ is a centered multivariate Gaussian distribution with covariance $\mathbf{R}$. Weights $\pi_i(t)$ are normalized to total one. We then proceed to a systematic resampling from the multinomial distribution defined by the particles $\{\mathbf{x}_i^f(t)\}$ and their corresponding weights $\{\pi_i(t)\}$. The analyzed state $\mathbf{x}^a(t)$ is typically computed as the sample mean

$$ \mathbf{x}^a(t) = \frac{1}{N} \sum_{i=1}^{N} \pi_i(t) \mathbf{x}_i^f(t), \quad (13) $$

but one may also consider the posterior mode as the filtered state.

In theory, particle smoothers may also be considered. Different strategies have been proposed in the past but they showed numerical instabilities in preliminary experiments with the considered analog forecasting operator. We do not further detail the considered implementation but discuss these aspects in section 6.

# 5. Numerical experiments

To evaluate the relevance and performance of the proposed analog data assimilation, we consider numerical experiments on dynamical systems extensively used in the literature on data assimilation: Lorenz-63 and Lorenz-96 models. The experiments for evaluating the effect of the size of the catalog, the impact of noisy catalogs, and catalogs with parametric model error are conducted using the Lorenz-63 model. To evaluate the global and local analog forecasting operators we use the Lorenz-96 model, an extended dynamical nonlinear system with 40 variables.

## a. Chaotic models

We first consider the chaotic Lorenz-63 system. From a methodological point of view, it is particularly
interesting because of its nonlinear chaotic behavior and low dimension. Several works have used this system (e.g., Miller et al. 1994; Anderson and Anderson 1999; Pham 2001; Chin et al. 2007; Hoteit et al. 2008 or Van Leeuwen 2010). The Lorenz-63 model is defined by

$$ 
\begin{aligned}
\frac{dx_1(t)}{dt} &= \sigma[x_2(t) - x_1(t)], \\
\frac{dx_2(t)}{dt} &= x_1(t)[\gamma - x_3(t)] - x_2(t), \\
\frac{dx_3(t)}{dt} &= x_1(t)x_2(t) - \beta x_3(t),
\end{aligned}
\quad (14) $$

and behaves chaotically for certain sets of parameters, such as ($\sigma = 10$, $\gamma = 28$, $\beta = 8/3$). Here, we use the explicit (4, 5) Runge-Kutta integrating method (cf. Dormand and Prince (1980)) with time step $dt = 0.01$ (nondimensional units). As in Van Leeuwen (2010) only the first variable of the Lorenz-63 system ($x_1$) is observed every 8 integration time steps (i.e., with $dt = 0.08$). Considering the analogy between the Lorenz-63 and atmospheric time scales, it is equivalent to a 6-h time step in the atmosphere.

The Lorenz-96 model is another chaotic model largely used for evaluating data assimilation techniques in geophysics (Anderson 2001; Whitaker and Hamill 2002; Ott et al. 2004; Anderson 2007, 2012; Hoteit et al. 2012). It is defined by

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

where $j = 1, \dots, n$ and the boundaries are cyclic [i.e., $x_{-1}(t) = x_{n-1}(t)$, $x_0(t) = x_n(t)$, and $x_{n+1}(t) = x_1(t)$]. The three right-hand side terms in (15) simulate an advection, a diffusion, and a forcing term, respectively. As in Lorenz (1996), we choose $n = 40$ and external forcing of $F = 8$ for which the model behaves chaotically. Equation (15) is solved using the Runge-Kutta fourth-order scheme with integration time step $dt = 0.05$, corresponding to a time step of 6 h in the atmosphere. Observations are taken from half of the state vector (20 observed components randomly selected) every 4 time steps (i.e., $dt = 0.20$).

## b. Experimental details

The considered experimental setting is as follows. To avoid divergence of the filtering methods, we use $N = 100$ members/particles for the Lorenz-63 and $N = 1000$ members/particles for the Lorenz-96 for both model-driven and data-driven strategies. We use the same covariance matrix $R$ with a noise observation variance set to 2. To avoid any spinup effect, the initial state conditions is chosen as the ground truth mean and a covariance matrix

$B$ with noise variance 0.1. To compare the technique performances, we use the root-mean-square error (RMSE) on all the components of the state vector and for all assimilation times. As training dataset for the catalog and test dataset for RMSE computation, we use $10^3$ and 100 Lorenz times, respectively.

The analog forecasting operator involves two free parameters, namely, $K$ the number of nearest neighbors and $\lambda$ the scale parameter of the Gaussian kernel in (4). Two strategies can be considered for $K$: either a pre-defined number of nearest neighbors, or a predefined threshold on distance $d_{th}$ to select the analogs that are closer than $d_{th}$. For the sake of simplicity, we consider in this work the first alternative and set $K$ to 50. Besides, we use for $\lambda$ the following adaptive rule: $\lambda[x(t)] = 1/\text{md}[x(t)]$, where $\text{md}[x(t)]$ is the median distance between the current state $x(t)$ and its $K$ analogs. Note that a cross-validation procedure could be used to optimize the choice of $K$ and $\lambda$. All analog forecasting operators are fitted for forecasting time horizon corresponding to the time step of the numerical simulations (i.e., $dt = 0.01$ for Lorenz-63 experiments and $dt = 0.05$ for Lorenz-96 experiments). Numerical experiments (not reported here) show that this parameterization provides on average the best forecasting performance with respect to the forecasting time horizon.

## c. Experiments with Lorenz-96 model

### 1) EXPERIMENT 1

The first numerical experiment consisted only in the application of analog forecasting (without assimilation) from a catalog. We build a database using Lorenz-96 equations, then we split the samples randomly to 2/3 for training the analog forecasting operators and 1/3 for test. Finally, we compare the RMSE w.r.t. ground truth data as a function of Lorenz-96 forecast time. For local analogs, we consider $\nu = 2$ the width of the considered component-wise neighborhood. Figure 3 shows the results of this experiment using the three choices for the analog forecasting operator $\mathcal{A}$. The locally linear approach outperforms the two other approaches confirming that its forecasts are with lower bias compared to the other approaches. However, it also involves more parameters, which increases the variance of the forecasts. This bias-variance trade-off supports the greater generalization capabilities of the locally linear operator, when the dynamics can well be approximated locally by a linear operator.

Figure 3 also compares local and global analog strategies. When using locally constant operator, local analogs are always better than global analogs. Searching for nearest neighbors on 40-dimensional vectors results
<figure><img src="image_53.png" /><figcaption>FIG. 3. Results of the analog forecasting performance as a function of the horizon. Different analog forecasting methods are plotted: locally constant (green), locally incremental (blue), and locally linear (red) analog operators with local (straight line) and global (dashed line) analog strategies. The black dashed line corresponds to a persistent prediction over time.</figcaption></figure>

most likely in irrelevant analogs. This affects heavily the locally constant operator more than the two other operators, since it computes a weighted mean of their associated successors. The locally constant operator also limits novelty creation in the dynamics by always dragging the forecast near the mean of the *K* successors, and, according to these experiments, it seems poorly adapted to complex and highly nonlinear systems. Regarding the locally incremental and locally linear strategies, local analogs are more relevant than global ones for prediction in a near future (less than 0.5 in Lorenz-96 time for locally linear operator and less than 0.25 in Lorenz-96 time for locally incremental).

## 2) EXPERIMENT 2

We conducted a second experiment for evaluating the impact of analog forecasting in data assimilation using the Lorenz-96 model. We run the AnEnKS with 1000 ensemble members, only 20 variables are observed every 0.20 time steps. Figure 4 shows analog data assimilation experiments with the locally linear forecasting method using the Lorenz-96 model. Figures 4a and 4b show the true state and the observations, respectively. The reconstructed state with global analogs is shown in Fig. 4c and the one with local analogs in Fig. 4d. The local analog data assimilation experiment clearly outperforms the global analog data assimilation experiment.

## 3) EXPERIMENT 3

A third experiment with the Lorenz-96 system was conducted. For the local analog strategy, we further

compare the proposed AnDA algorithms, namely, AnEnKF, AnPF, and the AnEnKS using 1000 ensemble members/particles, in Table 1. Two main conclusions can be drawn: (i) EnKF algorithms outperform the particle filter and (ii) the locally linear analog forecasting operator gives the best reconstruction performance. We noticed that the AnPF suffers in the 40-dimensional Lorenz-96 system from sample impoverishment and degeneracy. Despite additional experiments with different settings, for instance, w.r.t. the number of ensemble members, the number of analogs as well as using jittering (i.e., perturbing the particles with a small noise), the AnPF still suffered from the aforementioned issues.

### d. Experiments with Lorenz-63 model

#### 1) EXPERIMENT 1

In the proposed AnDA, the size of the catalog is expected to be a critical parameter. For Lorenz-63 dynamics, we conducted different AnDA experiments varying the size of the catalog $S = \{10^1, 10^2, 10^3, 10^4\}$ in Lorenz-63 times. We consider the same setting as in Tandeo et al. (2015a) where the locally constant method with a Gaussian sampling was used for the AnEnKF, then we compare the three AnDA algorithms using 100 ensemble members/particles. As reported in Fig. 5, the RMSE decreases when the size of the catalog increases for all AnDA algorithms. Regarding filtering-only (i.e., no smoothing) AnDA algorithms, the AnPF (blue) outperforms the AnEnKF (green). This is an expected result since particle filters handle better nonlinear models and non-Gaussian probability distributions, although at a high cost in terms of computational complexity and execution time. The AnEnKS (red) clearly gives the lowest RMSE. This supports the additional benefit of the smoothing step performed by the AnEnKS. The zoom shown in the right panel of Fig. 5 highlights how the smoothing step corrects the piecewise effects resulting from the filtering step.

#### 2) EXPERIMENT 2

Modeling uncertainty is a critical source of error in data assimilation. In this experiment we evaluate whether AnDA can manage a situation in which the catalog is composed by multiple numerical simulations, which may have parametric model error. In (14), parameters $\gamma$ and $\beta$ define the center of the two attractors whereas $\sigma$ controls the shape of the trajectories. In Fig. 6, we depict trajectories using three sets of parameters with different values for $\sigma$: $\theta_1 = (10, 28, 8/3)$ (red), $\theta_2 = (7, 28, 8/3)$ (blue), and $\theta_3 = (13, 28, 8/3)$ (green). We generate three catalogs with Lorenz-63 trajectories