Yin Yang
email: iamyangyin@gmail.com

Keywords: Expectation-Maximization algorithm, Stochastic modeling, Parameter estimation, Data assimilation

Estimating the parameters of geophysical dynamic models is an important task in Data Assimilation (DA) technique used for forecast initialization and reanalysis. In the past, most parameter estimation strategies were derived by state augmentation, yielding algorithms that are easy to implement but may exhibit convergence difficulties. The Expectation-Maximization (EM) algorithm is considered advantageous because it employs two iterative steps to estimate the model state and the model parameter separately. In this work, we propose a novel ensemble formulation of the Maximization step in EM that allows a direct optimal estimation of physical parameters using iterative methods for linear systems. This departs from current EM formulations that are only capable of dealing with additive model error structures. This contribution shows how the EM technique can be used for dynamics identification problem with a model error parameterized as arbitrary complex form. The proposed technique is here used for the identification of stochastic subgrid terms that account for processes unresolved by a geophysical fluid model. This method, along with the augmented state technique, are evaluated to estimate such subgrid terms through high resolution data. Compared to the augmented state technique, our method is shown to yield considerably more accurate parameters. In addition, in terms of prediction capacity, it leads to smaller generalization error as caused by the overfitting of the trained model on presented data and eventually better forecasts.

Introduction

Parameter estimation is an important part in Data Assimilation (DA) technique. Traditionally the same estimation scheme is employed both for state and parameter estimation through the so-called augmented state technique. This technique consists in encapsulating the whole set of unknown parameters into an augmented state vector, which is then jointly estimated. Such works have been pioneered by [START_REF] Navon | Practical and theoretical aspects of adjoint parameter estimation and identifiability in meteorology and oceanography[END_REF] in the context of variational methods and extended by [START_REF] Anderson | An ensemble adjustment kalman filter for data assimilation[END_REF] to ensemble methods.

One of the crucial parameter for the forecast or analysis of geophysical flow is related to the inescapable lack of representativeness of the dynamics model. Those errors come from approximations made in the analytical constitution steps of the models and are dictated by pragmatic computational constraints. They may also ensue from a lack of physical information on the complex small-scale processes involved. The action on the large-scale (resolved) components of those unresolved processes must be properly modeled; they are usually qualified under the term of subgrid model. These models generally do not rely on rigorous theoretical grounds and are often defined to mimic qualitative descriptions of turbulence or of the small scales effect.

Recently, [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] and [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF] proposed to introduce a stochastic dynamics approach to model the unresolved small-scale components. The basic idea is to consider the stochastic nature of flow from the start in the physical modeling process to directly derive the governing physical dynamical system driven by stochastic processes. In these approaches, recovering the subgrid-scale errors boils down to the calibration of a variance tensor representing the fast-changing fluctuation associated to turbulence. Nowadays, observation snapshots are reaching a resolution that is much finer than the model grid resolution. As a consequence, the question on the ability to estimate the subgrid parameters from high resolution observation immediately arises.

In this study, we investigate such an estimation problem. We rely on an reconstruction basis formulated by coupling a stochastic modeling process [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] and relevant data through an optimal estimation scheme. Recently, [START_REF] Yang | High-resolution data assimilation through stochastic subgrid tensor and parameter estimation from 4DEnVar[END_REF] showed an implementation of the augmented state strategy by applying an ensemble optimal control (EnVar) method to this state/parameter estimation problem. The authors found that although the state estimated from the joint distribution can be quite accurate due to the enlarged ensemble spread induced by the stochastic dynamics representation, the estimated parameter field, on the contrary, tends to be much less accurate with an underestimation of the reference values. As a result, the parameter field requires to be retuned by adhoc inflation techniques to avoid model blow-up and/or covariance collapse. Model blow up results from strong uncontrolled numerical instabilities of the dynamical model, while ensemble collapsing correspond to ensemble members with no variabilities. An insufficient spread of the ensemble prevents in particular any accurate estimation of the probability density functions involved in the Bayesian inference process.

In this paper, we devise an Expectation-Maximization (EM) algorithm to overcome the above shortcomings associated to augmented-state-based techniques. Our motivation for using EM algorithm is twofold. Within a Bayesian framework, we want both to propose an elegant and flexible solution for estimating separately the state and the model parameters, and to provide a numerically stable formulation for maximizing the likelihood. Eventually, highly accurate stochastic parameters associated with the stochastic modeling process can be obtained, leading to better reconstruction and prediction of the state.

Before describing precisely the proposed approach, we briefly review the principles of the EM technique. In the statistical literature, the EM algorithm [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF]) is an efficient parameter estimation algorithm aiming at finding the maximum likelihood solution of a system containing latent variables. The latent variables (or hidden states) refer to non-observed variables that are usually introduced to facilitate the modeling process. The maximum likelihood (ML) estimation of parameter, θ, given observations, y, consists in finding θ so that the conditional probability, p(y|θ), reaches its maximum. For systems without latent variables, such estimation can be directly solved by employing Newton-type methods.

In terms of dynamical systems with latent variables, x, the maximum likelihood function, p(y|θ), is given by marginalizing the complete joint likelihood function, p(y, x|θ), over x: (y, x|θ).

ln p(y|θ) = ln ∑ x p
(1)

Maximizing the likelihood can be difficult, instead, it appears easier to maximize its lower bound.

Introducing an arbitrary probability distribution q(x):

ln p(y|θ) = ln ∑

x q(x) p(y, x|θ) q(x)

≥ ∑ x q(x) ln p(y, x|θ) q(x) L (q, θ), (2)
where the inequality showing that the last term L (q, θ) is the lower bound of ln p(y|θ) follows from Jensen's inequality. EM's principle is then built from a two stages coordinate ascent procedure: at iteration l +1, in the first step, referred to as Expectation step (E-step), it can be proven that when setting q l+1 (x) to the posterior conditional distribution p(x|y, θ l ), the functional L (q, θ l )

is maximized. This step can be interpreted as achieving the fixed lower bound L (q, θ l ) predetermined through θ l . Interested readers may refer to [START_REF] Jordan | An Introduction to Probabilistic Graphical Models[END_REF]; [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF] for a proof based on Kullback-Leibler (KL) divergence; in the second step, referred to as Maximization step (M-step), we start by substituting q l+1 (x) = p(x|y, θ l ) into L (q, θ):

L (q l+1 , θ l+1 ) = ∑ x p(x|y, θ l ) ln p(y, x|θ l+1 ) -∑ x p(x|y, θ l ) ln p(x|y, θ l ). (3)
Since the second term is independent of θ l+1 , only the first term, the expectation of the log joint distribution with respect to the posterior conditional distribution, needs to be considered in the maximization step. Consequently, we search for

θ = arg max θ Q(θ, θ l ), (4)
where Q is defined by

Q(θ, θ l ) = ∑ x p(x|y, θ l ) ln p(y, x|θ).
(5)

This step can as well be interpreted as finding θ l+1 to increase the lower bound L (q l+1 , θ l+1 ).

Thus we iterate over the two consecutive steps so as to increase the original likelihood function gradually. This iterative procedure is pursued until a convergence criterion is satisfied -usually set to a small constant threshold of the change of θ between two consecutive iterations. Although the EM algorithm is guaranteed to increase iteratively the likelihood function, approximate computation in the E-step and/or M-step instead of exact solutions may hinder in practice/theory global convergence.

Several comments can be made by comparing the above EM formulation to most data assimilation approaches. First, we notice that the latent variables, corresponding to the internal states, are imperative and crucial in the dynamics modeling process because they reflect intrinsic mechanisms of the studied phenomenon. Second, typical DA schemes seek to estimate the posterior conditional distribution of the latent variable x. This indeed corresponds to the Expectation step in EM algorithm. Lastly, the augmented-state technique treats both the parameters and state on an equal footing as random variables, -which is totally correct in the Bayesian context -yet, state and parameters are usually of different natures. Indeed, the hidden state is generated by the dynamical model while the dynamics itself is known only up to the parameters' value. On the other hand, the EM algorithm allows to build two separated estimation procedures for the state and the parameters respectively. Yet, this is done in a consistent Bayesian framework. Here the term 'consistent' is used to describe the synchronous and compatible relationships between two objective distributions for estimating the state and the parameters respectively from a time-sequence of data.

From a statistical perspective, the augmented-state technique tries to find the state/parameters by maximizing the posterior joint conditional probability function p(x, θ|y) while the EM algorithm works on the marginalized likelihood instead. Searching for the augmented state is conceptually easy to implement but can pose practical optimization problems because the state and the parameter may differ by orders of magnitude. This generally requires the use of a proper scaling technique in the optimization. The joint optimization procedure can be decoupled into two separate optimization procedures for the state and the parameters. This decoupled strategy, forming the dual-filter approach, despite its simplicity, is however likely to produce inconsistent state and parameter estimations. Also, augmented-state technique fails in estimating the stochastic parameters (e.g. model error covariance matrix Q) [START_REF] Delsole | State and parameter estimation in stochastic dynamical models[END_REF]. This is related to the parameter identifiability problem due to the lack of correlation between the state and the parameters. Augmented state techniques have their own advantages in easily building the covariance information.

However, this is also their weakness. When the considered parameter has a spatial structure, the size of the covariance matrix increases quadratically with the size of the parameter. Whereas in the EM algorithm, the parameter covariance is not needed. Due to the complexity of the optimization procedure (e.g. associated with a low-rank covariance matrix), the augmented-state method is known to exhibit instability issues yielding model blow-up or ensemble collapse phenomenons discussed earlier. Even if the inference problem can be solved smoothly, the added number of degrees of freedom associated with the parameters may lead to a so-called over-fitted model, which is too sensitive to the observation noise.

Recent advances on integrating EM algorithm into data assimilation feature the work of [START_REF] Dreano | Estimating model-error covariances in nonlinear state-space models using Kalman smoothing and the expectation-maximization algorithm[END_REF]; [START_REF] Pulido | Stochastic parameterization identification using ensemble kalman filtering combined with maximum likelihood methods[END_REF]. Their work focused on the estimation of a low-rank model error covariance terms associated to a reduced order noisy dynamics. The form of the noise term has been set empirically. Such empirical form for the randomization of the dynamical system, which appears to be innocent at first glance, may indeed be problematic. As a matter of fact, the energy intake due to the increase of variance brought by the noise model is not compensated by any dissipation mechanism. This imbalance may lead to numerical instability and to the break of physical conservation laws as the variance may increase in an uncontrolled way. In addition, a nonphysical random forcing (i.e. with a form that does not ensue from a physical conservation law) may lead to random dynamics whose low-noise limit differs significantly from the deterministic original system [START_REF] Chapron | Large-scale flows under location uncertainty: a consistent stochastic framework[END_REF]. In other words the random dynamics converges poorly to the deterministic system for low noise.

In this work, we investigate a different objective that consists in estimating a model parameter ensuing from a rigorous random representation of geophysical fluid dynamics. We show in particular how the estimation of such parameters is indeed tractable using an iterative ensemble-based EM algorithm. To that end, we begin by a brief review on how the effect of small scale oscillations exerted on large scale dynamics is naturally introduced through a stochastic transportation principle, followed by a formal introduction of the EM algorithm in the context of DA. We later present the EnRTS smoother that is used in the state estimation stage of EM. Then we center on the ensemble-based EM algorithm that allows the physical parameters to be explicitly estimated.

A detailed comparative study of our proposed EM-based parameter estimation scheme with the EM used in [START_REF] Dreano | Estimating model-error covariances in nonlinear state-space models using Kalman smoothing and the expectation-maximization algorithm[END_REF] and the ML used in [START_REF] Pulido | Stochastic parameterization identification using ensemble kalman filtering combined with maximum likelihood methods[END_REF] is presented in section 6. In this section, we also compared our method with two closely related approaches: one by [START_REF] Delsole | State and parameter estimation in stochastic dynamical models[END_REF] where the generalized maximum likelihood approach combined with the dual-Kalman-filter is employed to estimate the state and the parameter separately, the other by [START_REF] Ueno | Bayesian estimation of the observation-error covariance matrix in ensemble-based filters[END_REF] where EM is used with a Gaussian mixture observation model represented by samples to estimate the observation error covariance matrix. We conclude with results of numerical assessments and related discussions.

Location uncertainty modeling related to scale discrepancy

Location uncertainty principle is primarily based on the decomposition of the fluid particle displacements as:

X t = X t 0 + t t 0 w(X s , s)ds + t t 0 σ(X s , s)dB s , (6)
with X t denoting the flow map at time t, w the smooth fluid particles velocity corresponding to the sought "large-scale" solution, and dB t a d-dimensional Wiener increment (d=2,3) function. The term σ(x,t)dB t represents the fast-changing fluctuation associated to turbulence; it has a much lower time-scale than w, the smooth component. This term is defined from an integral operator of the Brownian function over the fluid domain Ω:

(σ (x,t)dB t ) i = Ω ∑ j σ i j (x, y,t)dB j t (y)dy (7)

Both σ and w require to be solved or modeled. Stochastic Navier-Stokes equations can be deduced by considering a stochastic representation of the Reynolds transport theorem (RTT) [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF], through variational principle [START_REF] Holm | Variational principles for stochastic fluid dynamics[END_REF] or homogenization techniques [START_REF] Cotter | Stochastic partial differential equations as a diffusive limit of deterministic Lagrangian multi-time dynamics[END_REF]. The stochastic representation of RTT introduces meaningful terms for turbulence modeling. It includes (i) a diffusion term that figures the mixing of the large-scale component by the small-scale random field; (ii) a modified large-scale advection describing the statistical action of the turbulence inhomogeneity on the large-scale drift; (iii) a multiplicative random term depicting the transport of the large-scale component by the small-scale random component. The variance increase brought by this latter term is in equilibrium with the amount of diffusion associated to the first term (Resseguier et al., 2017a). The inhomogeneous diffusion term can be seen as a generalization of classical eddy-viscosity subgrid models. The modified advection term corresponds to a noise induced velocity term, which drives the fluid particles from the turbulent region of high variance toward more quiescent areas with minimum variance. And finally, the multiplicative advection random process describes well the shape of the nonGaussian probability distribution function associated to turbulent flows, with fat tail emerging from turbulence intermittency [START_REF] Kraichnan | Small-scale structure of a scalar field convected by turbulence[END_REF][START_REF] Kraichnan | Anomalous scaling of a randomly advected passive scalar[END_REF][START_REF] Majda | Simplified models for turbulent diffusion:theory, numerical modelling, and physical phenomena[END_REF][START_REF] Sura | Multiplicative noise and non-gaussianity: A paradigm for atmospheric regimes?[END_REF].

The modeling under location uncertainty has been used to derive stochastic expression of Geophysical flows (Resseguier et al., 2017a,b,c), efficient reduced order flow models [START_REF] Chapron | Large-scale flows under location uncertainty: a consistent stochastic framework[END_REF]Resseguier et al., 2017d), and large eddies simulation models for various types of flows [START_REF] Chandramouli | Coarse large-eddy simulations in a transitional wake flow with flow models under location uncertainty[END_REF][START_REF] Harouna | Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling[END_REF].

In this article we focused on the shallow water model under location uncertainty. In the following we give the 1D stochastic shallow water directly in terms of the depth h, the free surface, η = h + b with b the bottom height and the velocity u. For a complete description of the model derivation, please refer to [START_REF] Yang | High-resolution data assimilation through stochastic subgrid tensor and parameter estimation from 4DEnVar[END_REF].

d t h + ∂ x h(u - 1 2 ∂ x a xx u ) - 1 2 ∂ x (a xx ∂ x h) dt + ∂ x h(σ dB t ) x = 0, (8a)
d t (hu) + ∂ x (hu u + 1 2 gη 2 )- 1 2 ∂ x a xx ∂ x (hu) + a xx ∂ x h∂ x u dt + u∂ x h(σ dB t ) x = 0, (8b)
The system (8) is constituted essentially by stochastic partial differential equations. Besides the fast-fluctuation part associated with the Wiener increment dB t , extra terms associated with the variance tensor, a, are involved. This quantity is directly related to the variance of the small-scale component as it is defined from the diagonal of the uncertainty covariance:

a(x,t)δ (t -t ′ )dt = E (σ(x,t)dB t ) σ(x,t ′ )dB t ′ T . (9)
3. An EM algorithm for model parameter estimation

The data assimilation problem can be formulated under a state-space model framework. A general continuous evolution model can be described by a discrete Markovian nonlinear state mapping operator:

x k = ϕ k (x k-1 , a) + ξ k , ( 10
)
where k is the time index. The unknown initial condition is modeled by:

x 0 = x f 0 + ω k , (11)
where the upperscript ' f ' refers to the forecast indicating the prior knowledge about the initial condition, combined with the observation equation:

y k = H(x k ) + ǫ k , (12)
which links the observation y ∈ R m and state x ∈ R n through a nonlinear operator H. The observations are assumed mutually conditionally independent. The noises may depend on the state but are assumed Gaussian when conditioned on the state: ξ |x ∝ N (0, Q), ω |x ∝ N (0, P 0 ) and

ǫ |x ∝ N (0, R).
For the remainder of the paper, we denote the length of observational time sequences by capital K and the parameter set (a, x f 0 , Q, P 0 , R) by θ.

Remarks on systems (10), ( 11) and ( 12):

• Note that even though any parameters contributing to ϕ k and H could also be identified by EM algorithm, in this paper we consider that all sources of errors in ϕ k are attributed to the variance tensor a. We will nevertheless show (section c) how to estimate the observation covariance R.

• We still assume unimodal Gaussian errors because such assumption leads to simple and efficient estimation algorithms. To represent arbitrary non-Gaussian error distributions, Gaussian mixture model can be employed in terms of either system 6,7 or 8 [START_REF] Sondergaard | Data Assimilation with Gaussian Mixture Models Using the Dynamically Orthogonal Field Equations. Part I: Theory and Scheme[END_REF]. They are not harder to deal with, although additional approximative estimation methods must be introduced to estimate the mixing coefficients.

To find the optimal θ, the EM algorithm proceeds as follows.

In the E-step, a state estimation problem is performed, as suggested in section 1, by computing the posterior conditional distribution p(x 0:K |y 1:K , θ l ). Consequently, in the M-step, we look for θ l+1 maximizing Q(θ l+1 , θ l ).

Since we are dealing with time-sequential data, we rewrite its definition as:

Q(θ l+1 , θ l ) = p(x 0:K |y 1:K , θ l ) ln p(y 1:K , x 0:K |θ)dx. (13)
A classical Bayesian manipulation of ln p(y 1:K , x 0:K |θ), with the Markov property of the dynamics and the conditional independence of the observations with respect to the state leads to:

ln p(y 1:K , x 0:K |θ l+1 ) = ln p(x 0 ) + K ∑ k=1 ln p(x k |x k-1 , θ l+1 ) + K ∑ k=1 ln p(y k |x k , θ l+1 ). (14)
Substituting this formula into equation ( 13), yields

Q = p(x 0 |y 1:K , θ l ) ln p(x 0 )dx 0 Q 1 + K ∑ k=1 p(x k , x k-1 |y 1:K , θ l ) ln p(x k |x k-1 , θ l+1 )dx k dx k+1 Q 2 + K ∑ k=0 p(x k |y 1:K , θ l ) ln p(y k |x k , θ l+1 )dx k Q 3 . ( 15
)
Introducing the state-space model equations ( 10), ( 11) and ( 12), we have:

Q 1 = - 1 2 ln |P 0 | -E x 0 |y 1:K ,θ l { 1 2 x 0 -x f 0 2 P 0 } (16) Q 2 = - K -1 2 ln |Q| -E x k,k-1 |y 1:K ,θ l { 1 2 K ∑ k=1 x k -ϕ k (x k-1 , a l+1 ) 2 Q } (17) Q 3 = - K 2 ln |R| -E x k |y 1:K ,θ l { 1 2 K ∑ k=0 y k -H(x k ) 2 R } (18)
where E X ( f ) denotes the expectation of f with respect to the law of X and we omit the constant terms. We noted the Mahalanobis norm

f 2 A = ( f , A -1 f ) in which A -1
is an inverse covariance matrix and f , g is the L 2 inner product. Note that closed form solutions can be found for parameters such as x f 0 , P 0 , R by setting the derivative of Q with respect to corresponding parameters to zero. Readers can refer to [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF] for solutions with linear models or to Dreano et al.

(2017) when ensemble approximations of non-linear models are introduced. Our approach differs because firstly, we want to obtain directly a point estimation of the physical parameter a, and not an estimation through some compensating error covariances; secondly, our approach can yield the re-estimation for the model error covariance Q, similarly to the standard EM approaches, but also extends the scenario to the case where Q is related to a. Due to the (usual) nonlinear nature of the flow map, no closed form solution exists. The standard procedure for maximizing the cost function with respect to parameters of nonlinear autoregressive models given one dimensional time-series is provided in [START_REF] Nelson | Nonlinear Estimation and Modeling of Noisy Time-series by Dual Kalman Filtering Methods[END_REF] and [START_REF] Haykin | Kalman Filtering and Neural Networks[END_REF]. Here, we reformulate the scheme to deal with parameters involved in complex nonlinear dynamic model of large spatial resolution and we devise an iterative formulation relying on an ensemble approximation of the parameter-related tangent linear model. For the rest of this section, we will exhibit simpler forms of Q 1 , Q 2 and Q 3 by applying the expectation operation inside the sum operator. The principle is to obtain forms only function of the smoothed state, its covariance and the tangent linear dynamical model, which can be effectively approximated through ensemble methods. To start with, let us distinguish at iteration l + 1 between the parameters to-be-sought (a l+1 , Q l+1 ) and the prior parameters (a l , Q l ).

After identifying that Q 1 is independent of (a l+1 , Q l+1 ), we can discard Q 1 and focus on Q 2 and Q 3 . By applying the expectation operator inside the sum operator1 , and using properties of the matrix trace (tr) (namely -for square matrices A and B -linearity: tr(aA

+ bB) = atr(A) + btr(B), tr(A T ) = tr(A), and tr(X T Y ) = tr(XY T ) = tr(Y T X ) = tr(Y X T ),tr(X T X ) = X 2
where X and Y

are n × m matrices), Q 3 can be expanded to first order as:

Q 3 = - K 2 ln|R|- 1 2 K ∑ k=1 y k -H(x k|K ) 2 R + tr (∂ x H ) T R -1 ∂ x H P k|K , (19)
where capital K after the conditioning delimiter '|' is a shortcut notation for y 1:K , and ∂ x H denotes the tangent linear model of H around x k|K . The terms x k|K and P k|K denote the conditional mean and covariance of x k respectively:

x k|K = E x k |y 1:K (x k ), (20)
P k|K = E x k |y 1:K (x k x T k ) -x k|K x T k|K . (21)
It can be easily verified that the terms x k|K and P k|K are independent of the parameter to-besought a l+1 , therefore Q 3 does not depend on a l+1 (at first order) and can be discarded in the minimization. Similarly, applying the expectation operation inside the sum operator of Q 2 gives,

Q 2 = - K -1 2 ln|Q| - 1 2 K ∑ k=1 x k|K -x - k|K 2 Q + tr Q -1 (P k|K + P - k|K -P † k|K -(P † k|K ) T ) (22)
where we have defined:

x - k = ϕ k (x k-1 , a l+1 ), (23)
x

- k|K = E x k,k-1 |y 0:K (ϕ k (x k-1 , a l+1 )), ≈ ϕ k (x k-1|K , a l+1 ), (24)
P - k|K = E x k,k-1 |y 0:K (x - k -x - k|K )(x - k -x - k|K ) T = E x k,k-1 |y 0:K (x - k x - k T ) -x - k|K (x - k|K ) T , (25)
P † k|K = E x k,k-1 |y 0:K (x - k -x - k|K )(x k -x k|K ) T = E x k,k-1 |y 0:K (x - k x T k ) -x - k|K x T k|K . ( 26
)
where x - k|K and P - k|K are the conditional mean and covariance of x - k , respectively. From now on, to simplify the notation, we will drop the dependences on y and x k,k-1 in the expectation. At first order, the term P - k|K can be rewritten in terms of the covariance P k-1|K :

P - k|K =E x {[ϕ k (x k-1 , a l+1 ) -ϕ k (x k-1|K , a l+1 )] [ϕ k (x k-1 , a l+1 ) -ϕ k (x k-1|K , a l+1 )] T } =E x {∂ x ϕ k (x k-1|K , a l+1 ) • (x k-1 -x k-1|K ) (x k-1 -x k-1|K ) T • ∂ x ϕ T k (x k-1|K , a l+1 )} =∂ x ϕ k (x k-1|K , a l+1 ) • E x {(x k-1 -x k-1|K ) (x k-1 -x k-1|K ) T } • ∂ x ϕ T k (x k-1|K , a l+1 ) =∂ x ϕ k (x k-1|K , a l+1 )P k-1|K ∂ x ϕ T k (x k-1|K , a l+1 ). (27)
In the second line, we introduced the tangent linear operator of the model ϕ k around x k-1|K and its adjoint. These operators are denoted by ∂ x ϕ k (x k-1|K , a) and ∂ x ϕ T k (x k-1|K , a), respectively. The same manipulations can be applied to the conditional cross-covariance of x k and x - k denoted by P † k|K :

P † k|K = E x {[ϕ k (x k-1 , a l+1 ) -ϕ k (x k-1|K , a l+1 )](x k -x k|K ) T } = E x {∂ x ϕ k (x k-1|K , a l+1 ) • (x k-1 -x k-1|K )(x k -x k|K ) T } = ∂ x ϕ k (x k-1|K , a l+1 ) • E x {(x k-1 -x k-1|K )(x k -x k|K ) T } (28)
After these manipulations, we see that, when searching for a l+1 , it is necessary to consider the covariance terms P - k-1|K and P † k|K that are indeed both functions of a l+1 . Indeed, this variance tensor parameter has a direct impact on the tangent linear operator ∂ x ϕ k (x k-1|K , a l+1 ). Note that if the dynamic model is linear, a closed form for a l+1 can be obtained by setting

∇ a Q 2 = 0. Looking for Q l+1 , the model error covariance, amounts to finding Q l+1 solution of ∇ Q Q 2 = 0. A closed
form solution generally exists regardless of the nature of the dynamics. However, as suggested in our stochastic formulation (8), estimating a and Q are not entirely independent of each other. So in section 5, we show a possible solution of this coupled problem within an ensemble assimilation framework for the physical parameter a as well as the model error covariance Q.

EnKF and Ensemble RTS Smoother in E-step

In the former section, we showed that in the E-step, a smoothing posterior distribution of the state conditioned on the whole data span and on the prior parameters must be found. The method for computing such smoothing state is quite versatile. Methods referred to as the Rauch-Tung-Striebel smoother (RTS), the two-filters smoothing or the joint smoothing enable to perform such estimation. Interested readers may refer to [START_REF] Cosme | Smoothing Problems in a Bayesian Framework and Their Linear Gaussian Solutions[END_REF] and the references therein for a complete review on those methods. In this study, we choose to rely on an ensemble version of the RTS smoother (EnRTS) similar to the one proposed in [START_REF] Raanes | On the ensemble Rauch-Tung-Striebel smoother and its equivalence to the ensemble Kalman smoother[END_REF], as this smoother can be performed efficiently and is compatible with the EnKF methods.

a. EnKF

The central idea underlying EnKF consists in approximating the estimation covariance matrix by a few samples. Therefore the intractable covariance evolution problem is transformed into a tractable sampling evolution followed by an empirical approximation. By denoting E k the ensemble state made of N sampled states at time k:

E k = (x 1,k , . . ., x l,k , . . ., x N,k ) ∈ R n×N , (29)
the covariance can be approximated through the sampling covariance matrix P e :

P k ≈ P e = 1 N -1 (E k -E k )(E k -E k ) T , (30)
where the columns of E k consist of the ensemble mean, (also referred to as the forecast ensemble):

x k = 1 N ∑ N l=1 x l,k . It is also convenient to define the ensemble anomaly matrix A k = 1 √ N-1 (E k -Ēk ) that leads to P k = A k A k T .
E f k = ϕ k (E a k-1 ). (31)
The forecast covariance matrix is implicitly inferred from the forecast ensemble

P f k = A f k A f k T .
The following analysis/update stage aims at estimating the correction of state given new observation y k at time t k . In stochastic EnKF approach, an observation ensemble Y k is built by perturbing the observation y k with additive noise η l ∝ N (0, R). The update scheme reads:

E a k = E f k + P f k ∂ T x H(∂ x H P f k ∂ T x H + R) -1 (Y k -H(E f k )). (32)
In practice, P, is never computed explicitly. An efficient implementation is formulated in terms of the ensemble anomaly matrix A:

E a k = E f k + A f k (∂ x H A f k ) T S -1 k D k , (33)
where we define given by:

S k = ∂ x H A f k (∂ x H A f k ) T +R and D k = Y k -H(E f k ).
E k|K = E a k + J k (E k+1|K -ϕ k (E a k )), (34)
with E K|K = E a K . In eqs (34), J k reads (Eq. 17 in Raanes ( 2016)):

J k = A a k (A f k+1 ) + , (35)
where the subscript ′ + ′ stands for Moore-Penrose inverse that can be computed through the singular value decomposition (SVD) given that the ensemble approximated covariance is of reducedrank. This formulation, therefore can be applied directly to the ensemble anomaly matrix as long as the forecast/analysis ensembles supplied by EnKF are stored in memory. Note that the ensemble Kalman smoother (EnKS) is also a viable choice to produce the smoothing state. Interested readers may refer to [START_REF] Raanes | On the ensemble Rauch-Tung-Striebel smoother and its equivalence to the ensemble Kalman smoother[END_REF] for a comparative study of EnRTS smoother and EnKS.

c. Localization in EnKF and EnRTS smoother

The standard localization schemes must be applied to the EnKF scheme to handle the unavoidable fictitious long-range correlations emerging in low-rank sampled correlation matrices of high dimensional system. In this paper, we extend the Schur product-based covariance localization scheme to EnRTS smoother. In this approach, we use an altered form of the global ensemble sampling matrix to compute J k :

J k = A a k (A f k+1 ) -1 , ( 36
)
where A is a localized version of A such that A A T = L⊙(AA T ). Therefore, the EnRTS smoother update scheme with a localized covariance remains exactly of the same form. It consists just to replace A by its localized counterpart A .

An ensemble formulation of M-step

In the previous section, we introduced EnKF and EnRTS smoother to solve the E-step at the (l + 1)th EM iteration. The ensemble E-step yields the smoothing ensemble state E k|K . A low rank approximation of the smoothing covariance P k|K as well as P - k|K and P † k|K in equation ( 22) are provided by the smoothing ensemble anomaly matrix A k|K :

x k|K = E k|K , (37)
P k|K = A k|K A T k|K , (38)
P - k|K = A - k|K (a l+1 )(A - k|K (a l+1 )) T , (39)
P † k-1,k|K = A - k|K (a l+1 )A T k|K , (40)
where A - k|K (a l+1 ) denotes the evolution of the smoothing ensemble anomaly matrix A k-1|K , expressed as the difference between the smoothing ensemble driven by the nonlinear model and its mean:

A - k|K (a l+1 ) ∂ x ϕ k (x k-1|K )A k-1|K ≈ 1 √ N -1 ϕ k (E k-1|K ) -ϕ k (E k-1|K ) . ( 41
)
Substituting the above ensemble approximations into Eq.( 22), and defining the cost function, J M = -Q 2 , allows us to transform the maximization problem into a minimization one with:

J M (a l+1 , Q l+1 ) = K -1 2 ln|Q l+1 | + 1 2 K ∑ k=1 x - k|K (a l+1 ) -x k|K 2 Q l+1 + tr (A - k|K (a l+1 ) -A k|K ) T Q -1 l+1 (A - k|K (a l+1 ) -A k|K ) . (42)
One remark on equation ( 42):

• In our stochastic dynamics (8), the random part depends on a therefore minimizing directly J M against a requires computing ∂ Q/∂ a. This can be computationally demanding for high dimensional state and parameter. As a result, we employ a decoupled estimation procedure in which J M is minimized with respect to a using the current estimate Q (that can be factored and hence plays only an indirect role in the minimization through the ensemble driven by the random dynamics and the Q-Wiener process) (section a), then the minimization of J M is led with respect to Q with J M evaluated using the current estimate a l+1 (section b).

a. Physical stochastic parameter estimation

To facilitate the ensemble formulation, the cost function (42) maybe reformulated (at first order) in terms of the increment δ a l+1 a l+1a l :

J M (δ a l+1 ) = 1 2 K ∑ k=1 N ∑ j=1 { 1 √ N -1 [ϕ k (E k-1|K , a l ) -ϕ k (E k-1|K , a l ) + ∂ a ϕ k (E k-1|K , a l ) -∂ a ϕ k (E k-1|K , a l ) δ a l+1 ] -A k|K } j 2 Q + ϕ k (x k-1|K , a l ) -x k|K + ∂ a ϕ k (x k-1|K , a l )δ a l+1 2 Q , ( 43
)
where the first term on the right hand side of eq.( 42) is discarded as it depends only on Q; and the term ∂ a ϕ k (x k-1|K , a l ) denotes the tangent linear model of ϕ k with respect to a:

lim α→0 ϕ k (x k , a + αda) -ϕ k (x k , a) α = ∂ a ϕ k (x k , a)da. (44)
In the same spirit as ensemble methods, where the increment is sought in the range of the ensemble anomaly matrix, we introduce here an additional ensemble induced by the perturbation of a. Such perturbation gives rise to an ensemble E a ∈ R p×M where p is the size of the parameter vector and M is the size of the parameter ensemble. By doing so, we fix the search space in the span of the column vectors of A a :

δ a l+1 = A a Υ = 1 √ N -1 (E a -Ēa )Υ. (45)
The cost function in terms of Υ now reads:

J M (Υ) = 1 2 K ∑ k=1 N ∑ j=1 { 1 √ N -1 [ϕ k (E k-1|K , a l ) -ϕ k (E k-1|K , a l )] -A k|K + 1 N -1 ∂ a ϕ k (E k-1|K , a l ) -∂ a ϕ k (E k-1|K , a l ) A a Υ} j 2 Q + ϕ k (x k-1|K , a l ) -x k|K + 1 √ N -1 ∂ a ϕ k (x k-1|K , a l )A a Υ
The terms ∂ a ϕ k (x k-1|K , a l )A a and ∂ a ϕ k (E k-1|K , a l )A a denote the evolution of the parameter ensemble perturbation through the tangent linear model with respect to the smoothed state x k-1|K .

These two terms can be approximated in a similar way as (41). The tricky part is to compute

∂ a ϕ k (E k-1|K , a l )A a , since the tangent linear term ∂ a ϕ k (E k-1|K , a l
) is a function of the state en- semble. Therefore evaluating this term consists in perturbing each state sample with the parameter ensemble:

∂ a ϕ k (E k-1|K , a l )(E a -Ēa ) ≈ [ϕ k (x 1 k-1|K , E a ) -ϕ k (x 1 k-1|K , E a ), • • • , ϕ k (x N k-1|K , E a ) -ϕ k (x N k-1|K , E a )]. (47)
Here ϕ k (x j k-1|K , E a ) ∈ R n×M , as ϕ k is applied to x j k-1|K and each column of E a . The above ensemble approximation yields an ensemble of ensembles with N × M state samples in total. Although the computational time needed for simulating the whole sample space trajectory is potentially larger, the ensemble implementation itself is simple and its high parallelism potential has not been weaken. Parallelization techniques can thus be introduced to reduce the computing time. Besides, memory space can be saved because we do not need to store all time series of the ensemble of N × M members as long as their contributions to the corresponding cost function term is computed and accumulated. Finally, all terms in (46) can be derived from the smoothing state ensemble E k|K and the parameter ensemble E a . The ensemble approximation transforms the original nonlinear least-square optimization problem to a (Gauss-Newton like) series of incremental linear systems (in terms of the anomaly ensemble coefficients ϒ) that can be solved efficiently using gradientbased algorithm such as the conjugate-gradients technique.

b. Model error covariance estimation

In the general case, the model error covariance term Q can be optimally estimated by simply setting to zero the derivative of formula Q 2 (equation 22) with respect to Q. At the lth M-step in EM iteration, the re-estimated Q reads:

Q l+1 = 1 K -1 K ∑ k=1 [E x (x k x T k ) + E x (x - k x - k T ) -E x (x - k x T k ) -E x (x k (x - k ) T )]. (48)
Using a similar ensemble formulations, we derive the following scheme in evaluating Q:

Q l+1 = 1 K -1 K ∑ k=1 [A - k|K (a l+1 ) -A k|K ][A - k|K (a l+1 ) -A k|K ] T + [x - k|K (a l+1 ) -x k|K ][x - k|K (a l+1 ) -x k|K ] T (49)
Note that A - k|K and x - k|K have to be recomputed using the re-estimated parameters a l+1 at the lth EM iteration. We summarize the steps of the proposed ensemble-based EM algorithm below.

Let us outline that the noise resampling procedure requires some cares. The noise sampling can be performed through two additive noises with diffusion matrices associated to the square roots A - k|K -A k|K and x - k|Kx k|K respectively. However these noises must be independent (in order to get the right covariance). To do so, the eigenspace of the first covariance matrix can be split into two complementary subspaces to sample the two independent noises. The first one is used to sample the first noise while a projection on the complementary subspace is used to sample the second one. Note that the first square root is likely to have a greater number of non zero singular values as it relies on the ensemble of ensembles -as opposed to the second one which is associated to the ensemble of trajectories.

Algorithm 1 Ensemble-based EM algorithm 1: Given observation sequence y from 0 to K, prior initial state x 0 l and prior parameters a l and Q l at l = 0: 2: while convergence criterion is not met do 3:

Expectation step:

4:

for k = 0 : K do 5:

Forecasting state ensemble update (31) with parameter a l and additive error resampled from Q l using SVD 6:

Filtering state ensemble update (33) 7:

Smoothing state ensemble update (34) 8:

Maximization step:

9:

Obtain a l+1 through ( 46) and recover δ a l+1 using ( 45) 10:

Update model error covariance Q l+1 using (49)

11:

Set l . = l + 1, a l . = a l+1 and Q l . = Q l+1 then return to step 2.

Discussions

In this section, we present a detailed discussion on the comparison between the EM-based parameter estimation scheme we proposed and other prominent methods of the literature. formulation corresponds only to the second line of ( 49) and can be interpreted as a first-order approximation of our complete formulation (49). Besides, our formulation for Q l+1 requires the re-estimation of the physical parameters, a, beforehand. For complex error terms, such as parametric error or multiplicative error, their approach may alter the error characterization by estimating an additive Gaussian error covariance compensating the error. [START_REF] Pulido | Stochastic parameterization identification using ensemble kalman filtering combined with maximum likelihood methods[END_REF] used the same EM algorithm as in [START_REF] Dreano | Estimating model-error covariances in nonlinear state-space models using Kalman smoothing and the expectation-maximization algorithm[END_REF] for evaluating the model error covariance. They restated the limitation of EM algorithm when dealing with physical parameters, as closed analytical forms for estimating the physical parameters generally do not exist. So, as an alternative mean to estimate the physical parameters with an EM technique, [START_REF] Pulido | Stochastic parameterization identification using ensemble kalman filtering combined with maximum likelihood methods[END_REF] proposed to rely on an augmented state technique, in which the state and the parameter are jointly estimated during the filtering/smoothing phase of the E-step, whereas the augmented model error covariance is estimated in the M-step. This approach is shown to perform well in terms of the one-scale Lorenz-96 system with a specific polynomial noise term whose coefficients are defined from a Gaussian process. In this way, the mean and the variance of the corresponding component, acting as the deterministic parameter and the stochastic parameter respectively, can be both estimated using the E-step and the M-step of this augmented-state based EM technique. Nevertheless the particular form of the model limits any generalizations of their approach. Let us note that this approach has the same restriction in terms of the parameters dimension as the augmented state technique and has been used

to estimate only parameters of very small dimension. Contrary to these strategies, our formulation is able to identify not only the unknown parameters, associated to any forms of the dynamics, but also the model error covariance. We will show in the experimental section that it is also beneficial to use a full EM algorithm to estimate physical parameters of much higher dimension.

Apart from the EM approach based on augmented state technique, [START_REF] Pulido | Stochastic parameterization identification using ensemble kalman filtering combined with maximum likelihood methods[END_REF] also resorted to a maximum likelihood approach that seeks to maximize the likelihood function of the observations given the parameters, p(y 1:K |θ). Their formulation relying on the Markov property of time-sequential data and Bayes rule led to a factorized form:

ln p(y 1:K |θ) = ln K ∏ k=1 p(y k |y 1:k-1 , θ). (50)
For comparison purpose, we recall here the marginalized full joint distribution over the state, x, used in EM:

ln p(y 1:K |θ) = ln ∑ x p(y 1:K , x 0:K |θ). (51)
One advantage of ML approach is that only the filter is required to get a solution. By introducing the state space formulation (Eqs. ( 10), ( 11), ( 12)) to p(y k |y 1:k-1 , θ), the Kalman filter solution is exactly recovered for linear models and additive Gaussian errors. For general nonlinear and/or non-Gaussian cases, Monte-Carlo methods must be employed to approximate this distribution.

The challenge in terms of the formulation ( 50) is that, in the case of non-additive model error, the gradient of the likelihood w.r.t θ can be hard to evaluate. To that end, [START_REF] Pulido | Stochastic parameterization identification using ensemble kalman filtering combined with maximum likelihood methods[END_REF] proposed a derivative-free optimization method to maximize ln p(y 1:K |θ). However, despite their practical interest, derivative-free methods are generally not among the most efficient optimization technique. The original likelihood form (51) in EM is at the opposite hard to manipulate. So a lower bound that is iteratively increased towards the original likelihood ( 51) is defined instead.

A great advantage of EM, compared to the ML approach, is that in order to maximize Q (Eq. ( 13)) in the M-step, only the complete joint distribution of the state and observation conditioned on the parameters is needed since the posterior conditional distribution term is known from the E-step. Indeed, in many scenarios, deriving the gradient of the complete joint distribution w.r.t the parameters is much easier than deriving the gradient of (50).

b. Differences compared to Dual-Kalman-Filter approach

Delsole and Yang (2010) advocated the use of generalized maximum likelihood approach to estimate the model state and the stochastic parameters. This approach is still based on the maximization of the posterior joint conditional distribution p(x, θ|y). But instead of using an augmented state technique (which can cause ensemble collapsing), the optimization process is decoupled into two alternate procedures, in which the parameters and the state are alternately optimized. The procedure can be implemented through two separated Kalman-like filters, namely a state-filter and a parameter-filter. This approach is an online method. Such a dual-filter technique can be readily extended to dual-ensemble-filter approach to handle non-linear dynamical operators.

Compared to this technique, our (offline) EM based approach maximizes the marginalized complete joint distribution through an alternate optimization. It provides in addition a coupled estimation of the smoothed state x 0:K and of the parameter θ given the set of observations from 0 to K.

Another major difference concerns the treatment of the error covariance term (noted P f in Delsole and [START_REF] Delsole | State and parameter estimation in stochastic dynamical models[END_REF] and Q in our manuscript). Note that P f in a filter context corresponds to the composition of the background error covariance's evolution P 0 and of the model error covariance Q. The filtering formulations in [START_REF] Delsole | State and parameter estimation in stochastic dynamical models[END_REF] is general, so that it can be applied to deterministic or stochastic parameter estimation problems. However, the authors focused on a particular case where the error covariance is a function of the stochastic parameter. The formula for updating the parameters in the parameter-filter stage requires the derivative w.r.t. the parameters of the covariance term P f . Due to the complexity of computing those derivatives, the authors proposed an evaluation technique relying on an additional parameter-induced ensemble, in the same way as we did. Formally in our stochastic shallow water model, the variance tensor, a plays a role both in the drift part and in the random part. The random term is a non-Gaussian error term, that can be hardly directly estimated. So we first estimate the variance tensor, a, as a deterministic physical parameter accompanied with a centered Gaussian error of covariance Q; and then, as depicted in algorithm 1, we consider a decoupled estimation approach by minimizing the cost function Q 2 (Eq.( 22)) w.r.t. a and Q separately in the M-step.

Apart from the fundamental differences discussed above, our formulation shares thus similarities with [START_REF] Delsole | State and parameter estimation in stochastic dynamical models[END_REF] in the way of approximating the gradient of the state w.r.t. the parameters using an ensemble of ensembles. we show here the estimation of R that ensues from our approach and the maximization of the Q 3 term (Eq. ( 18)). By solving ∇ R Q 3 = 0 with the ensemble approximation, we get,

R l+1 = 1 K K ∑ k=1 [ y k -H(x k|K ) 2 + ∂ x H A k|K (∂ x H A k|K ) T ], (52)
which can be computed using the smoothed ensemble state. Comparing this expression with formula (35) in [START_REF] Ueno | Iterative algorithm for maximum-likelihood estimation of the observation-error covariance matrix for ensemble-based filters[END_REF], we observe that the EM approach in [START_REF] Ueno | Iterative algorithm for maximum-likelihood estimation of the observation-error covariance matrix for ensemble-based filters[END_REF] does not involve any state inference problem in the E-step as their updating formula ( 35) corresponds only to an innovation term (i.e. the first term of ( 48)) related to the filtered state. Their approach requires thus only an EnKF procedure before proceeding to the EM algorothm. Due to this trait, their method can be implemented online. However, it cannot be used to estimate the state and the model error covariance Q. This approach has been extended in [START_REF] Ueno | Bayesian estimation of the observation-error covariance matrix in ensemble-based filters[END_REF] to incorporate a prior model on R and to proceed to a Bayesian a posteriori estimation.

EnVar-EM with augmented state

An alternative approach consists in combining the EM algorithm with EnVar [START_REF] Yang | High-resolution data assimilation through stochastic subgrid tensor and parameter estimation from 4DEnVar[END_REF] to estimate the parameters by an augmented state technique. The iterative procedure in EM is in the same spirit as the iterative resampling strategy used in the EnVar technique proposed in [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF]. This technique, built from a Gauss-Newton nested-loop variational framework, contains an update step of both the prior/posterior ensemble and the associated ensemble anomaly matrix. The cost function (in non-incremental form) in terms of the augmented state vector

s k = [x k , a] T reads: J(s 0 ) = 1 2 s 0 2 B f s + 1 2 K ∑ k=0 H(ϕ k (s 0 )) -y k 2 R , (53)
where the augmented error covariance matrix B s is given by: 

  B xx B T ax B ax B aa    . (54)
Note that the EnVar method relies on a strong dynamical constraint, thus no explicit model error term is considered here. However, between two consecutive outer loops, the covariance matrix B aa along with B xx are re-estimated following the joint update of the augmented state/parameter vector. This strategy, which is similar to the one used in [START_REF] Pulido | Stochastic parameterization identification using ensemble kalman filtering combined with maximum likelihood methods[END_REF], can be interpreted as the Maximization step in the EM algorithm. For the experimental assessments, presented in the next section, we applied both the Aug-EnVar-EM method and the proposed EnKF/EnRTS-EM method to the same dynamical models and data sets. The goal here is to evaluate numerically the difference of behavior between the augmented state techniques and the EM algorithm in retrieving unknown state/parameter fields.

Experiments and results

a. Experimental settings

We have devised two different scenarios using a similar experimental setup as in [START_REF] Yang | High-resolution data assimilation through stochastic subgrid tensor and parameter estimation from 4DEnVar[END_REF]. In all cases, we use the 1D stochastic shallow water system (8). The 1D domain has length L x = 6000 km with the initial surface height h(x, 0) given by:

h(x, 0) = H 0 - fU 0 x g + Aξ (x), (55)
where H 0 = 5000m, f = 1.03 × 10 -4 s -1 is the Coriolis parameter, U 0 = 40m/s and g is the gravity acceleration. An additional noise is considered for the initial condition with A = 10 as its amplitude. The Gaussian random field ξ is simulated through the spectral method discussed in Evensen ( 2003) (with de-correlation length equals to 20%L x ). The initial velocity field is inferred from the geostrophic relation.

The first case corresponds to a typical state and parameters estimation problem, for which the true stochastic parameters distribution is known. The background state is run on the same resolution (26 grid points) as the true state simulation. The objective of this case is to evaluate quantitively the performance of different approaches in the recovering of the parameters as well as the state.

Case II corresponds to the same configurations as the one used in [START_REF] Yang | High-resolution data assimilation through stochastic subgrid tensor and parameter estimation from 4DEnVar[END_REF] where no true model of the parameters is prescribed. Instead, the resolution-induced error is intrinsically defined and parameterized according to the stochastic dynamics. In that case, the quality of the estimated parameters can only be implicitly evaluated through the accuracy of the state recovery.

To that end, the coarse-resolution background state with the same initial condition as in case I is simulated with the stochastic model on 26 grid points, while the true state used to construct the fine-resolution observation is obtained by integrating the standard shallow water equation with the initial condition (55) on a grid of finer resolution (401 grid points). Meanwhile, the time step ∆t is set to 15s for the true state model and to 30s for the coarse-resolution model to satisfy the different CFL conditions.

In both cases, the synthetic observations are extracted from the true state every 150s by adding an i.i.d Gaussian noise to the true state. The magnitude of the noise level is proportional to the background errors.

The EnKF and the EnRTS smoother are performed to produce the optimal state estimation given the observations. The parameter estimation resulting from the EM algorithm is compared with the Aug-EnVar-EM approach implementing an augmented state technique. In all ensemble-related computations, we have used 32 members combined with a covariance localization technique based on a Shur product with the Gaspari-Cohn kernel [START_REF] Gaspari | Construction of correlation functions in two and three dimensions[END_REF] (see section 4.3).

b. Case I: EM numerical assessment with known true parameter

In this first case, we try to recover both the noisy background state and the unknown parameter field. Because the true parameter field is invariant in time, a single assimilation window containing all 5 observations is employed for the augmented EnVar approach. As the EnKF/EnRTS Smoother is sequentially implemented, the 5 observations are successively assimilated. The analysis is started at time t = 600s to maintain the balance of the state variables in the ensemble.

Depending on our EM algorithm formulation, whether or not the M-step can find the optimal parameter hinges on the quality of the posterior state given by the E-step. In figure 1, we plot RMSE curves of different posterior states (h and u for the top and bottom figures respectively) after consecutive EnKF/EnRTS-EM iteration. It can be noticed that immediately after the 1st Estep, the smoothed state at each observational time is already of excellent quality. The subsequent iterations improve progressively further the results in making the RMSE curves become flatter.

The flatness of the analysis RMSE curve reflects the accuracy of the estimated variance tensor parameter.

To examine the effectiveness of the maximization step of EnKF/EnRTS-EM, we begin by plotting the RMSE of the posterior parameter with respect to the reference true parameter against 10 EM iterations in figure 2. In current scenario, the RMSE result indicates that the EnKF/EnRTS-EM converges fairly fast in 3 iterations.

We also plot the exact spatial forms of the estimated posterior parameter field (black line) yielded at each EM iteration along with the current prior parameter (red line), the background parameter (green line) and the reference true parameter (blue line) in terms of EM iterations on the left column of figure 3. Note that at the 1st iteration, the prior parameter equals the background parameter.

We can verify that the true parameter field is successfully recovered by the EnKF/EnRTS-EM after the 3rd iteration although the fourth and the fifth iterations continue improving the accuracy of the posterior parameter. The right column of figure 3 depicts the evolution of the posterior parameter using Aug-EnVar-EM in terms of the outer loop iteration. The parameter evolution converges after the 3rd iteration as well. Nevertheless the sought parameter field only partly coincides with the true parameter field. In terms of the accuracy of the estimated parameter, the EnKF/EnRTS-EM method greatly outperforms the Aug-EnVar-EM.

Finally, we compare the state RMSE after the final E-step of EnKF/EnRTS-EM and Aug-EnVar-EM in figure 4. Compared to the smoothed trajectory provided by EnKF/EnRTS-EM, the initial analysis state of Aug-EnVar-EM is more accurate. However the inaccuracy of the parameter estimation interferes finally with the dynamics and makes the trajectory diverged from the reference trajectory, leading eventually to bad forecast. This is exactly an expression of the overfitting phenomenon. The analysis state yielded by the Aug-EnVar-EM corresponds too closely to the observations, but the parameters estimated from the augmented state approach, even combined with the EM strategy, leads to large generalization errors reflecting bad predictive performance of the trained model.

c. Case II: EM numerical assessment with unknown parameter

Finally we tested different parameter estimation approaches on minimizing the errors due to scale discrepancy. Since the variance tensor, a, to-be-sought is time varying, it is necessary to introduce a sliding assimilation window technique, in which the parameters are assumed invariant during one window and are then allowed to vary when shifted to the next window. For Aug-EnVar-EM approach, we employed 3 cycling windows for 6 observations. We also consider an extreme case where one window only contains one observation and the analysis time is set to the observation time. In that case, the Aug-EnVar-EM reduces to the Aug-EnKF-EM.

A similar sliding window technique is also applied to the EnKF/EnRTS-EM approach to incorporate the fast changing variance tensor. Here we implemented an adapted version of the off-line ensemble-based EM scheme shown in algorithm 1. In the original algorithm, the parameters are kept invariant during the full observation sequence from y 0 to y K . In the adapted version, we introduced K subwindows. Each subwindow started successively from the 1st to the (K -1)th observation and ended all at the final Kth observation. By doing so, we are able to update the variance tensors when a new observation occurs but also to take account of the future observations.

The analysis time of the first window starts at 2400s. In terms of the EnKF/EnRTS-EM method, the global RMSE trajectory recovered by the EM technique is clearly of better quality. The state errors at the observational time are slightly higher than the corresponding errors of Aug-EnVar/Aug-EnKF. Nevertheless the evolution of state errors between two consecutive observations is characterized by a U shape curve, which suggests that the variance tensor is reasonably recovered by the EnKF/EnRTS-EM method. The overfitting phenomenon can thus be controlled or even avoided using this marginalized parameter estimation scheme. The quality of the variance tensor estimation is clearly improved with the EnKF/EnRTS-EM method.

The distinct RMSE evolution in terms of the different approaches can also be attributed to the estimated parameter profiles plotted in figure 6. Note that the true parameter distribution is here unknown. We observed that the spatial distribution of the variance tensor for the Aug-EnKF method is characterized by high frequency fluctuations. This finding corresponds to our expectations since the smoothing procedure in both Aug-EnVar and EnKF/EnRTS-EM tends to filter out those small scale variations. The variance tensor obtained from the Aug-EnVar exhibits higher values at the boundary. This helps introducing extra forcing into the system and eventually reduces the state estimation error. The variance tensor resulting from the EnKF/EnRTS-EM, on the other hand, shows a smaller magnitude and a longer correlation, both of which enables achieving flatter RMSE curve with slowly-growing errors. This is likely due to large magnitude of the variance tensor as well as of its gradient, which may trigger numerical instability associated to a strong noise. However, the forecast quality starting at 3150s is not clearly improved by the EM approach. To improve the forecast, some kind of prediction of the parameter should probably be introduced. Learning procedures for the parameter's dynamics would be worth studying; however it is clearly a research topics by itself and we leave it for future exploration.

Conclusions

Recovering the exact parameters helps compensating for errors due to approximations made in the analytical and/or numerical constitution steps of the geophysical dynamics models. To address this issue, a new parameter estimation scheme is proposed in this paper. This scheme is based on an ensemble formulation of both the E-step and M-step of the EM algorithm. Unlike previous studies, our ensemble formulation of EM enables to identify the physical subgrid parameter of fluid dynamics model. This formulation extends the power of EM algorithm to complex dynamical model identification problems that are quite common in geophysical dynamics modeling. Under these circumstances we are more interested in finding the explicit parameter (with a possible spatial distribution) rather than the partially informative model error covariance that is not only cumbersome to compute but also far less representative of the real physics. Besides, direct parameter estimation obviates the auxiliary resampling procedure used to add random forcing during simulation. In the traditional setting, the EM algorithm is hardly applicable to such problems as closed form expressions for the parameter updating do not exist in general. A key in the proposed scheme is the use of an ensemble of ensemble in the Maximization step to compute the sensitivity matrix of the error discrepancy with respect to the physical parameters.

Boosted by this ensemble trick, the proposed scheme is able to tackle the estimation of parameters of high dimension. Our ensemble-based EM algorithm can, as well, be used to estimate the error covariance terms in a similar way to other EM algorithms but with an extended criterion. We have indeed shown that the ensemble-based covariance update scheme proposed in Dreano et al.

(2017), corresponds to an approximation of our scheme.

In this paper, we focused on the identification from high resolution observations of the model errors associated to unknown small-scale physical processes. The effect of those unresolved process on the large-scale (resolved, coarse resolution) components must be properly parameterized through a subgrid-scale stress model. Within the stochastic flow dynamics framework proposed in [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF], such subgrid-scale stress model takes the form of a diffusion term together with a corrective stochastic advection term. The subgrid model parameters, are entirely controlled by the noise variance tensor, which is thus the quantity to be inferred from the data.

Therefore our goal focused on evaluating the performance of the proposed EM algorithm in retrieving as accurately as possible such subgrid model parameters from high resolution observations. This technique has been compared to an augmented state technique, build from an EnVar method [START_REF] Yang | High-resolution data assimilation through stochastic subgrid tensor and parameter estimation from 4DEnVar[END_REF]. This latter technique supplied an accurate enough reconstruction of the state, but provided a far less accurate parameter field with the need of a further retuning by an adhoc inflation technique to avoid model blow-up and/or ensemble collapse.

The proposed EnKF/EnRTS-EM scheme has shown to overcome the shortcomings associated to augmented-state-based techniques. The twin experiment designed here consisted in two cases: the first case corresponds to a typical state and parameter estimation problem where the true stochastic parameter distribution is known, the second case corresponds to the estimation of an optimum stochastic parameterization to represent the resolution-induced error associated to the scale discrepancy between the model and the data. For the first case, the proposed EnKF/EnRTS-EM scheme has been able to retrieve with a substantial higher accuracy the parameter profile compared to the augmented state techniques. For both cases, the proposed EnKF/EnRTS-EM scheme outperformed the augmented state techniques in terms of the forecast errors. It yielded better quality of both the estimated state and parameters. Note that for the second case, the variance tensor was varying in time. To answer to this problem, a quasi-online version of the EnKF/EnRTS-EM scheme that has been implemented.

Although the proposed EnKF/EnRTS-EM scheme is very promising, concerns can be raised regarding the computational efforts associated with the introduction of an ensemble of ensemble in the M-step. However, the introduction of such mechanisms allowed us to compute an accurate ensemble based evaluation of the cost function gradient with respect to the unknown parameter.

Highly parallelized implementations as well as more efficient optimization algorithm can be further explored to reduce the computational burden and keep the same accuracy in the model error characterization.

LIST OF FIGURES

  b. EnRTS smootherFrom a statistical perspective, a filtering problem, seeks to estimate x k by maximizing the posterior conditional probability function p(x k |y 0:k ), given a sequence of past data up to the current time k. Smoothing, on the other hand, aims to estimate the posterior probability function p(x k |y 0:K ) conditioned on the whole data range from the initial time to the final time K. The filtering state ensemble E k|k of EnKF serves as input of the EnRTS smoother. The smoothing is performed recursively backward in time from the final time to the initial time. The smoothing ensemble state is

  a. Differences compared to[START_REF] Dreano | Estimating model-error covariances in nonlinear state-space models using Kalman smoothing and the expectation-maximization algorithm[END_REF] and[START_REF] Pulido | Stochastic parameterization identification using ensemble kalman filtering combined with maximum likelihood methods[END_REF][START_REF] Dreano | Estimating model-error covariances in nonlinear state-space models using Kalman smoothing and the expectation-maximization algorithm[END_REF] proposed an EM-based technique for estimating the model error covariance, Q, followed by a resampling of the model error, ξ k , from the current covariance, Q l+1 , to simulate the optimal state trajectories. In their setup the estimation of Q has a closed form formulation based on an ensemble approximation (equations 18 and 44 in[START_REF] Dreano | Estimating model-error covariances in nonlinear state-space models using Kalman smoothing and the expectation-maximization algorithm[END_REF]). This

  c. Differences compared to[START_REF] Ueno | Iterative algorithm for maximum-likelihood estimation of the observation-error covariance matrix for ensemble-based filters[END_REF] Ueno and Nakamura (2014) proposed estimating the observation error covariance, R, in a ML framework. They exploited the equivalent formulation between the sampling-based likelihood function and the Gaussian mixture model. Based on this link, the authors employed the EM algorithm to estimate iteratively R, and the latent mixture components. For comparison purpose,

Figure 5

 5 Figure5shows the analysis state RMSE curves of the different approaches. Let us examine first

  The EM technique[START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF] offers an elegant and flexible solution for estimating separately the state and the model parameter. This method has been extensively used in many research domains, such as machine learning or system identification. Studies on integrating EM algorithm into data assimilation framework have been considered only very recently.[START_REF] Dreano | Estimating model-error covariances in nonlinear state-space models using Kalman smoothing and the expectation-maximization algorithm[END_REF];[START_REF] Pulido | Stochastic parameterization identification using ensemble kalman filtering combined with maximum likelihood methods[END_REF], among others, proposed applying the EM algorithm on the estimation of a low-order model error covariance terms associated to a reduced order noisy dynamics.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Case I: Comparison of state RMSE evolution for EnKF/EnRTS-EM. Black solid/dashed/dotted/dash-dotted line denotes the analysis trajectory yielded after the 1st/2nd/3rd/5th E-step iteration out of 10 total iterations respectively. . . . . . . . . Fig. 2. Case I, Parameter RMSE against EnKF/EnRTS-EM iterations . . . . . . . . . . Fig. 3. Case I: Estimated parameter profile w.r.t to iterations of EnKF/EnRTS-EM (left columniterations #1(fig. a); #2(fig. c); #3(fig. e); #4(fig. g); # 5(fig. i)) compared to Aug-EnVar-EM (right column -iterations #1(fig. b); #2(fig. d); #3(fig. f); #4(fig. h); # 5(fig. j)). . . . .

Fig

  Fig. 4. Case I, Comparison of state RMSE evolution after the final final EM iterations. Black curve: the analysis trajectory yielded by EnKF within the E-step. Red cross points/Black square points: the forecast/analysis state yielded by EnKF in the E-step. Red star points: the smooth state yielded by EnRTS Smoother at observational time. Black dotted curve: the smooth trajectory yielded by integrating the initial smooth state using final estimated parameters. Red curve: the analysis trajectory yielded by Aug-EnVar-EM. . . . . . . . . . .

  FIG. 2: Case I, Parameter RMSE against EnKF/EnRTS-EM iterations

FIG. 5 :

 5 FIG.5: Case II, Comparison of state RMSE evolution after the final EM iterations. Black curve: the smooth trajectory yielded by integrating the initial smooth state after final E-step using final estimated parameters after final M-step. Red curve: the analysis trajectory yielded by Aug-EnVar-EM. Red dashed curve: the analysis trajectory yielded by Aug-EnKF-EM.

  . 4. Case I, Comparison of state RMSE evolution after the final final EM iterations. Black curve: the analysis trajectory yielded by EnKF within the E-step. Red cross points/Black square points: the forecast/analysis state yielded by EnKF in the E-step. Red star points: the smooth state yielded by EnRTS Smoother at observational time. Black dotted curve: the smooth trajectory yielded by integrating the initial smooth state using final estimated parameters. Red curve: the analysis trajectory yielded by Aug-EnVar-EM. . . . . . . . . . . Fig. 5. Case II, Comparison of state RMSE evolution after the final EM iterations. Black curve: the smooth trajectory yielded by integrating the initial smooth state after final E-step using final estimated parameters after final M-step. Red curve: the analysis trajectory yielded by Aug-EnVar-EM. Red dashed curve: the analysis trajectory yielded by Aug-EnKF-EM. . . . Fig. 6. Case II: Estimated parameter profile yielded from window starting at 3000s . . . . . . Case I: Comparison of state RMSE evolution for EnKF/EnRTS-EM. Black solid/dashed/dotted/dash-dotted line denotes the analysis trajectory yielded after the 1st/2nd/3rd/5th E-step iteration out of 10 total iterations respectively.

						Height
	10	1
	RMS(h)
	10	0
		600	700	800	900	1000 1100 1200 1300
						t
						Velocity
						Background
	10	-1			Observation
	RMS(u)					Analysis 1st Ite Analysis 2nd Ite Analysis 3rd Ite Analysis 5th Ite
	10	-2
		600	700	800	900	1000 1100 1200 1300
						t
	FIG. 1:
						41

  ). Case I, Comparison of state RMSE evolution after the final final EM iterations. Black curve: the analysis trajectory yielded by EnKF within the E-step. Red cross points/Black square points: the forecast/analysis state yielded by EnKF in the E-step. Red star points: the smooth state yielded by EnRTS Smoother at observational time. Black dotted curve: the smooth trajectory yielded by integrating the initial smooth state using final estimated parameters. Red curve: the analysis trajectory yielded by Aug-EnVar-EM.

					Height
					Height
	10 10	1 1
	RMS(h) RMS(h)
	10 10	0 0
	600	700	800	900	1000 1100 1200 1300
	t 1000 1100 1200 1300 t Velocity 900 Background Trajectory 700 800 Observation Analysis Trajectory Forecast EnRTS Analysis Smoothed Trajectory Aug-EnVar-EM Analysis t 2400 2500 2600 2700 2800 2900 3000 3100 3200 3300 600 10 -2 10 -1 t 10 -1 Velocity Background Observation EnKF/EnRTS-EM Analysis FIG. 4: 2400 2500 2600 2700 2800 2900 3000 3100 3200 3300 RMS(u) Aug-EnVar Analysis RMS(u) Aug-EnKF Analysis
					44

We omit the prior θ l in the expectation operator notation.

Q (46)

Acknowledgments.