Pierre Tandeo
email: pierre.tandeo@imt-atlantique.fr

Pierre Ailliot

Marc Bocquet

Alberto Carrassi

Takemasa Miyoshi

Manuel Pulido

Yicun Zhen

Joint Estimation of Model and Observation Error Covariance Matrices in Data Assimilation: a Review

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Data Assimilation (hereinafter denoted DA) for geosciences is generally formulated in terms of nonlinear state-space models with additive and Gaussian errors for the dynamical and observation equations. This is statistically convenient and representative of a lot of DA prob-lems, see e.g. [START_REF] Carrassi | Data Assimilation in the Geosciences: An overview on methods, issues and perspectives[END_REF]. The errors on the dynamics and observations are assumed to be zero-mean Gaussian vectors with covariance matrices Q and R. Using the discrete time index k from 1 to K for the sake of simplicity, it is assumed that

x(k) = M (k -1, x(k -1)) + η(k), (1) y
(k) = H (k, x(k)) + ε(k), (2)
with M the time dependent dynamical model, H the time dependent transformation operator from the hidden state Generated using v4.3.2 of the AMS L A T E X template 1

x to the noisy observations y, Gaussian white noise errors η(k) ∼ N (0, Q(k)) and ε(k) ∼ N (0, R(k)). We suppose that the initial state condition at k = 0 is a Gaussian vector with mean x b and variance B and that η and ε are mutually independent. However, in some situations it may be relevant to consider cross-correlation between these errors (see [START_REF] Berry | Correlation between system and observation errors in data assimilation[END_REF] for more details). DA algorithms are used to estimate sequentially the state of the system x conditionally to the observations y. When current and past observations are used, the estimation is referred as filtering or analysis and when future observations are also used, it is referred as smoothing or reanalysis. The outcome of the analysis/filter or of the reanalysis/smoother highly depends on the uncertainty associated to observations and to the model state, which have to be as realistic as possible. Using the formulation given in Eqs. [START_REF] Blanchet | A comparison of adaptive kalman filters for a tropical pacific ocean model[END_REF](2) these uncertainties are Q and R. In practice, the observation error covariance matrix R in Eq. ( 2) can be determined empirically by estimating the instrument noise and the representativeness error between the state and the observation space, but a correct estimation of the latter is often challenging (see [START_REF] Janjić | On the representation error in data assimilation[END_REF] for more details). During the dynamical model evolution from k -1 to k in Eq. ( 1), the model state is contaminated by two sources of uncertainty, the error in the state at k -1 and the model uncertainty itself which is represented in state-space models via the additive model error term η. Determining the covariance matrix Q model error is difficult because it accounts for the model deficiencies to represent the underlying physics, the cumulative effects of errors in the parameters, the numerical schemes, the unresolved scales and the fact that in geosciences, we usually have much less observations that those needed to estimate the entries of Q (see e.g. [START_REF] Daley | Estimating Model-Error Covariances for Application to Atmospheric Data Assimilation[END_REF] and [START_REF] Dee | On-line Estimation of Error Covariance Parameters for Atmospheric Data Assimilation[END_REF]). When using either variational or ensemble-based DA methods, the quality of the reconstructed state vector highly depends on the relative amplitudes between the assumed observation and model errors. For instance, in Kalman filter-like, the ratio ||Q||/||R|| impacts the filter gain that gives the relative weights of the observations against the model forecasts. [START_REF] Desroziers | Diagnosis and adaptive tuning of observation-error parameters in a variational assimilation[END_REF] also studied this ratio in variational DA. Unfortunately, in real DA frameworks, the impact of Q, R and ||Q||/||R|| on the reconstruction of the state is not easy to evaluate. This is due to the complexity and size of the dynamical models, the effect of forcing terms and the huge variety of observations.

The importance of estimating error covariance matrices in state-space models can be illustrated using a simple example with linear dynamics. Suppose that we aim at tracking a scalar state x governed by an autoregressive AR [START_REF] Blanchet | A comparison of adaptive kalman filters for a tropical pacific ocean model[END_REF] model in Eq. ( 1) defined by

x(k) = 0.95x(k -1) + η(k), (3)
with η ∼ N (0, Q t ) where the superscript t means "true" and Q t = 1. Furthermore, observations y of the state are contaminated with another independent additive zeromean and unit-variance Gaussian noise (i.e. R t = 1) in Eq. ( 2) with H (x) = x. The goal is to reconstruct x from the noisy observations y at each time step. The AR [START_REF] Blanchet | A comparison of adaptive kalman filters for a tropical pacific ocean model[END_REF] model defined by Eq. ( 3) has an autoregressive coefficient close to one and thus represents a process which evolves slowly in time. The linear dynamical model evolves stochastically and the measurement process also introduces a noise at each time step. Although, the knowledge of these two sources of noise is crucial for the estimation problem, in practice identifying them is not an easy task. Given that the dynamical model is linear and the error terms are additive and Gaussian in this simple example, the Kalman smoother provides an exact algorithm to compute the smoothing distribution (see Sect. 2 for more details). To evaluate the impact of badly specified Q and R errors on the reconstructed state with the Kalman smoother, different experiments were conducted using values of {0.1, 1, 10} for the ratio Q/R. Figure 1 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 line). We also report the Root Mean Squared Error (RMSE) of the reconstruction as well as the so called "coverage probability" or percentage of x falling in the 95% confidence intervals (defined as the mean ±1.96 the standard deviation in the Gaussian case). In this synthetic experiment, the best RMSE and coverage probability obtained using the Kalman smoother with true Q t = R t = 1 are respectively 0.71 and 95%. When using a low model error variance Q = 0.1Q t in Fig. 1(a), it gives an important weight to the forecasts given by quasi-persistent autoregressive dynamical model. On the other hand, when using a low observation error variance R = 0.1R t in Fig. 1(b), too much weight is given to the observation and the reconstructed state is too 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. 1(c), the mean reconstructed state vector and so its RMSE are similar to Fig. 1(b). In the same way, overestimated observation error variance as in Fig. 1(d) gives similar mean reconstruction than Fig. 1(a). These two last 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 Fig. 1(e) and Fig. 1(f) the case where Q/R ratio is good (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 the one 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 FIG. 1. Example of a univariate AR(1) process generated using Eq. ( 3) with Q t = 1 (red line), noisy observations using R t = 1 (black dots) and reconstructions using a Kalman smoother (black lines and gray 95% confidence interval) with different values of Q and R.

in Fig. 1(e) and overestimated in Fig. 1(f) with respective coverage probability of 36% and 100%.

Finally, let us consider again the Kalman smoother where one of the variances Q or R is erroneous. This time, we compensate the error in the wrong prescribed variance by optimizing the other free variance using the maximum likelihood estimation method given in [START_REF] Shumway | An Approach to Time Series Smoothing and Forecasting using the EM Algorithm[END_REF], see Sect. 4. Results presented in Fig. 2 show that best optimal RMSE (0.71) and coverage probability (95%) are reached close to the optimal variance noises Q t = R t = 1. Results also indicate that a compensation of bad variances is possible but is not optimal. For instance, when fixing Q to a bad value like 0.25 in Fig. 2(a) and Fig. 2(b), the maximum likelihood estimator of R is 1.57 and corresponding RMSE and coverage probability are respectively 0.86 and 80% (out of range of the color bar). These two skill metrics are extremely important to evaluate the quality of the reconstructed state. Nevertheless, often only the RMSE is presented in research papers. The coverage probability is a measure of the DA capability to quantify the uncertainty, a problem that we believe to be of increasing relevance for the DA community in the coming years. Indeed, the reconstructed state error variance may vary significantly as shown in Fig. 1(e) and Fig. 1(f). Here, we use a simple synthetic example but for large dimensional and highly nonlinear dynamics as in the DA in geosciences, such underestimation or overestimation of uncertainty may have a strong impact and may make filters to collapse. In this linear and Gaussian example, the use of RMSE and the probability of coverage are sufficient, but for nonlinear and more realistic DA cases, we should also consider the rank histograms and the proper scores.

Since the 90's, a significant number of works have dealt with the error covariances in state-space models. The first ones who mentioned the importance of noise covariance matrices Q and R in DA were [START_REF] Ghil | Data assimilation in meteorology and oceanography[END_REF] in their Sect. 4.1, as well as [START_REF] Daley | Atmospheric data analysis[END_REF] in his Sect. 4.9. [START_REF] Daley | Estimating Model-Error Covariances for Application to Atmospheric Data Assimilation[END_REF] clarified the difference between "predictability error" and "model error", the two components of the forecast error, denoted as P f in modern DA. As illustrated in Fig. 3, the first error is due to imperfect initial conditions and the second one is caused by model imperfections represented by Q. [START_REF] Dee | On-line Estimation of Error Covariance Parameters for Atmospheric Data Assimilation[END_REF] proposed a maximum likelihood estimator for parameterized versions of Q and R using the innovation likelihood criterion. Dee et al. (1999a) extended this online method to the estimation of the mean of the innovations, which depends on the biases in the forecast and in the observations, and later applied to realistic cases in Dee et al. (1999b). These initial studies clearly impulsed the treatment of this topic in modern DA literature and several works have appeared thereafter on the joint estimation of model and observation errors. However, authors like [START_REF] Todling | A lag-1 smoother approach to system-error estimation: sequential method[END_REF] pointed out that using only the current innovation is not enough to dis-FIG. 2. Estimation results of R and Q for respectively fixed Q and R in the case of a Kalman smoother with the univariate AR(1) process given in Eq. ( 3). Here, we consider the maximum likelihood estimator given by [START_REF] Shumway | An Approach to Time Series Smoothing and Forecasting using the EM Algorithm[END_REF].

tinguish the impact of Q and R in the Kalman equations, and to estimate them independently is challenging in this case. Thus, they proposed various alternatives to tackle this issue.

An history of what have been, in our opinion, the most relevant contributions and the key milestones for covariance estimation in geophysical systems is sketched in Fig. 4 and is discussed in this review with a summary given in Table 1. We distinguish four methodologies and among them, one could classify the approaches whether they rely upon the innovations or the likelihood. The innovations are defined as the difference between the observations and state estimates transformed to the observational space, with both the forecast and the analysis. The use of their corresponding statistics in the observation space has been initiated by [START_REF] Desroziers | Diagnosis of observation, background and analysis-error statistics in observation space[END_REF]. Then, this approach has been used extensively for the calibration of inflation of the forecast covariance with various implementations including additive, relaxation-to-prior and the multiplicative inflation case, see Li et al. (2009a) and [START_REF] Miyoshi | The Gaussian Approach to Adaptive Covariance Inflation and Its Implementation with the Local Ensemble Transform Kalman Filter[END_REF]. Instead of working on different innovations at a given time, [START_REF] Berry | Adaptive ensemble Kalman filtering of non-linear systems[END_REF] as well as [START_REF] Harlim | An ensemble Kalman filter for statistical estimation of physics constrained nonlinear regression models[END_REF] suggested to use lag-innovations or innovation between consecutive times. At the same time, methods based on likelihood functions and their maximization using statistical approaches appeared. Bayesian inference techniques with the use of prior distributions and hyperparameters as in [START_REF] Stroud | Sequential State and Variance Estimation within the Ensemble Kalman Filter[END_REF] or [START_REF] Stroud | A Bayesian adaptive ensemble Kalman filter for sequential state and parameter estimation[END_REF] are typical examples. Finally, [START_REF] Ueno | Iterative algorithm for maximumlikelihood estimation of the observation-error covariance matrix for ensemble-based filters[END_REF], [START_REF] Dreano | Estimating model error covariances in nonlinear state-space models using Kalman smoothing and the expectationmaximisation algorithm[END_REF] and [START_REF] Pulido | Stochastic parameterization identification using ensemble Kalman filtering combined with maximum likelihood methods[END_REF] proposed to maximize the total likelihood of the state-space model using iterative expectation-maximization algorithms.

The four methods mentioned above are detailed in this review and are factually described as they appear in the literature. We consider both online and offline estimations, for which the computational cost highly varies. In the online or adaptive approaches, we try to estimate a time-dependent Q(k) and R(k) at the same time as the state vector, using filtering methods. When considering offline or batch approaches, averaged Q and R are estimated using all the observations on a given time interval, using smoothing methods. Moreover, offline procedures are iterative, meaning that the procedures are repeated until convergence according to a given criterion, for instance the likelihood. Finally, in some methods presented here, additional tuning parameters are needed and have to be carefully chosen for practical implementations. We discuss this point in this review.

Note that other review papers on parameters estimation in state-space models appeared in the statistical and signal processing communities by [START_REF] Mehra | Approaches to Adaptive Filtering[END_REF], [START_REF] Kantas | On Particle Methods for Parameter Estimation in State-Space Models[END_REF] and [START_REF] Duník | Noise covariance matrices in state-space models: A survey and comparison of estimation methods-Part I[END_REF]. The [START_REF] Mehra | Approaches to Adaptive Filtering[END_REF] paper is a concise review which accounts for linear dynamical models using the classic Kalman filter while Kan-FIG. 3. Sketch of sequential data assimilation algorithms in the observation space where H is omitted for simplicity reasons. The ellipses represent the error covariances computed from the Kalman-based equations, i.e. the forecast P f and analysis P a , as well as the model error Q and observation error R, unknown entries of the state-space model in Eqs. [START_REF] Blanchet | A comparison of adaptive kalman filters for a tropical pacific ocean model[END_REF](2). This is a modified figure based on Fig. 1 from [START_REF] Carrassi | Data Assimilation in the Geosciences: An overview on methods, issues and perspectives[END_REF]. tas et al. (2015) and [START_REF] Duník | Noise covariance matrices in state-space models: A survey and comparison of estimation methods-Part I[END_REF] focus on the challenging covariance estimation in non-Gaussian Monte Carlo methods, and include various implementations of offline and online methods based on maximum likelihood and Bayesian approaches, but not methods based on innovation statistics. They evaluate the methods on onedimensional and easy-to-simulate nonlinear models, using a large number of members or particles. The [START_REF] Kantas | On Particle Methods for Parameter Estimation in State-Space Models[END_REF] and [START_REF] Duník | Noise covariance matrices in state-space models: A survey and comparison of estimation methods-Part I[END_REF] review papers address more theoretical aspects and asymptotic results. The current review is addressed to the DA community, so that we focus on methods which can be implemented for nonlinear operators and high-dimensional systems, in particular using Kalman formulation.

The review is organized as follows. Section 2 presents the filtering and smoothing DA algorithms used in this present work. Then, we present the main families of methods used in the literature to jointly estimate error covariance matrices Q and R. First, innovation-based methods are presented in Sect. 3. Then, we describe in Sect. 4 the likelihood-based approaches. We also list other alternatives in Sect. 5 with methods used in the past and not exactly matching the scope of this review, or diagnostic tools to check the accuracy of Q and R. Finally, in Sect. 6, we summarize this review and discuss possible perspectives for future works and remaining challenges to be faced in this domain.

Filtering and smoothing algorithms

For the overall discussion of the methods and for introduction of the notation, we present in this section a short description of the extended version of the Kalman equations for nonlinear dynamical systems and observation operators. Here we use time dependent linearizations M and H of the nonlinear operators M and H defined in Eqs. [START_REF] Blanchet | A comparison of adaptive kalman filters for a tropical pacific ocean model[END_REF](2). We have chosen to based the discussion on the Extended Kalman Filter and Smoother (EKF/EKS) in this review, compared to Ensemble Kalman Filter (EnKF) for instance, to avoid the overburdening notations introduced by the ensemble members. However, the methods are also straightforward to apply in stochastic and squareroot EnKFs.

Note that the most natural algorithms to solve the statespace model given in Eqs. (1-2) are the Particle Filter and Smoother (PF and PS) from [START_REF] Gordon | Novel approach to nonlinear/non-Gaussian Bayesian state estimation[END_REF] and firstly reviewed in DA by van [START_REF] Van Leeuwen | Particle filtering in geophysical systems[END_REF]. These methods converge to the true posterior distributions for a large number of particles in theory. However, we focus in this review on Gaussian additive errors η and ε in Eqs. (1-2), and EKF/EKS perform generally well in this situation. Moreover, the current PF and PS implementations are subject to the curse of dimensionality [START_REF] Snyder | Obstacles to High-Dimensional Particle Filtering[END_REF]) and are not suitable for high dimensional systems, although recent implementations appear to shed some light on these contention points (e.g., [START_REF] Atkins | Implicit Particle Methods and Their Connection with Variational Data Assimilation[END_REF] and [START_REF] Zhu | Implicit equalweights particle filter[END_REF]). However, because of PF relies on a state-space model formulation as in Eqs. (1-2), the model error covariance specification is an essential requirement in the definition of the transition density. Therefore, the estimation of Q with Gaussian Kalman-based methods may give a useful constraint or parameterization setup of the model error covariance matrix for PF (see [START_REF] Zhu | Estimating model error covariances using particle filters[END_REF]).

Kalman-based algorithms assume a Gaussian prior distribution p (x(k)|y( 1

: k -1)) ∼ N x f (k), P f (k) .
Then, filtering and smoothing estimates are corresponding to the Gaussian posterior distributions Forecast step (forward in time):

p (x(k)|y(1 : k)) ∼ N (x a (k), P a (k)) and p (x(k)|y(1 : K)) ∼ N (x s (k), P s (k)) of
x f (k) =M (x a (k -1)) (4) P f (k) =M(k)P a (k -1)M(k) + Q(k) (5)
Analysis step (forward in time):

d(k) =y(k) -H x f (k) (6) K f (k) =P f (k)H(k) H(k)P f (k)H(k) + R(k) -1 (7) x a (k) =x f (k) + K f (k)d(k) (8) P a (k) = I -K f (k)H(k) P f (k) (9)
Reanalysis step (backward in time):

K s (k) =P a (k)M(k) P f (k + 1)
-1

(10)

x s (k) =x a (k) + K s (k) x s (k + 1) -x f (k + 1) (11) P s (k) =P a (k) -K s (k) P f (k + 1) -P s (k + 1) K s (k) (12) P s (k, k + 1) =P s (k + 1)K s (k) (13)
Here, P f (k), P a (k), P s (k) and P s (k, k + 1) denote respectively the covariance matrices of the forecasted state x f (k), the filtered state x a (k), the smoother state x s (k) and the pair {x s (k), x s (k + 1)}. Finally, note that K f and K s are the filter and smoother Kalman gains and the innovation is denoted as d.

Innovation-based methods

The importance of the innovation statistics has been emphasized in the DA community by [START_REF] Daley | Estimating Model-Error Covariances for Application to Atmospheric Data Assimilation[END_REF] and [START_REF] Dee | On-line Estimation of Error Covariance Parameters for Atmospheric Data Assimilation[END_REF]. The "classic innovation" d, difference between the observations and the forecasted states in the observation space, defined in Eq. ( 6), implicitly takes into account the Q and R covariances. Unfortunately, as explained in [START_REF] Blanchet | A comparison of adaptive kalman filters for a tropical pacific ocean model[END_REF], by using only current observations, their individual contributions cannot be easily disentangled. Thus, the approaches using only the classic innovations are not studied in this review. Two main approaches were proposed in the literature to tackle this issue. They are based on the idea of producing multiple equations involving Q and R. The first one uses different innovation statistics in the observation space. The second one is based on lag-innovations or differences between consecutive innovations. From a statistical point of view, the innovationbased methods are "methods of moments", where we construct a system of equations which links various moments of the innovations with the parameters and then replace theoretical moments by the empirical ones in these equations.

a. Innovation statistics in the observation space [START_REF] Desroziers | Diagnosis of observation, background and analysis-error statistics in observation space[END_REF] proposed to examine various innovation statistics in the observation space. This method is now popular in the DA community. It is based on different innovation statistics between observations, forecasts and analysis, and all of them defined in the observation space: namely, d o-f (k) = y(k) -H x f (k) as in Eq. ( 6) and d o-a (k) = y(k) -H (x a (k)). We remark that another diagnostic using the difference between analysis x a (k) and reanalysis x s (k) has been proposed by [START_REF] Todling | A lag-1 smoother approach to system-error estimation: sequential method[END_REF] and [START_REF] Bowler | On the diagnosis of model error statistics using weak-constraint data assimilation[END_REF] in the case of sequential and variational DA respectively to estimate the covariance Q alone. In theory, in the linear and Gaussian case, the Desroziers innovation statistics should verify the equalities:

     E d o-f (k)d o-f (k) = H(k)P f (k)H(k) + R(k) (14) E d o-a (k)d o-f (k) = R(k) (15)
with E the expectation operator. In this approach, we do not estimate Q directly which is implicitly taken into account in P f . Instead, the approach attempts to compensate in P f for the lack of knowledge of Q as well as the systematic variance underestimation. This method is referred to as "covariance inflation". In practice, when using for instance EnKF with a small ensemble size, the spread is most of the time underestimated and this leads to filter divergence (see e.g. [START_REF] Carrassi | Data Assimilation in the Geosciences: An overview on methods, issues and perspectives[END_REF], their appendix A). Thus, covariance inflation can be required in a perfect model scenario (i.e. without η), because of sampling errors. For imperfect models, both sampling errors and an inappropriate representation of model errors lead to an underestimation of forecast ensemble spread and thus to filter divergence, see [START_REF] Raanes | Adaptive covariance inflation in the ensemble Kalman filter by Gaussian scale mixtures[END_REF].

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, see [START_REF] Anderson | A Monte Carlo Implementation of the Nonlinear Filtering Problem to Produce Ensemble Assimilations and Forecasts[END_REF]. Adaptive procedures to estimate this coefficient have been proposed by [START_REF] Wang | A Comparison of Breeding and Ensemble Transform Kalman Filter Ensemble Forecast Schemes[END_REF], Li et al. (2009a), [START_REF] Miyoshi | The Gaussian Approach to Adaptive Covariance Inflation and Its Implementation with the Local Ensemble Transform Kalman Filter[END_REF] and [START_REF] Bocquet | Ensemble Kalman filtering without the intrinsic need for inflation[END_REF] in the case of innovation statistics in the observation space. In the additive case, the diagonal of the forecast and/or analysis empirical covariance matrices is increased [START_REF] Mitchell | An adaptive ensemble Kalman filter[END_REF], [START_REF] Corazza | Use of the breeding technique to estimate the structure of the analysis errors of the day[END_REF], [START_REF] Whitaker | Ensemble Data Assimilation with the NCEP Global Forecast System[END_REF] and [START_REF] Houtekamer | Model Error Representation in an Operational Ensemble Kalman Filter[END_REF]). In the relaxation-to-prior case, [START_REF] Zhang | Impacts of Initial Estimate and Observation Availability on Convective-Scale Data Assimilation with an Ensemble Kalman Filter[END_REF] blended the forecast and analysis ensemble perturbations whereas [START_REF] Whitaker | Evaluating Methods to Account for System Errors in Ensemble Data Assimilation[END_REF] multiplied the analysis ensemble spread to relax the reduction of the spread, without blending perturbations. Finally, [START_REF] Bocquet | Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems[END_REF], [START_REF] Ying | An adaptive covariance relaxation method for ensemble data assimilation[END_REF] and [START_REF] Kotsuki | Adaptive covariance relaxation methods for ensemble data assimilation: Experiments in the real atmosphere[END_REF] proposed methods to adaptively estimate the relaxation parameters using innovation statistics. Adaptive covariance inflations are online estimation methods directly plugged to classic filtering method (like EKF here), with almost no additional computational cost. In practice, the use of this technique does not necessarily imply an additive error term η in Eq. ( 1). Thus, it is not a direct estimation of Q but an inflation applied to P f in order to compensate model uncertainties and sampling errors in EnKFs (see [START_REF] Raanes | Adaptive covariance inflation in the ensemble Kalman filter by Gaussian scale mixtures[END_REF], Sect. 4 and appendix C for more details). Several DA systems work with an inflation method and used it for its simplicity, low-cost and efficiency.

Here, we focus on the straightforward online estimation of a multiplicative inflation factor λ of the badly scaled P f (k) so that the corrected forecasted covariance is given by P f (k) = λ (k) P f (k). The estimate of the inflation factor is given by taking the trace of Eq. ( 14):

λ (k) = E d o-f (k) d o-f (k) -Tr (R(k)) Tr H(k) P f (k)H(k) . ( 16
)
The use of temporal smoothing for the online estimation of λ (k) is crucial in operational procedures and [START_REF] Miyoshi | The Gaussian Approach to Adaptive Covariance Inflation and Its Implementation with the Local Ensemble Transform Kalman Filter[END_REF] proposed augmenting the state vector with the inflation factor whose evolution is governed by a random walk equation. In this case, we need to specify an additional parameter for the variance term of this random walk, denoted by σ 2 λ . This parameter has to be carefully tuned to avoid the divergence of the filter. Then, at each time step k, when sufficient observations are available, an estimate of R(k) is directly given by Eq. ( 15). Note that Li et al. (2009a) proposed to estimate each component of a diagonal R matrix, and also suggested to use an offline procedure to compute the average of these variance terms.

b. Lag-innovation between consecutive times

Another way to estimate error covariances is to use multiple equations involving Q and R exploiting crosscorrelations between lag-innovations, i.e. the current d(k) and past classic innovations d(k -1), . . . , d(kl). For instance, considering the lag-zero and lag-one innovations, the following equations are satisfied in the linear and Gaussian case:

         E d(k)d(k) = H(k)P f (k)H(k) + R(k) = Σ(k) (17) E d(k)d(k -1) = H(k)M(k)P f (k -1)H(k) -H(k)M(k)K f (k -1)Σ(k -1). ( 18
)
Lag-innovations were introduced by [START_REF] Mehra | On the Identification of Variances and Adaptive Kalman Filtering[END_REF] in order to simultaneously recover the error covariance matrices for a Gaussian and linear state-space model. Mehra established analytic exact relations between Q and R, and the probabilistic expectations of d(k)d(kl) for linear systems in steady state. Then, [START_REF] Bélanger | Estimation of noise covariance matrices for a linear time-varying stochastic process[END_REF] extended these results to the case of time-varying linear stochastic processes, taking d(k)d(kl) as "observations" of Q and R and using a secondary Kalman filter to update them iteratively. As pointed out in [START_REF] Bélanger | Estimation of noise covariance matrices for a linear time-varying stochastic process[END_REF], this method would no longer be analytically exact if the error matrices are updated adaptively at each time step. Later, [START_REF] Dee | An Efficient Algorithm for Estimating Noise Covariances in Distributed Systems[END_REF] proposed a computationally cheaper algorithm for the Bélanger's method. More recently, authors focused on high-dimensional and nonlinear systems using the EKF and EnKF: [START_REF] Berry | Adaptive ensemble Kalman filtering of non-linear systems[END_REF] proposed a fast algorithm based on Mehra's method and Harlim et al. ( 2014) followed the original Bélanger's algorithm. [START_REF] Zhen | Adaptive error covariances estimation methods for ensemble Kalman filters[END_REF] proposed a modified version of Bélanger's method and compared it to the [START_REF] Berry | Adaptive ensemble Kalman filtering of non-linear systems[END_REF] approach.

Here, we briefly introduce the algorithm of Berry and Sauer (2013) using lag-one innovations. It is based on the online (or adaptive) estimation of Q(k) and R(k), which satisfy the following relations in the linear and Gaussian case:

P(k) =M(k) -1 H(k) -1 d(k)d(k -1) H(k) - + K f (k -1)d(k -1)d(k -1) H(k) - (19) Q(k) = P(k) -M(k -1)P a (k -1)M(k -1) - (20) R(k) =d(k -1)d(k -1) -H(k)P f (k -1)H(k) -. (21)
In this online procedure, joint estimations of Q(k) and R(k) can abruptly vary over time. Thus, the temporal smoothing of the covariances being estimated becomes crucial. As suggested by [START_REF] Berry | Adaptive ensemble Kalman filtering of non-linear systems[END_REF], an exponential smoothing between current and past estimates is a reasonable choice,

Q(k + 1) = Q(k) + ( Q(k) -Q(k))/τ, (22)
R(k + 1) = R(k) + ( R(k) -R(k))/τ (23)
with τ the smoothing parameter and started from Q(0). When τ is small, weight is given to the current estimate Q and when τ is larger it gives a smoother sequence Q. As pointed out by [START_REF] Zhen | Adaptive error covariances estimation methods for ensemble Kalman filters[END_REF], usually a large value of τ is chosen to avoid numerical instability.

It is worth pointing out that in the case of sparse observations, the estimate of P in Eq. ( 19) might be underdetermined, even if the system is observable. This is attributed to the use of only one lag-innovation. Theoretically, all components of Q should be identifiable if the system is observable and more lag-innovations are used. But in practice, using more lag-innovations implies increased computational cost and does not necessary lead to accurate estimates. [START_REF] Zhen | Adaptive error covariances estimation methods for ensemble Kalman filters[END_REF] compared the modified version of Bélanger's method with different choices of maximal lags and found that a maximal lag of 4 is optimal in a specific numerical example on Lorenz-96 defined in [START_REF] Lorenz | Predictability: A problem partly solved[END_REF].

Likelihood-based methods

The likelihood approaches were put forward in the DA community by [START_REF] Dee | On-line Estimation of Error Covariance Parameters for Atmospheric Data Assimilation[END_REF], [START_REF] Blanchet | A comparison of adaptive kalman filters for a tropical pacific ocean model[END_REF] as well as [START_REF] Mitchell | An adaptive ensemble Kalman filter[END_REF] where it was proposed to maximize the likelihood of the innovation, i.e. p (y(k)|y(k -1)), defined by the mean vector d(k) computed in Eq. ( 6) and covariance matrix Σ(k) introduced in Eq. ( 17) and also used in the computation of the Kalman filter gain in Eq. ( 7). Unfortunately, they reach the same conclusions than for the innovation-based methods, i.e. the joint estimation of Q and R is not straightforward if we use only the current observations. To tackle this issue, several methods have been proposed recently. The first one is to write the estimation problem using a Bayesian framework, and jointly estimate prior distributions of Q and R parameters with the innovation likelihood. The second one is to maximize the so-called "total likelihoods", i.e. taking into account the innovation likelihoods of each time step or taking into account the global structure of the sate-space model for all the time steps.

a. Bayesian inference

In a Bayesian approach, we assume that the elements of Q and R covariance matrices have a priori distributions which are controlled by some hyperparameters. In practice, it is difficult to have a prior distribution 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, and we denote the corresponding parameters as θ . Then, we jointly and adaptively estimate the state x and parameters θ using Bayes' theorem:

p (x(k), θ (k)|y(1 : k)) = p (x(k)|y(1 : k), θ (k)) p (θ (k)|y(1 : k)) . (24)
In Eq. ( 24), p (x(k)|y(1 : k), θ (k)) is given by filtering DA algorithms and we approximate recursively p (θ (k)|y(1 : k)) using the the likelihood of the innovations p (y(k)|y(1 : k -1), θ (k)) as

p (θ (k)|y(1 : k)) ∝ p (y(k)|y(1 : k -1), θ (k)) p (θ (k)|y(1 : k -1)) . (25)
Bayesian approaches have been applied in the atmospheric chemistry community and reviewed by [START_REF] Michalak | Maximum likelihood estimation of covariance parameters for Bayesian atmospheric trace gas surface flux inversions[END_REF] and [START_REF] Wu | Hyperparameter estimation for uncertainty quantification in mesoscale carbon dioxide inversions[END_REF]. [START_REF] Purser | A Bayesian technique for estimating continuously varying statistical parameters of a variational assimilation[END_REF] introduced the Bayesian approach in variational DA for the estimation of two statistical parameters, controlling the magnitude of the variance and the spatial dependencies in Q, assuming that R is known and using a univariate model. Then, [START_REF] Stroud | Sequential State and Variance Estimation within the Ensemble Kalman Filter[END_REF] used a similar approach combined with EnKF in the Lorenz-96 model for the estimation of a common multiplicative scalar parameter for predefined matrices Q and R. In that case, the scalar parameter affects simultaneously the Q and R matrices. Considering the experiments about the importance of ||Q||/||R|| ratio presented in Fig. 1, we can guess that this approach is maybe not optimal. Then, other works have applied similar Bayesian approaches for the estimation of parameters governing the shape of R only: [START_REF] Frei | Sequential State and Observation Noise Covariance Estimation Using Combined Ensemble Kalman and Particle Filters[END_REF] in the Lorenz-96 system, [START_REF] Winiarek | 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[END_REF][START_REF] Winiarek | 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[END_REF] assimilating nuclear pollutants using a regional atmospheric model (in this case, R partially accounts for model error), [START_REF] Ueno | Bayesian estimation of the observation-error covariance matrix in ensemble-based filters[END_REF] using two linear shallow-water equations to assimilate satellite altimetry. By contrast, [START_REF] Solonen | Estimating model error covariance matrix parameters in extended Kalman filtering[END_REF] proposed a Bayesian approach for the estimation of Q only, assuming that the R matrix is known, in a two-layer quasi-geostrophic model. Finally, [START_REF] Stroud | A Bayesian adaptive ensemble Kalman filter for sequential state and parameter estimation[END_REF] tested their estimation method on different spatio-temporal systems with a joint estimation of Q and R. Their experiments on the Lorenz-96 system assumed no model error Q.

The Bayesian inference approach is an online estimation procedure with joint estimation of the state of the system and hyperparameters θ (k), controlling shape parameters of the Q(k) and R(k) error covariance matrices. In terms of the hidden state, this corresponds to a hierarchical Bayesian approach so that an ensemble of filters may be required to determine the posterior distribution in Eq. ( 24). In practice, this method may required a large number of Monte-Carlo simulations to estimate correctly p (θ (k)|y(1 : k)) defined in Eq. ( 25) or iterative procedures as in [START_REF] Ueno | Bayesian estimation of the observation-error covariance matrix in ensemble-based filters[END_REF]. Alternatively, Scheffler et al. ( 2018) assume a Gaussian distribution for p (θ (k)|y(1 : k -1)) in Eq. ( 25) and use two nested EnKFs, reducing the computational cost of this hierarchical Bayesian procedure. In principle, the Bayesian approach is able to estimate time dependent hyperparameters, however recent works assume that θ (k) is constant, not depending in time. The potential of this framework to estimate time dependent covariances is an interesting topic to be addressed. The joint estimation of parameters controlling separately Q and R still remains a challenge. However relevant parametric shapes of covariance matrices, such as the Matérn covariance model for R, have been proposed as in [START_REF] Stroud | A Bayesian adaptive ensemble Kalman filter for sequential state and parameter estimation[END_REF].

b. Maximization of the total likelihoods

The innovation likelihood at time k is defined by p (y(k)|y(1 : k -1), θ (k)) in Eq. ( 25). Maximizing this likelihood at each time step has been proposed by various authors: Dee (1995), [START_REF] Blanchet | A comparison of adaptive kalman filters for a tropical pacific ocean model[END_REF], [START_REF] Zheng | An adaptive estimation of forecast error covariance parameters for Kalman filtering data assimilation[END_REF], [START_REF] Mitchell | An adaptive ensemble Kalman filter[END_REF] and [START_REF] Liang | Maximum likelihood estimation of inflation factors on error covariance matrices for ensemble Kalman filter assimilation[END_REF]. But the maximization of the innovation likelihood has two main issues. Firstly, the innovation covariance matrix Σ(k) = H(k)P f (k)H(k) + R(k) is mixing the information of R and Q, contained in P f . When using only time k, it is difficult to identify the model and observation error covariances and in practice, authors only estimated one of them. Secondly, the number of observations at each time step is in general limited and as pointed out by [START_REF] Dee | On-line Estimation of Error Covariance Parameters for Atmospheric Data Assimilation[END_REF], available observations should exceed "the number of tunable parameter by two or three orders of magnitude". For these reasons, a reasonable alternative is to use a batch of observations distributed in time and assume θ being constant in time. The resulting total likelihood expressed sequentially through conditioning is given by

p (y(1 : K)|θ ) = K ∏ k=2 p (y(k)|y(1 : k -1), θ ) . ( 26
)
This likelihood is said to be "incomplete" because it only depends on the observations of the state-space system, not the hidden state. But since this is an integration innovation likelihood over a long period of time, this gives more information to estimate Q and R or related parameters. This likelihood is said to be marginal since it results from marginalizing the hidden state at a given observation time. The incomplete total likelihood is also a useful tool to evaluate the quality of model forecasts or how well they match the observations, considering both model and observation uncertainties. [START_REF] Hannart | DADA: data assimilation for the detection and attribution of weather and climate-related events[END_REF] and [START_REF] Carrassi | Estimating model evidence using data assimilation[END_REF] used it for model evidence, when various models are in competition. The maximization of the incomplete total likelihood given in Eq. ( 26) has been applied to linear and nonlinear systems in the estimation of deterministic and stochastic parameters (related to Q) in Delsole and Yang (2010) using a direct sequential optimization procedure. Then, for nonlinear dynamics, Ueno et al. ( 2010) used a grid-based procedure to estimate noise levels and spatial correlation lengths of Q and a level noise for R. This grid-based method used predefined sets of covariance parameters and test the different combinations to find the one that maximizes the likelihood criterion. [START_REF] Brankart | Efficient Adaptive Error Parameterizations for Square Root or Ensemble Kalman Filters: Application to the Control of Ocean Mesoscale Signals[END_REF] proposed also a method using the same criterion additionally with an initial information on scale and correlation length parameters of Q and R. This information is only given at the first time and progressively forgotten with time using a decreasing exponential factor. The marginalization of the hidden state in Eq. ( 26) is considering all the previous observations and in practice it requires the use of a filter. The maximization of the total incomplete likelihood using the EnKF to estimate model error covariance Q was conducted in [START_REF] Pulido | Stochastic parameterization identification using ensemble Kalman filtering combined with maximum likelihood methods[END_REF] where they used a gradient-based optimization technique.

Authors also proposed to work on the maximization of the total likelihood using the marginalization of the whole trajectory of the hidden state from k = 0 to K. In that case, we talk about the "complete" total likelihood or joint density of the observations and the hidden state, expressed as

p (y(1 : K), x(0 : K)|θ ) = p (x(0)) K ∏ k=1 p (x(k)|x(k -1), θ ) K ∏ k=1 p (y(k)|x(k), θ ) (27)
where the three terms on the rhs are related to the initial state, the state equation in Eq. ( 1) and the observation equation in Eq. ( 2). In practice, the marginalization of the full hidden state from k = 0 to K is not possible, so that the complete total likelihood cannot be evaluated directly, see explanations in [START_REF] Pulido | Stochastic parameterization identification using ensemble Kalman filtering combined with maximum likelihood methods[END_REF]. Therefore, [START_REF] Shumway | An Approach to Time Series Smoothing and Forecasting using the EM Algorithm[END_REF] proposed to use an iterative procedure, requiring the use of a smoother, to maximize the likelihood criterion given in Eq. ( 27). They used the Expectation-Maximization algorithm (hereinafter noted EM, see [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF] for more details) to ensure the convergence to the maximum likelihood estimator, and applied it to estimate Q and R in the case of linear dynamics. In DA, the EM algorithm has been implemented for estimating only R in [START_REF] Ueno | Iterative algorithm for maximumlikelihood estimation of the observation-error covariance matrix for ensemble-based filters[END_REF] for the Zebiak and Cane (ZC) model and satellite altimetry, for the background covariance matrix B and R in [START_REF] Liu | Uncertainty quantification of pollutant source retrieval: comparison of Bayesian methods with application to the Chernobyl and Fukushima Daiichi accidental releases of radionuclides[END_REF] in case of accidental pollutant source retrieval, for both Q and R with an orographic subgrid-scale nonlinear observation operator in [START_REF] Tandeo | Offline parameter estimation using EnKF and maximum likelihood error covariance estimates: Application to a subgrid-scale orography parametrization[END_REF], and for the Lorenz-63 in [START_REF] Dreano | Estimating model error covariances in nonlinear state-space models using Kalman smoothing and the expectationmaximisation algorithm[END_REF]. Recently, [START_REF] Pulido | Stochastic parameterization identification using ensemble Kalman filtering combined with maximum likelihood methods[END_REF] used the EM algorithm and compared with a gradient ascent optimization of the incomplete total likelihood to estimate Q along with the deterministic and stochastic physical parameters in the one and two scales Lorenz model described in [START_REF] Lorenz | Optimal Sites for Supplementary Weather Observations: Simulation with a Small Model[END_REF]. Finally, [START_REF] Chau | An efficient particle-based method for maximum likelihood estimation in nonlinear state-space models[END_REF] combined conditional particle filters and the EM algorithm for the joint estimation of Q and R and show the improvement compared to Kalmanbased filters.

The EM algorithm considers the total period of time k = 0 to K to maximize complete the total likelihood given in Eq. ( 27). Thus, it leads to an offline estimation of the error covariance matrices or related parameters (i.e. constant, non-adaptive). In the expectation step, we evaluate the expected likelihood given in Eq. ( 27) conditionally to the previous estimates of Q and R as well as the total observations {y(1), ..., y(K)}. This leads to the use of Kalman smoother procedures to estimate x s and P s at the different times. In the maximization step, Q and R are updated using the following estimators:

Q = 1 K K ∑ k=1 {(x s (k) -M x s (k -1)) (x s (k) -M x s (k -1)) + M(k)P s (k -1)M(k) + P s (k) -P s (k -1, k)M(k) -M(k)P s (k -1, k) } (28) R = 1 K K ∑ k=1 {(y(k) -H(k)x s (k)) (y(k) -H(k)x s (k)) + H(k)P s (k)H(k) }. (29)
In practice, the total and complete likelihood given in Eq. ( 27) cannot be evaluated exactly. Thus, to evaluate the performance at each iteration of the EM procedure, we compute the incomplete total likelihood given in Eq. ( 26) and this criterion increases along the iterations of the EM algorithm, see [START_REF] Wu | On the convergence properties of the EM algorithm[END_REF]. The computational cost of the EM algorithm is large because it requires the use of Kalman-based filter and smoother at each iteration. Nevertheless, this EM method does not require any additional parameter and is robust to initial conditions, i.e. the values given for Q and R matrices in the first EM iteration. Additionally, note that the maximization of the total likelihood allows the estimation of the initial state vector x b and covariance matrix B as discussed in [START_REF] Tandeo | Offline parameter estimation using EnKF and maximum likelihood error covariance estimates: Application to a subgrid-scale orography parametrization[END_REF] and [START_REF] Dreano | Estimating model error covariances in nonlinear state-space models using Kalman smoothing and the expectationmaximisation algorithm[END_REF]. In that case, we should write p (x(0)|θ ) instead of p (x(0)) in Eq. ( 27).

Other methods

In this section, we describe other methods which have been used in the past or methods that are more diagnostic tools than direct estimations of error covariance matrices. We also include some relevant references on inverse problems in environmental data.

a. State augmentation

State augmentation was first proposed in [START_REF] Schmidt | Applications of state space methods to navigation problems[END_REF] and is known as the Schmidt-Kalman filter. Then, [START_REF] Jazwinski | Stochastic processes and filtering theory[END_REF] proposed some extensions. The idea is to augment the state vector in order to estimate both the state of the system and additional parameters among which the bias, forcing terms, physical parameters, and finally error covariances as in [START_REF] Zupanski | A General Weak Constraint Applicable to Operational 4DVAR Data Assimilation Systems[END_REF] and [START_REF] Tremolet | Model-error estimation in 4D-Var[END_REF]. In these works, authors create cross-correlation between the state of the system and the additional parameters. The method works only for parameters strongly related to the state of the system, such as physical parameters, see [START_REF] Ruiz | Estimating Model Parameters with Ensemble-Based Data Assimilation: A Review[END_REF]. However, [START_REF] Stroud | Sequential State and Variance Estimation within the Ensemble Kalman Filter[END_REF] as well as [START_REF] Delsole | State and parameter estimation in stochastic dynamical models[END_REF] formally demonstrated that augmentation methods fail for variance parameters and thus for Q and R. Indeed, a critical aspect of this approach is that one needs to define an evolution model for the covariance parameters. This is a difficult task, and often persistence is used, which means that the estimate and the associate error variance only change at analysis times, and the estimated variance is thus bound to decrease in time: it is reduced during the forecast step using a persistence model. This makes the use of random walk or inflation mandatory, or a change in the parametric error dynamics such as in Carrassi and Vannitsem (2011a).

b. Analysis increment approach

Analysis increments refers to statistical methods that study the relationship between two consecutive times of a dynamical system. The use of regression methods has been firstly proposed by [START_REF] Lorenz | An Experiment in Nonlinear Statistical Weather Forecasting[END_REF] and then by [START_REF] Leith | Objective Methods for Weather Prediction[END_REF] to learn error statistics of dynamical models in meteorology (e.g. bias and covariance Q of the error η). Then, this approach was first discussed in the context of DA by Li et al. (2009b), and it was then further expanded by [START_REF] Carrassi | Accounting for Model Error in Variational Data Assimilation: A Deterministic Formulation[END_REF] in the context of variational DA with time correlated model error. The same reanalysis increment approach has been used in Carrassi and Vannitsem (2011b) to estimate model error due to unresolved scale and later applied in the context of a deterministic EnKF by [START_REF] Mitchell | Accounting for model error due to unresolved scales within ensemble Kalman filtering[END_REF].

c. Covariance matching

The covariance matching method has been introduced by [START_REF] Fu | Fitting Dynamic Models to the Geosat Sea Level Observations in the Tropical Pacific Ocean. Part II: A Linear, Wind-driven Model[END_REF]. It consists in matching sample covariance matrices to their theoretical expectations. It corresponds to a method of moments applied to different innovations. Thus, it is very similar to [START_REF] Desroziers | Diagnosis of observation, background and analysis-error statistics in observation space[END_REF] method, except that covariance matching is performed on a set of historical observations and numerical simulations, not online in a DA scheme. It has been extended by [START_REF] Menemenlis | Error Estimates for an Ocean General Circulation Model from Altimeter and Acoustic Tomography Data[END_REF] to time-lagged innovations as in [START_REF] Bélanger | Estimation of noise covariance matrices for a linear time-varying stochastic process[END_REF] and they also relaxed the assumption of independence between the true state and model simulation errors.

d. The χ 2 , cross-validation and whiteness of laginnovations

These methods are not direct estimation of error covariance matrices but diagnostic tools. The first one is a statistical test that examines the variance of the normalized innovations that follow in theory a χ 2 distribution with a given number of degrees of freedom. As pointed out by respectively by [START_REF] Michalak | Maximum likelihood estimation of covariance parameters for Bayesian atmospheric trace gas surface flux inversions[END_REF] and [START_REF] Rayner | Reconstructing the recent carbon cycle from atmospheric CO2, δ 13C and O2/N2 observations[END_REF], lots of combination of Q and R errors can lead to the acceptation of this test and this method cannot guide the relative allocation of error between the 2 ones. Crossvalidation is a classic bootstrap strategy to test the accuracy in statistical modeling. It is based on the repetition of validations between a learned model from a random training dataset and test on the rest of the observations, see [START_REF] Wu | Hyperparameter estimation for uncertainty quantification in mesoscale carbon dioxide inversions[END_REF] for more details. The third one was pointed by [START_REF] Jazwinski | Stochastic processes and filtering theory[END_REF] and consists in evaluating the properties of lag-innovations. They are supposed to be white in time in case of optimal filtering, i.e. using appropriates Q and R matrices.

Summary, conclusions and perspectives

In this review paper, we presented different methods to estimate the error covariances in data assimilation. As usually stated in data assimilation, we assume that model and observation errors are additive and Gaussian with centered null mean and covariance matrices noted Q and R. The individual and joint impacts of badly calibrated covariances were firstly presented using a linear toy-model. The experiments clearly showed that for reaching reasonable results of the filter, in terms of pure reconstruction using root mean squared error, the joint estimation of both components is a crucial point. We also highlighted the impact on the coverage probability, related to the estimated covariance of the reconstructed state and thus to the uncertainty quantification in data assimilation.

Summary of existing methods

We focused on four main methodologies for the joint estimation of Q and R used in data assimilation. They are summarized in Table 1. We first dealt with methods based on innovations, i.e. the difference between observations and forecasted state, and the use of empirical and theoretical moments, also known as method of moments in statistics. We presented the method of innovation statistics in the observation space by [START_REF] Desroziers | Diagnosis of observation, background and analysis-error statistics in observation space[END_REF]. This approach is often associated with inflation methods where model error covariance Q is not necessary needed and where an estimated inflation factor artificially increases the forecast error covariance P f . This online method is low-cost and adaptive but the inflation factor controls the global amplitude of the covariance and cannot adapt to specific parts of it. Nevertheless, the covariance inflation is used in lot of operational data assimilation systems and is very common with lot of implementations like the multiplicative one, see Li et al. (2009a) and [START_REF] Miyoshi | The Gaussian Approach to Adaptive Covariance Inflation and Its Implementation with the Local Ensemble Transform Kalman Filter[END_REF] for instance. Another approach received recently specific attention in the data assimilation context: the laginnovation exploiting the autocorrelation of the innovation between consecutive times. It has been introduced by [START_REF] Mehra | On the Identification of Variances and Adaptive Kalman Filtering[END_REF] and [START_REF] Berry | Adaptive ensemble Kalman filtering of non-linear systems[END_REF] extended this lag-one innovation method to the nonlinear case in data assimilation, whereas [START_REF] Harlim | An ensemble Kalman filter for statistical estimation of physics constrained nonlinear regression models[END_REF] implemented a lag-l innovation method following [START_REF] Bélanger | Estimation of noise covariance matrices for a linear time-varying stochastic process[END_REF] idea. These lag-innovation techniques are adaptive, online and plugged into any filter with a moderate additional computational cost. However, results are very sensitive to a tuning parameter τ, used to smooth the estimated covariance matrices along time and avoid the method from breaking down. Lag-innovation methods have been tested on linear and Lorenz-96 systems, not yet on real data assimilation schemes. At this stage, authors considered constant Q and R but in practice, when τ parameter is correctly tuned, lag-innovations can deal with time varying matrices.

The two last methods summarized in Table 1 are based on the maximization of the likelihood criterion. [START_REF] Dee | On-line Estimation of Error Covariance Parameters for Atmospheric Data Assimilation[END_REF] firstly pointed out the importance of maximizing the innovation likelihood, but the estimation of both Q and R using only current observations is limited. Thus, authors then proposed Bayesian inferences to jointly maximize the innovation likelihood and the likelihood of parameters of the error covariance matrices. These parameters are assumed to follow prescribed prior distributions that have to be carefully chosen. Moreover, large Monte-Carlo simulations are needed to estimate the hyperparameters of these distributions. The joint estimation of Q and R has not been [START_REF] Stroud | Sequential State and Variance Estimation within the Ensemble Kalman Filter[END_REF] where they used the same inflation parameter for both matrices. These approaches have been tested with the Lorenz-96 system in [START_REF] Frei | Sequential State and Observation Noise Covariance Estimation Using Combined Ensemble Kalman and Particle Filters[END_REF] and with realistic models and observations in [START_REF] Winiarek | 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[END_REF][START_REF] Winiarek | 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[END_REF] as well as [START_REF] Ueno | Bayesian estimation of the observation-error covariance matrix in ensemble-based filters[END_REF]. Bayesian approaches are online methods and able to deal with time varying covariances. The last method is based on the maximization of the complete total likelihood as proposed in linear statespace by [START_REF] Shumway | An Approach to Time Series Smoothing and Forecasting using the EM Algorithm[END_REF]. Given its stability and robustness, we choose to detail the expectationmaximization procedure that iteratively aims at maximizing the complete total likelihood, even in nonlinear statespace models as shown in [START_REF] Dreano | Estimating model error covariances in nonlinear state-space models using Kalman smoothing and the expectationmaximisation algorithm[END_REF]. In comparison to the previous three methods, this one is offline and not adaptive, assuming that Q and R are constant over a batch period. It has a high-computational cost but it does not require any additional tuning parameter. So far, recent works jointly estimate model and observation error matrices in toy-models like Lorenz-63 and Lorenz-96 with one and two scales, respectively in [START_REF] Dreano | Estimating model error covariances in nonlinear state-space models using Kalman smoothing and the expectationmaximisation algorithm[END_REF] and [START_REF] Pulido | Stochastic parameterization identification using ensemble Kalman filtering combined with maximum likelihood methods[END_REF]. [START_REF] Ueno | Iterative algorithm for maximumlikelihood estimation of the observation-error covariance matrix for ensemble-based filters[END_REF] applied the EM algorithm for the estimation of R in a realistic coupled atmosphere-ocean model and [START_REF] Liu | Uncertainty quantification of pollutant source retrieval: comparison of Bayesian methods with application to the Chernobyl and Fukushima Daiichi accidental releases of radionuclides[END_REF] applied EM algorithm to a radionuclide transport model using real observations.

Remaining challenges

There are still remaining challenges for the four methods detailed in this review. The first concerns the improvements of online and adaptive techniques regarding additional parameters that control the variations of Q and R estimates in time. Instead of using fixed values for these parameters, for instance τ in lag-innovations or σ 2 λ in inflation methods, we suggest to use time-dependent adaptations. This will avoid the problems of instabilities close to the solution. Another option could be to adapt these online procedures to the offline case, working with very stable parameters values (τ high, σ 2 λ low) and iterate the procedures on a batch of observations as in the EM algorithm. This was suggested and rapidly tested in [START_REF] Desroziers | Diagnosis of observation, background and analysis-error statistics in observation space[END_REF] with encouraging results. To the best of our knowledge, it has not been yet tested with lag-innovation methods.

The second challenge concerns the Bayesian approach where joint estimation of Q and R seems problematic. As pointed out by [START_REF] Berry | Correlation between system and observation errors in data assimilation[END_REF], correlation between model and observation error terms is in practice highly probable in real data assimilation problems. In this case, instead of using independent prior distributions, the use of joint prior distributions for parameters of Q and R might physically constraint the optimization procedure.

A third challenge concerns the offline EM algorithm using the total likelihood over a large batch of observations. This procedure can be adapted to account for time varying error covariances. A simple way is to work on small independent sets of observations and apply various EM procedures. Then, the Q and R could be smoothed in time. Another way is to apply online EM algorithms (see for instance [START_REF] Cappé | Online Expectation-Maximisation. Mixtures: Estimation and Applications[END_REF]) with the likelihood averaged locally in time. Note also that EM algorithm, whether for the online or offline case, can be coupled with direct optimization methods like Newton-Raphson to speedup convergence.

From our point of view, the last challenge concerns the estimation of other statistical parameters of the state-space model given in Eqs. (1-2) and associated filters. Indeed, the initial condition x b and B are crucial for certain satellite retrieval problems and have to be estimated, principally in offline cases where filtering and smoothing are repeated on various iterations. Finally, estimation methods should also consider the estimation of systematic or time varying biases, corresponding to the deterministic part of η and ε. It has been initially proposed by Dee et al. (1999a) and tested in Dee et al. (1999b) in the case of the maximization of the innovation likelihood and re-cently adapted with a Bayesian formulation in [START_REF] Liu | Uncertainty quantification of pollutant source retrieval: comparison of Bayesian methods with application to the Chernobyl and Fukushima Daiichi accidental releases of radionuclides[END_REF] as well as [START_REF] Berry | Correcting Biased Observation Model Error in Data Assimilation[END_REF].

Perspectives

Beyond these possible improvements in techniques, we also discuss prospects for future work. A first perspective concerns the combination of methods, especially regarding the estimation modes, i.e. online or offline. In our opinion, a great advantage of the expectationmaximization algorithm is the absence of additional parameter and its robustness, due to the number of observations used to approximate the total state-space likelihood, potentially large when using an important batch period. Thus, it is a useful tool to get a first estimate and averaged form for Q and R matrices, i.e. to infer mean amplitudes and parametric shapes including bloc structures and spatial dependencies. Since the computational cost for this offline method is important, this calibration part has to be done once during the configuration stage of the assimilation system. Then, low-cost online methods, e.g. based on lag-innovations, can use the robust offline estimates as initial conditions and adapt to slow variations of Q(k) and R(k) using relevant parametric shapes of these matrices.

In a following work, we plan to evaluate offline and online methods to different cases, considering both constant and time varying Q and R matrices. First, we will compare the different methodologies to the linear and unidimensional model used in this review in order to evaluate their performance and convergence in the asymptotic case. We will use the root mean squared error and the coverage probability to measure the inference performance on mean and covariance. Then, we will test the estimation methods on the chaotic Lorenz-96 model and evaluate their performance and robustness varying the number of available observations. We will finally implement methodologies to a more realistic case for operational data assimilation, applying for instance to a mid-complexity general circulation model and real or simulated satellite data. In this case, the number of observations will be limited compared to the size of the state. To reduce the degrees of freedom, the use of parametric shapes for error covariances will be necessary to tackle the rank-deficient observations. Moreover, the use of deterministic ensemble filters with localizations will be necessary in such realistic data assimilation problems.

Finally, the Gaussian and additive formulation of the error terms generally stated in data assimilation as in Eqs (1-2) is extremely convenient. In practice, it allows the application of Kalman-based algorithm and greatly simplify the use of method of moments and maximum likelihood approaches as detailed in this review. But is it able to compensate for non-additive sources of errors? Indeed, in realistic data assimilation problems, errors are multiplicative or introduced in the model by misparametrization and/or parameter evolution. The presented methods have to be tested on those configurations to evaluate whether or not the additive and Gaussian formulation with covariances Q and R is robust enough and identify its limitations.

  the state conditionally to past/present observations and past/present/future observations respectively.Here, we briefly remind the equations of the EKF and EKS based on the Rauch-Tung-Striebel (RTS) solution detailed in[START_REF] Cosme | Smoothing Problems in a Bayesian Framework and Their Linear Gaussian Solutions[END_REF]. They are divided in three main steps:

TABLE 1 .

 1 Comparison of several methods to estimate error covariances matrices Q and R in data assimilation.

	Described	Criteria	Estimation method	Estimation mode	Estimation of	Additional
	method			(computation cost)	covariance Q?	parameters?
	Covariance	Innovation statistics in	Method of moments	Online (low)	No (inflation λ	Yes (inflation
	inflation	the observation space			instead)	variance σ 2 λ )
	Lag-	Lag-innovation between	Method of moments	Online (moderate)	Yes	Yes (temporal
	innovation	consecutive times				smoother τ)
	Bayesian	Innovation & hyperpa-	Likelihood methods	Online (high)	No (or joint pa-	Yes (prior dis-
	inference	rameter likelihoods			rameter with R)	tributions for θ )
	Expectation-	Total likelihoods of the	Maximum likelihood Offline (high)	Yes	No
	maximization	state-space model
	evaluated, except in

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). A.C. has been funded by the project REDDA (#250711) of the Norwegian Research Council.