Yin Yang
email: iamyangyin@gmail.com

Etienne Mémin

E M Émin

High-resolution data assimilation through stochastic subgrid tensor and parameter estimation from 4DEnVar

In this paper we explore a dynamical formulation allowing to assimilate high resolution data in a largescale fluid flow model. This large-scale formulation relies on a random modeling of the small-scale velocity component and allows to take into account the scale discrepancy between the dynamics and the observations. It introduces a subgrid stress tensor that naturally emerges from a modified Reynolds transport theorem adapted to this stochastic representation of the flow. This principle is used within a stochastic shallow water model coupled with an 4DEnVar assimilation technique to estimate both the flow initial conditions and the inhomogeneous time-varying subgrid parameters. The performance of this modeling has been assessed numerically with both synthetic and real world data. Our strategy has shown to be very effective in providing a more relevant prior/posterior ensemble in terms of the dispersion compared to other tests using the standard shallow water equations with no subgrid parameterization or with simple eddy viscosity models. We also compared two localization techniques. The results indicate the localized covariance approach is more suitable to deal with the scale discrepancy related errors.

Introduction

The numerical simulation of all the scales of geophysical fluid flows, up to the smallest dissipation scales, remains nowadays completely intractable due to the dimension of the problem to handle (which scales as the cube of the Reynolds number). For atmospheric flows for instance, the spatial scales range from several millimeters to thousands of kilometers. The temporal scale is also huge and goes from decades to seconds. Many of the physical processes related to the submeso-scale and microscale of the geophysical systems are left unresolved. However, due to nonlinear interactions between scales, it is important to understand the actions of those unresolved components to take them into account in a large-scale representation of the dynamics.

For several years there has been a growing interest in modeling the unresolved processes through random terms. In that respect, the simplest method consists in introducing random forcing terms into the flow dynamics equations [START_REF] Buizza | Stochastic representation of model uncertainties in the ecmwf ensemble prediction system[END_REF][START_REF] Leith | Stochastic backscatter in a subgrid-scale model: plane shear mixing layer[END_REF][START_REF] Mason | Stochastic backscatter in large-eddy simulations of boundary layers[END_REF][START_REF] Shutts | A kinetic energy backscatter algorithm for use in ensemble prediction systems[END_REF][START_REF] Slingo | Uncertainty in weather and climate prediction[END_REF]. This kind of methods is often coupled to the usual subgrid-scale (SGS) parameterization used in deterministic large-scale approaches where the statistically averaged flow on a coarse grid are solved directly and the effects of the small unresolved scale on the large resolved scales are modeled thereafter. Those models of the Reynolds stress reflect primarily the dissipation due to turbulence while stochastic subgrid closure terms have been proposed to model also the energy backscattering from small to large scales [START_REF] Frederiksen | Subgrid modelling for geophysical flows[END_REF].

Another way of considering the stochastic nature of the neglected scales is to study directly the governing differential equations driven by stochastic processes. To that end, a stochastic dynamics modeling is proposed in [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF]. In this approach, which constitutes the fundamental cornerstone to this article, the author proposed to introduce the stochasticity from the very beginning in the dynamics physical derivation.

The basic idea is built on the assumption that the Lagrangian fluid particle displacement results from a smooth velocity component and a fast oscillating stochastic velocity component uncorrelated in time,

Xt = Xt 0 + t t 0
w(Xs, s)ds + t t 0 σ(Xs, s)dBs, (1) c 0000 Tellus with Xt denoting the flow map at time t, w the smooth fluid particles velocity corresponding to the sought "large-scale" solution, and dBt a d-dimensional Wiener increment (d=2,3). The term σ(x, t)dBt represents the fast-changing fluctuation associated to turbulence; it has a much lower time-scale than w, the smooth component. Both σ and w require to be solved or modeled. Stochastic calculus principles provides then the rules to derive the expressions describing the evolution of functions of such flows.

Large-scale models are usually assessed and calibrated by comparing simulated results and data. Opposite to the sparse heterogeneous network of in situ measurements, satellite image data exhibit a spatial resolution that is finer than the resolved scales at which the state variables are solved through a filtered dynamics. This still increasing scale discrepancy due to the rapid technological progress of satellite sensors advocates the setup of methods enabling to calibrate the numerical flow model from high-resolution data. The data assimilation (DA) approach can indeed serves this purpose. Note that image data assimilation is a research topic drawing more and more attention since the last decades and several strategies have been proposed to deal with these difficulties, [START_REF] Beyou | Weighted ensemble transform kalman filter for image assimilation[END_REF][START_REF] Bereziat | Coupling dynamic equations and satellite images for modelling ocean surface circulation[END_REF][START_REF] Chabot | Accounting for observation errors in image data assimilation[END_REF][START_REF] Corpetti | Pressure image assimilation for atmospheric motion estimation[END_REF][START_REF] Cuzol | A stochastic filtering technique for fluid flow velocity fields tracking[END_REF][START_REF] Huot | Optimal orthogonal basis and image assimilation: Motion modeling[END_REF][START_REF] Papadakis | Variational Assimilation of Fluid Motion from Image Sequence[END_REF][START_REF] Souopgui | Assimilation d'images pour les fluides géophysiques[END_REF][START_REF] Titaud | Assimilation of image sequences in numerical models[END_REF].

Traditionally, data assimilation techniques are widely used to recover as exactly as possible the initial condition given a dynamics model, a certain span of time called assimilation window and some observations across this window [START_REF] Talagrand | Assimilation of observations, an introduction[END_REF][START_REF] Blum | Data Assimilation for Geophysical Fluids[END_REF]. The initial condition is crucial for Numerical Weather Prediction (NWP) systems but the inference of the models' parameters is also tractable through the data assimilation process. In a broader sense, the assimilation technique discussed in this article can be viewed as a learning process of a dynamic model on the basis of the observations available by adjusting the model's parameters.

Data assimilation techniques are classically categorized into two groups. The first group, referred to as variational assimilation techniques, is formulated from the optimal control theory and variational calculus [START_REF] Lions | Optimal control of systmes governed by PDEs[END_REF][START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF]). The second group of "ensemble methods" originates from Monte Carlo Bayesian filtering principles [START_REF] Evensen | Sequential data assimilation with a non linear quasi-geostrophic model using Monte Carlo methods to forecast error statistics[END_REF][START_REF] Gordon | Novel approach to nonlinear/non-gaussian bayesian state estimation[END_REF]. The parameter estimation problem can also be divided in the same spirit into those two groups.

In the former approach, the estimation consists in considering the unknown parameters as a control parameter as proposed by [START_REF] Navon | Practical and theoretical aspects of adjoint parameter estimation and identifiability in meteorology and oceanography[END_REF]. Consequently, the objective of the variational system is to find the optimal initial condition and the associated optimal parameters for which a system trajectory fits best the available observations. Several works in this vein can be found in [START_REF] Cacuci | Computational Methods for Data Evaluation and Assimilation[END_REF][START_REF] Navon | Practical and theoretical aspects of adjoint parameter estimation and identifiability in meteorology and oceanography[END_REF][START_REF] Zhu | Impact of parameter estimation on the performance of the fsu global spectral model using its full-physics adjoint[END_REF][START_REF] Pulido | Gravity wave drag estimation from global analyses using variational data assimilation principles. i: Theory and implementation[END_REF][START_REF] Kazantsev | Sensitivity of a shallow-water model to parame-ters[END_REF]. The latter approach is a direct extension of standard ensemble Kalman filter methods. Within this framework, one can choose freely the specific ensemble filtering schemes in terms of specific applications [START_REF] Zupanski | Model error estimation employing an ensemble data assimilation approach[END_REF][START_REF] Tong | Simultaneous estimation of microphysical parameters and atmospheric state with simulated radar data and ensemble square root kalman filter. part i: Sensitivity analysis and parameter identifiability[END_REF][START_REF] Yang | Using the ensemble Kalman filter to estimate multiplicative model parameters[END_REF][START_REF] Bocquet | Joint state and parameter estimation with an iterative ensemble Kalman smoother[END_REF][START_REF] Ruiz | Estimating model parameters with ensemble-based data assimilation: Parameter covariance treatment[END_REF][START_REF] Tandeo | Offline parameter estimation using EnKF and maximum likelihood error covariance estimates: Application to a subgrid-scale orography parametrization[END_REF][START_REF] Sawada | Assimilating high-resolution winds from a Doppler lidar using an ensemble Kalman filter with lateral boundary adjustment[END_REF].

In this article, we propose to use an ensemble-based variational scheme which gathers the advantages of both the variational assimilation schemes and the ensemble-based sequential methods [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF]. In addition to the obvious opportunity of not requiring any adjoint dynamical models, ensemblebased schemes provide other benefits especially when considering stochastic dynamics. As a matter of fact ensemble prediction systems are becoming more and more crucial in meteorology and oceanography. In that respect, under-dispersion of the ensemble greatly jeopardize forecasting quality [START_REF] Slingo | Uncertainty in weather and climate prediction[END_REF]. Such underestimation of the ensemble spread is usually related to the unawareness of the uncertainties associated with the dynamic system. Recent studies show that the quality of the ensemble spread can be significantly improved by not only perturbing the initial condition but also some parameters related to the dynamics [START_REF] Wei | The impact of initial spread calibration on the relo ensemble and its application to lagrangian dynamics[END_REF][START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF]. Those facts advocate the setup of physically relevant stochastic dynamical models coupling the random production terms and the dissipation term.

Within this paper, we focus on the specification through highresolution image data assimilation of the dynamics model uncertainties. A fully justified argument on the derivation of the stochastic fluid flow dynamics is not the main purpose here. This framework is briefly reviewed in section 2, then we introduce on this basis the stochastic shallow water system. In section 3, we investigate methods for estimating a central object called the quadratic variation tensor from the data. Model and experimental settings, results and discussions with some concluding remarks are provided in the following sections.

Stochastic modeling

Analogously to the standard deterministic case, the derivation procedure from the physical conservation laws of the Navier-Stokes equations is based primarily on the Reynolds transport theorem (RTT). In our context the RTT consists of a stochastic expression of the rate of change of a given quantity (mass, momentum and energy) in a material volume. A comprehensive derivation process of this theorem [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] is presented in appendix A for reader convenience.

Applied to scalar function, q, the expression of the RTT within a volume transported by a stochastic flow (1) with an incompressible small-scale velocity component (∇ • σ = 0), is given by:

d V(t) qdx = V(t) [dtq + ∇ • q(w - 1 2 ∇ • a w ) dt - 1 2 ij ∂ ∂xi aij ∂ ∂xj q dt + ∇q • σdBt]dx. (2)
This modified RTT involves the time increment of function q (the differential of q at a fixed point) instead of the time derivative. A diffusion operator emerges also naturally. For clarity's sake, this term is designated as "subgrid stress tensor" following the protocols of large eddies simulation (LES). However, its construction is quite different. It is not based on Boussinesq's eddy viscosity assumption nor on any structural turbulence models [START_REF] Sagaut | Large eddy simulation for incompressible flows: an introduction[END_REF] but arises from stochastic calculus rules. It expresses the mixing process exerted by the fast oscillating velocity component. This diffusion term is directly related to the small-scale component through the variance tensor, a, defined from the diagonal of the uncertainty covariance:

a(x)δ(t -t ′ )dt = E (σ(x, t)dBt) σ(x, t ′ )dB t ′ T .
The variance tensor corresponds to an eddy viscosity term (with units in m 2 s -1 ).

A corrective advection term, w = w-1/2∇•a, appears also in the stochastic RTT formulation. This correction expresses the influence of the small scales inhomogeneity on the large ones. A drift is engendered from the regions associated with maximal variance (maximal turbulent kinetic energy -TKE) toward the area of minimum variance (e.g. minimal TKE). For homogeneous noise (i.e. within an homogeneous turbulence assumption), the variance tensor is constant and this corrective advection does not come into play.

Through this modified RTT, stochastic versions of the mass and momentum conservation equations can be (almost) directly derived. Incompressibility conditions can for instance be immediately deduced from the RTT applied to q = 1 and the flow jacobian (J):

d V(t) dx = V(t) ∇ • w(x, t)dtdx = V(t 0 ) d(J(Xt(x), t))dx = V(t 0 ) (J∇ • w) (Xt(x), t)dtdx. (3)
Together with the incompressibility of the random term, the incompressibility condition reads thus:

∇ • σ = 0 and ∇ • w = 0. (4)
Imposing, a priori, an incompressibility condition to the largescale component, w, simplifies the previous constraint as:

∇ • σ = 0 and ∇ • w = ∇ • (∇ • a) T = 0. (5)
Note that for an homogeneous unresolved velocity, the last condition is naturally satisfied, since this unresolved velocity component is associated with a constant variance tensor, a. It should be stressed that condition (4) corresponds to a weaker form of the incompressibility constraint in the sense that it does not impose an incompressibility condition to the large-scale component but only to the effective drift. It has the advantage not to include a cumbersome condition on the variance tensor.

In the following section we will briefly present how this principle can be applied to derive a stochastic expression of the shallow water equations. Those expressions are derived from the ap-plication of the stochastic RTT to the linear momentum conservation principle. Interested readers may consult [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] for a thorough development or [START_REF] Resseguier | Geophysical flows under location uncertainty, part I random transport and general models[END_REF]b;[START_REF] Leith | Stochastic backscatter in a subgrid-scale model: plane shear mixing layer[END_REF] for the derivation and the analysis of general geophysical flow models and their quasigeostrophic approximations.

The stochastic 2D nonlinear shallow water model

Considering the usual assumptions underlying the constitution of the Shallow water model: the depth is shallow compared to the horizontal dimension, the friction forces and the vertical shear on the velocity horizontal components can be neglected, the external force is only due to gravity, we get the following 2D momentum equations for the horizontal velocity components (where superscript "h" stands for the horizontal coordinates)

∂u h ∂t + ∇ h • (uu T ) - 1 2 (ij) h ∂x i ∂x j (a h ij u) = - 1 ρ ∇ h p, (6)
u h = u h - 1 2 ∇ • a h (7) ∇ h dp = -ρ(w h ∇ h )(σdBt) h . ( 8
)
The last equation arises from a balance of the advected Brownian term and the random pressure, dp, associated to the random small-scale component. This equation stems from a separation principle between the smooth terms and the Brownian terms. It can be noted that in those equations the noise components are assumed to depend only on the horizontal components. Thus the diffusion tensor as well as the variance tensor do not depend on depth. These equations are complemented by the three incompressibility constraints (5). The noise component (and its related variance tensor) as well as the horizontal velocity component can be interpreted as depth-averaged components of the 3D velocity and uncertainty components respectively. Integrating the incompressibility constraints (∇ • w = 0) and (∇ • σ = 0) along the depth from the bottom h b to the surface hs we get from the effective drift expression (7):

∇ h • u(hs -h b ) = w| hs -w| h b , and
(σdBt) z | hs = (σdBt) z | hu .
We may noticed that for (∇ • ∇ • a = 0) we recover the traditional balance between the divergence of the horizontal averaged velocity components and the vertical motions at the interfaces. In the same way as in the traditional setting, considering the difference of the kinematic boundary conditions at the surface and the still bottom, and introducing the previous relation for the depth variable, h = hs -h b , the shallow-water continuity equation reads:

dth+ ∇ h • h u h - 1 2 (i,j) h ∂x i (aij∂x j h) dt+∇ h •(hσdBt) = 0. (9)
It should be remarked that this continuity equation involves the horizontal components of the (divergence free) effective advection term:

u h = u h -1 2 ∇ h •a h . The term (h∇ h •∇ h •a h
) plays the role of an additional forcing condition that ensues from the small-scale inhomogeneity. For smooth enough fluctuation variances this forcing is negligible.

Introducing the hydrostatic balance (which is well respected in the shallow water model due to a negligible vertical acceleration) the pressure (after depth integration and with a boundary condition provided by the atmospheric pressure at the water-air interface) and its gradient reads respectively p|z = ρg(hs -z) + p|atm, and ∇p|z = ρg∇hs (10)

We get hence the following horizontal linear momentum conservation equation:

∂u ∂t + ( u h • ∇ h )u h - 1 2 (ij) h ∂x i (aij∂x j u h ) = -g∇ h hs.
(11) It is important to outline that the free surface acts here as a random forcing term on the momentum equations. It is driven by the stochastic evolution of the height variable (9). For a deterministic height evolution such as the one associated to the mean depth evolution conditionally to a deterministic velocity field, the momentum equation would provide the velocity fields attached precisely to this particular (mean) instance of the evolving free-surface.

The uncertainty formalism allows us to obtain the 2D shallow water equation in the form of (9, 11). It can be observed that the dissipative terms emerged both in the continuity equation and in the momentum conservation equation. It is similar in spirit to the eddy viscosity assumption used in LES and Reynold average numerical simulation (RANS). Nevertheless, our subgrid-scale model, takes its origin in the stochastic representation of the uncertainty associated to the model itself. No extra assumption is needed. However, their similarities indeed bring us another perspective for the interpretation of the subgrid-scale effect. For instance, it can be shown [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] that setting the variance tensor as:

a h = 2ν h I h , (12)
with the eddy viscosity coefficient

ν h = C∆x∆y ∂xu 2 + ∂yv 2 + 1 2 ∂yu + ∂xv 2 , ( 13
)
corresponds, for smooth enough rate of strain tensor norm, to the famous [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. the basic experiment*[END_REF] subgrid model.

Conservative form

The conservative form of shallow water equations is usually more useful when implementing the finite-volume method. It can be inferred directly by combining equation ( 11) and equation (9). Alternatively, one can deduce such a form by integrating equation ( 9) along the depth from the bottom h b to the surface hs. Finally, we have, (note we dropped in the following the subscript h for sakes of clarity; all the vectorial variables and spatial derivative operators represent planar vectors and derivations w.r.t. horizontal components respectively)

dt(hu)+∇• (hu u T ) + 1 2 gη 2 I dt - (i,j) 1 2 ∂x i aij∂x j (hu) -aij ∂x j h∂x i u dt + u∇• hσdBt = 0, (14)
As previously indicated, if the focus is only on the ensemble average depth evolution (denoted by h) or if the noise is assumed to live on iso-height surfaces (i.e σdBt • ∇h = 0) , the above equation greatly simplifies as it boils down to a deterministic system:

∂ h ∂t + ∇ • ( h u) - 1 2 (i,j) ∂x i aij∂x j h = 0, (15a)
∂( hu) ∂t +∇• ( hu u T ) + 1 2 gη 2 I - (i,j) 1 2 ∂x i aij∂x j (hu) -aij∂x j h∂x i u = 0. (15b)
This system, as opposed to the usual deterministic Shallow water model, still includes additional subgrid terms that depend on the variance of the small-scale random component.

1D stochastic shallow water equation

The tests shown in the next sections are carried out with a 1D stochastic shallow water equation since it is less timeconsuming but remains representative of the considered dynamics. The 1D stochastic shallow water system simply reads,

dth + ∂x h(u - 1 2 ∂xaxx u - 1 2 ∂x(axx∂xh) dt + ∂xh(σdBt)x = 0, (16a) dt(hu)+ ∂x(hu u + 1 2 gη 2 ) - 1 2 ∂x axx∂x(hu) -axx∂xh∂xu dt + u∂xh(σdBt)x = 0, (16b)
and the mean depth and its associated horizontal fields satisfy,

∂t h + ∂x( h u) - 1 2 ∂ 2 ∂x 2 (axx h) = 0, (17a) ∂t( hu) + ∂x( hu u + 1 2 gη 2 ) - 1 2 ∂x axx∂x( hu) -axx∂x h∂xu = 0. ( 17b
)
As stated before, one of the main motivations of our approach is to propose methods for the assimilation of fine resolution image data. An interesting idea consists thus in estimating the diagonal quadratic variation tensor, a, or to model the (plain) diffusion tensor, σ, from the statistical variations of small scale velocity measurements extracted from image data [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF].

In the following experiments, we will restrain our interests to the horizontal mean flow fields described in Eqs. ( 17) and (15). However the use of the full stochastic model with the noise forcing terms is also possible alternative [START_REF] Miller | Data assimilation into nonlinear stochastic models[END_REF] that will be explored in a future study for data assimilation methods expressed in a stochastic filtering context.

Estimation of the quadratic variation tensor from data

Estimation by explicit formulation

Several immediate and simple strategies can be used to define the variance tensor. As detailed previously, an explicit form of this tensor can be given so as to recover classical subgrid models. However the specific statistical nature of this tensor allows us to specify a variety of empirical forms for this tensor [START_REF] Kadri Harouna | Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling[END_REF]. For instance, it can be defined from local velocity covariances matrices estimated on local spatial or temporal neighborhood. The previous estimation methods are easy to implement, but constitute only rough approximations as they rely either on the resolved components or on the flow velocity measurements, which are in practice largescale observations [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF]. A more accurate way of estimating the parameters is through a data assimilation process.

Estimation through data assimilation process

Data assimilation techniques consist in estimating the evolution of a given state variables, denoted x ∈ R n , by associating measurements, Y ∈ R m , of this system with a dynamics model. The evolution model will be denoted in a general way as

∂tx(t, x) + M(x(t, x), a(t, x)) = 0, (18)
x(t0,

x) = x 0 f (x) + ω(x). (19)
It involves a non autonomous differential operator M -usually nonlinear -that depends in our case on the variance tensor, a.

The second equation defines the initial state, x(t0, x), which is assumed to be known up to a zero Gaussian noise, ω, of covariance B. The integration of the initial forecast guess, x 0 f , is called the forecast trajectory, x f (t, x). It is usually fixed from a previous assimilation result on an anterior temporal window. Note that depending on the context, x f is also referred as the background. The observations, Y, and the state variable x are linked together through an observation operator H

Y(t, x) = H(x(t, x)) + ǫ(t, x). ( 20
)
The error, ǫ, between the measurements and the state is assumed to be a zero mean Gaussian random field with a covariance tensor, R. This noise is assumed to be uncorrelated from the state variable.

Optimal parameter can be found through data assimilation techniques along with the system state. The strategies are formulated as the minimization with respect to control variables (here the initial condition and the variance tensor) of an energy functional encoding a distance between the data and a state variable trajectory with low values of the parameters norm.

In our approach, we express the cost function in an incremental form [START_REF] Courtier | A strategy for operational implementation of 4d-var, using an incremental approach[END_REF], with a Mahalanobis norm f 2 A = (f, A -1 f ) where A -1 is an inverse covariance matrix and (f, g) is the L 2 inner product :

J(δx0, δa) = 1 2 δx0(x) 2 Bx + 1 2 δa 2 Ba + 1 2 t f t 0 ∂xH δx(t, x, a) -D(t, x, a) 2 R dt. ( 21
)
The initial analysis filed (denoted by subscript a) are related to the increments to be sought through

x 0 a = x 0 f + δx0, aa = a f + δa. ( 22
)
The innovation vector D(t, x, a) is defined as:

D(t, x, a) = Y(t, x) -H ϕt(x f , a f ) , (23)
where ϕt(x0, a0) is the flow map (viewed here as a function of the initial condition and of the initial parameters). In addition the increment δx(t, x, a) is driven by the tangent-linear model linearized around the current trajectory x f ,

δx(t, x, a) = ∂xϕt(x f , a f )δx0.
We are interested in applying a particular DA method, namely a group of ensemble variational scheme (4DEnVar), to our state/parameter estimation problem. 4DEnVar [START_REF] Liu | An Ensemble-Based Four-Dimensional Variational Data Assimilation Scheme. Part I: Technical Formulation and Preliminary Test[END_REF]; [START_REF] Buehner | Intercomparison of Variational Data Assimilation and the Ensemble Kalman Filter for Global Deterministic NWP. Part I: Description and Single-Observation Experiments[END_REF]) features a 4D ensemble-based covariance fields in a variational system. The strategy described here is a direct extension of the scheme discussed in [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF]. The extension of 4DEnVar for parameter estimation is done by the so-called augmented state vector technique. If we assume that the parameter variables are random variables, then the extension of state-space based methods to parameter-space is straightforward by encapsulating the parameters. The state augmentation technique has shown to be effective in diverse applications [START_REF] Tong | Simultaneous estimation of microphysical parameters and atmospheric state with simulated radar data and ensemble square root kalman filter. part i: Sensitivity analysis and parameter identifiability[END_REF]; [START_REF] Kang | variable localization in an ensemble kalman filter: Application to the carbon cycle data assimilation[END_REF]; Simon et al. (2015); [START_REF] Sawada | Assimilating high-resolution winds from a Doppler lidar using an ensemble Kalman filter with lateral boundary adjustment[END_REF]). The parameter estimation is also related to the estimation of model error from ensemble-based approach. The model error can take simple forms of model bias [START_REF] Dee | Bias and data assimilation[END_REF][START_REF] Baek | Local ensemble Kalman filtering in the presence of model bias[END_REF] or functional forms through parameterization.

Augmented 4DEnVar algorithm

The constitution of a corresponding joint estimation algorithm is illustrated below within an augmented state context. First we need to define s0 as the initial augmented state vector [x0, a0] ∈ R (n+p) where p is the degree of freedom of a. In an incremental framework, we are looking for δs0 defined by s0s f . The cost function in terms of the augmented state increment vector δs0 reads:

J(δs0) = 1 2 δs0 2 B f s + (24) 1 2 t f t 0 ∂xH ∂xϕ t(sf )δs0 -D(t, x, a) 2 R dt,
where the symbol • denotes the augmented version of the tangent linear model or the tangent observation model respectively. Firstly, the extended tangent linear model operator ∂xϕ t is composed of the parameter evolution model as well as the tangent linear dynamic model operator ∂xϕt. Secondly, in a fully and directly observable case, the observation operator corresponds to the identity matrix. In our case, the validity of a hinges on the quality of the analysis trajectory driven by the model, therefore, rather than to simply fit to the noisy empirically observed a values with few samples, we will assume that a is non-observable. Let us note that a prior information on the variance tensor a0 is included in the background regulation term. The observation operator or its linear tangent expression takes the form of (Im×n, 0m×p) T . The augmented state ensemble, denoted by S f ∈ R (n+p)×N , is obtained by perturbing the background s f ; N is the number of ensemble members. We defined the ensemble anomaly matrix as,

A ′ s := 1 √ N -1 (S 1 f -S f , . . . , S N f -S f ), A ′ s ∈ R (n+p)×N
(25) where the operator S takes the ensemble mean and the jth perturbed ensemble member is noted

S j f = [x j , a j ]. The background covariance B f
s is approximated by the ensemble anomaly matrix:

B f s = A ′ s A ′ s T = Bxx B T ax Bax Baa , (26)
where the sub-matrix Bxx and Baa correspond to the autocovariance matrix and Bax contains the cross covariance between the state and the parameter. As suggested by several studies, a localization procedure enables to robustify the construction of ensemble-based covariance to sampling errors. Due to the parameters' spatial structure, it is necessary to introduce similar localization to the parameter-related covariances. In this study we have tested two different localization techniques referred to as 'localized covariance" and "local analysis" respectively. The former relies on a Schur product between the global covariance matrix and an isotropic correlation function,

P f s = Cs ⊙ A ′ s A ′ s T = Cxx C T ax Cax Caa ⊙ Bxx B T ax Bax Baa , (27)
to filter out the long range spatial correlations. The spatial correlation matrix C•• is defined from a polynomial approximation of a Gaussian function with compact support and a hard cut-off [START_REF] Gaspari | Construction of correlation functions in two and three dimensions[END_REF]. Same values of the cutoff distance are adopted for all C••. The latter consist in bounding -or localizing -directly the region of influence of a given observation on the grid to a restricted area. The 4DEnVar consists in applying a change of variable with the square-root of the error covariance matrix, by setting δs0 = Uδz0, where δz0 is the control weight vector and U is the square root of P f s calculated through a spectral decomposition [START_REF] Buehner | Ensemble-derived stationary and flow-dependent background-error covariances: Evaluation in a quasi-operational NWP setting[END_REF] or simply A ′ s in the case of no localization. We can define the following scheme: U = A ′ s D. The control variable transformation is also relevant as a preconditioning technique. The cost function in terms of control variable reads,

J(δz0) = 1 2 δz0 2 + 1 2 t f t 0 ∂xH B 1 2 t δz0 -D(t, x, a) 2 R dt.
where B

1 2 t = ∂xϕ t(sf )U. The B 1 2
t contains a set of samples based on an expression of the ensemble integrated in time. Note that for the Shur-product based localization of this empirical matrix, we implemented ( ∂xϕ t(sf )A ′ s )D rather than ∂xϕ t(sf )(A ′ s D) (which comes to consider a time invariant localization).

In order to handle high nonlinearity either related to the model or the observations, many authors advocated the use of an iterative approach [START_REF] Gu | An Iterative Ensemble Kalman Filter for Multiphase Fluid Flow Data Assimilation[END_REF][START_REF] Kalnay | Accelerating the spin-up of Ensemble Kalman Filtering[END_REF][START_REF] Sakov | An Iterative EnKF for Strongly Nonlinear Systems[END_REF]Bocquet and Sakov, 2013a). The algorithm in [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF] is characterized by an iterative resampling strategy of the prior ensemble and an update of the ensemble sensitivity matrix within the increment minimization steps. It adopts the form of a Gauss-Newton scheme with nested loops that shows the relevance of this approach in particular in the case of a badly known background information and incomplete observations. In our current problem, with the background information bounded by the coarse resolution and the absence of observed parameter fields, it is indeed appropriate to introduce such iterative approach. Here, the iterative update concerns two parts: on the one hand, the innovation vector of the current outer loop is updated by the analysis trajectory based on the analysis state of the previous outer loop driven by the nonlinear model; on the other hand, the ensemble anomaly matrix is updated through the posterior ensemble fields. The ensemble update can either be derived from the posterior ensemble built on perturbed observations or by a direct transformation of the prior perturbations.

Perturbed observation with global analysis

In the first approach, at the kth outer loop iteration, a parallel realization of minimization with regard to each member of the initial ensemble is conducted,

δS (j),k 0 = UδZ (j),k 0 , j = 1, . . . , N. J(δZ (j),k 0 ) = 1 2 δZ (j),k 0 2 + (28) 1 2 t f t 0 ∂xH ∂xϕ t(s k 0 )UδZ (j),k 0 -D (j),k (t, x) 2 R dt,
where D (j) denotes the jth innovation vector calculated from perturbed observation field and Z0 denotes the control weight matrix. Finally, the updated initial ensemble field and its perturbation matrix read:

S (j),k+1 0 = S (j),k 0 + U δZ0 (j),k , j = 1, . . . , N. A ′k+1 s = 1 √ N -1 S (1),k+1 0 -S k+1 0 , . . . , S (N),k+1 0 -S k+1 0 . ( 29
)
3.2.1.2. Ensemble transformation with local analysis The square root transformation approach corresponds to a mean preserving transformation as used in Ensemble Transform Kalman filter (ETKF). The updated background ensemble anomaly matrix reads,

A ′k+1 s = A ′k s {I+ t f t 0 B 1 2 T t (∂xH) * R -1 ∂xH B 1 2 t dt} -1 2 ,k V.
(30) The transformation matrix is the square root of inverse Hessian1 . The arbitrary orthogonal matrix V is used to center the posterior ensemble on the updated initial condition/analysis. In this approach a single minimization process with respect to the background state is sufficient since the ensemble does not need to be updated explicitly. Finally the updated initial ensemble fields are,

S (j),k+1 0 = s k+1 0 + √ N -1A ′(j),k+1 s , j = 1, . . . , N (31)
where s k+1 0 corresponds to the updated augmented vector,

s k+1 0 = s k 0 + A ′k s δz0 k . ( 32
)
Since the minimization is done in the ensemble space, we use here the local analysis approach. This computationally efficient approach proposed in [START_REF] Ott | A local ensemble Kalman filter for atmospheric data assimilation[END_REF] and [START_REF] Hunt | Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter[END_REF] consists to lead the analysis step through an EnKF procedure in a local spatial region. The minimization and the error covariance update steps are done around each grid point considering the flow representation and the observations within a region of small size. Within the local domain, observation error is gradually increased when further away from the analysis grid point [START_REF] Greybush | Balance and Ensemble Kalman Filter Localization Techniques[END_REF]. Note that the ensemble forecast step must be done globally with the full non-linear dynamic model.

Deterministic global analysis

We also tested another strategy capable of generating global deterministic ensemble anomaly matrix without perturbing the observations. This formulation is close in spirit with the DEnKF method proposed in [START_REF] Sakov | A deterministic formulation of the ensemble Kalman filter: an alternative to ensemble square root filters[END_REF]. In this approach, the posterior ensemble is obtained by using half of the Kalman gain therefore preventing the level of posterior error from underestimation in the context of unperturbated observations. Here we use the same format for the analysis ensemble which can be found in [START_REF] Fairbairn | A comparison of 4DVar with ensemble data assimilation methods[END_REF]:

S k+1 0 ← S k+1 0 + 1 2 (S k 0 -S k 0 ) + 1 2 (S k+1 0 -S k+1 0 ). ( 33
)
Let us note, that the iterative schemes described above are also used in the ensemble cycling procedure to move from one assimilation window to the next.

Parameter evolution and uncertainties

It is essential to realize that the parameters involved in geophysical models can be global or local in time and/or space. The parameters are generally time-dependent and their values may change when a new observation occurs. However, in this work we will make the assumption that the parameters exhibit only slow temporal variation compared to the model state variation. In DA approaches, it is even custom to assume a constant value of parameters during the assimilation window. This consists in forecasting the analysis state, at the end of the ith assimilation window, to the beginning of the i+1th cycle with the optimum parameter a i a :

x i+1 f = ϕi+1(x i a , a i a ), (34)
and to propagate the parameters using a persistent model,

a i+1 f = a i a . (35)
However, this persistent model can be a source of instability to the model integration with limited ensemble numbers. [START_REF] Gottwald | A mechanism for catastrophic filter divergence in data assimilation for sparse observation networks[END_REF] explained this instability by the contribution of the unrealistic large covariance terms. In order to damp the rapid temporal change of model parameter, we set,

a i a = a i f + βi δa i , (36)
where δa i is the optimizer obtained at the ith cycle and βi is a damping coefficient that needs be tuned in order to maximize the ensemble spread while not causing filter divergence based on trial-and-error principle.

In the framework of the augmented state technique, the parameters are treated as random variables. So their distribution must be set implicitly or explicitly. In our case we can use the estimators proposed in section 3.1. as prior values. The errors associated with those prior values of a are assumed to be centered Gaussian fields, with variance given by an uncertainty representation of the Smagorinsky model:

var(a) = C 1 2 ij ||∂x i δw j + ∂x j δw i ||I, ( 37
)
note δw is defined as variation of high resolution velocity wwk on local spatial region, C(x k ), centered on the coarse grid points indexed by k.

Ensemble Inflation

Ensemble methods typically require an explicit variance inflation scheme in order to maintain a relevant ensemble spread. Although such scheme is widely used in practice by either additive or multiplicative inflation formulations, it reflects a lack of knowledge on how to address properly the model errors responsible for the maintenance of the ensemble spread.

Here we adopt the multiplicative inflation [START_REF] Anderson | A monte carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts[END_REF] by multiplying the ensemble anomaly matrix with a tunable factor γ ranging from (1.0, 1.3) for each cycle. This strategy is applied to the standard model as a reference run but not to the stochastic model. Because one advantage of our stochastic model is to provide an inflation-free DA scheme. More precisely, our combined approach discussed in section 2. and 3.2. can be seen as an implicit inflation scheme where the corresponding scale discrepancy errors are dealt by the stochastic dynamics.

However, as shown in next section, the ensemble transformation approach with local analysis is particularly unable to produce sufficient ensemble spread leading to the filter divergence. To tackle such bad behavior, we employed a relaxation scheme in which part of the posterior ensemble perturbation is relaxed back to the prior perturbations [START_REF] Zhang | Impacts of initial estimate and observation availability on convective-scale data assimilation with an ensemble kalman filter[END_REF].

S k+1 0,relaxed = S k+1 0 + (1 -α)(S k 0 -S k 0 ) + α(S k+1 0 -S k+1 0 ). (38)
with 0 α 1. This approach is similar to the conditional covariance inflation method used in [START_REF] Aksoy | Ensemblebased simultaneous state and parameter estimation in a twodimensional sea-breeze model[END_REF].

Parameter estimation effect

It can be emphasized that although the estimation of the initial condition and the parameter proceed simultaneously, the estimation process of the initial condition is nevertheless not directly related to the auxiliary covariance terms (Baa, Bax) partly associated to the parameter a. The a, however, interferes in the quality of the analysis trajectory in two ways: firstly, it explicitly manifests itself on the model integration in the form of the estimated parameter a; secondly, it implicitly affects the analysis state through the ensemble spread from which the propagation of the ensemble perturbation matrix is calculated. Depending on the context, either of the two factors can exhibit large impact to the result [START_REF] Ruiz | Parameter estimation using ensemblebased data assimilation in the presence of model error[END_REF][START_REF] Koyama | Reducing forecast errors due to model imperfections using ensemble kalman filtering[END_REF]. Note that in our ensemble-based 4DVar approach, the parameter is not only related to the forecast-analysis cycling process, but also can affect the ensemble fields generated between two consecutive outer loops.

Parameter identifiability

Before we carry out the experiments, it is important to assess the parameter identifiability. This is especially pertinent in our case as the quadratic variation tensor a is non-observable. This concept, following the definition in [START_REF] Navon | Practical and theoretical aspects of adjoint parameter estimation and identifiability in meteorology and oceanography[END_REF], can be viewed as a criteria to decide whether or not the parameter of interest can be inferred from the data. With an ensemble-based method, an efficient way of examining the parameter identifiability consists in calculating the correlation coefficient [START_REF] Tong | Simultaneous estimation of microphysical parameters and atmospheric state with simulated radar data and ensemble square root kalman filter. part i: Sensitivity analysis and parameter identifiability[END_REF].

cor(x, a) = cov(x, a) var(x)var(a) . (39)
A strong correlation suggests that any changes in the parameter space will heavily affect the state space through the crosscovariance part of the error covariance matrix and vice versa.

Model and experimental settings

In order to evaluate the methods proposed, we have first carried out experiments on a one-dimensional stochastic shallow water model (Eqs.( 17)) with synthetic data. We extended then this procedure on a 2D stochastic shallow water model (Eqs.( 15)) with both synthetic and real image data.

Model numerical scheme

The systems studied here concentrates on the large scale evolution of an average depth. We note that the subgrid-scale model associated with the quadratic tensor a emerges both in the continuity equation and the momentum conservation equation. The finite volume method is directly applicable to these equations [START_REF] Bradford | Finite-volume model for shallow-water flooding of arbitrary topography[END_REF]. Indeed, the numerical flux formulation of the new subgrid-scale terms is analog to the viscous terms.

Experimental settings

Several experiments have been conducted to evaluate the performance of the stochastic model. We first tested the 1D stochastic shallow water equation. The 1D domain is of length Lx = 6000 km with the initial surface height h(x, 0).

h(x, 0) = H0 - f U0x g + Aξ, (40)
where H0 = 5000m, f = 1.03 × 10 -4 s -1 is the Coriolis parameter, U0 = 40m/s and g is the gravity acceleration. Let us remark that though the Coriolis force is not described in the original formulation for sake of clarity, it can be added without any difficulty to the model (see [START_REF] Resseguier | Geophysical flows under location uncertainty, part I random transport and general models[END_REF]b;[START_REF] Leith | Stochastic backscatter in a subgrid-scale model: plane shear mixing layer[END_REF]). An additional noise is considered with A = 10 as its amplitude. And ξ is a random Gaussian covariance field computed as the spectral methods discussed in [START_REF] Evensen | The ensemble Kalman filter: Theoretical formulation and practical implementation[END_REF] (with decorrelation length equals to 20%L). The initial velocity field is inferred from the geostrophic relation. The true state used to construct the finer observation is obtained by integrating the standard shallow water equation with the initial condition (40) on 401 grid points. The coarse-resolution background state with the same initial condition is simulated with the stochastic model on only 101 grid points. The time step ∆t is set to 37.5s for the true state model and to 150s for the coarse-resolution model in order to satisfy the different CFL conditions. The synthetic observation is extracted from the true state every 600s by adding i.i.d Gaussian noise to the true state. Each assimilation window is set as 1800s therefore containing 3 observations. In order to maintain the balance of state variables for ensemble fields, the analysis started at 6000s. For the 2D case, a very similar setup is implemented. The length of the domain is the same as Lx in 1D case. The width is Ly = 2400 km. The initial height surface is a 2D field with the extra error term as a 2D random Gaussian field. The background state resolution is 21×51 and the observation resolution is 81× 201. The time step ∆t is set to 10s for the true state model and to 50s for the coarse-resolution model.The synthetic observations are sampled every 150s from the true state. Each assimilation window is set as 450s therefore also containing 3 observations.

For both cases, we tested the cycling windows strategies in order to adapt the temporal variation of parameters. Here the parameters values change at the beginning of every assimilation window. In terms of the iterative scheme, we found it necessary to introduce similar parameter evolution model between two outer loops. For the 1D case, we firstly carried out a run of 30 assimilation windows in order to show the stability of the scheme. Then we adopted 6 windows along with a forecast stage whose duration equals 3 assimilation windows to show the details of error growth in the first few cycles. Similarly, 5 windows are chosen for the 2D case. Note that in these synthetic cases, the coarse grid nodes of the background state coincide with those of fine synthetic observations. Finally, we applied this stochastic model to those Kinectcaptured images. Kinect is a depth sensor developed by Microsoft as the motion sensing input device for game consoles. Those images related to the work of [START_REF] Combes | A particle filter to reconstruct a free-surface flow from a depth camera[END_REF] are captured by a Kinect sensor to reflect the evolution of a unidirectional wave generated by an initial free surface height difference ∆h = 1 cm on a grid of 222 × 88 pixels. The flow is bounded in a rectangle flat tank of size Lx × Ly = 250 cm × 100 cm. The flow surface was located between 680 and 780 mm from the device. When the attenuation coefficient of the liquid is larger than 113 m -1 , the Kinect sensor displays a mean measurement error of 0.5 mm with standard deviations of about 0.5 mm for both flat and sinus-like surfaces. The sensor captures successfully sinus-like varying elevations with spatial periods smaller than 20 mm and amplitudes smaller than 2 mm. However, these raw data suffer large areas of missing information. Therefore, the same image pre-processing as in [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF] has been done in order to fix the issues related to the observation errors' spatial distribution and the background state, etc. For the observation errors: in terms of a point located in the unobserved region, the observation error is set to be dependent on the distance. The longer the distance, the larger the error. The observation error is bounded by a maximal value of 60% of the height difference ∆h. Within the observed region, we set the observation error homogeneously to the instrument error σo = 0.5%∆h. The simulation is executed on a 124 × 49 rectangle grids mapped on the same domain. In this case, the initial background was completely unknown, hence it was set to a filtered observation with interpolated values on the missing data regions at the initial time. Since the initial background quality is bad, the system is designed to assimilate repeatedly the first few observations in order to obtain a faster 'spin-up'. So we chose five overlapping assimilation windows over nine observations with each window containing 5 observations. In all cases, the observation error covariance matrix R is assumed to be diagonal and all the observations are local in the sense that all the observations are taken at a certain spatial location. We employed different strategies to construct the ensemble in the synthetic case and in the real case. For the synthetic cases, as there is no explicit initial error associated with the background, the ensemble was initialized as a function of the parameters. And the prior parameters are determined from an empirical estimation of the variance tensor on local spatiotemporal regions, C(x k , t), centered on the coarse grid points indexed by k

a(ti, x k ) = C δt n 2 -1 i,l∈C (w(ti, x l ) -w(t, xC) ) (w(ti, x l ) -w(t, xC) ) T , (41)
where δt is the characteristic time step used to regularize the unit of a and • denotes the spatial-temporal mean. The error uncertainty is computed by equation (37). For the real cases, parameters are perturbed in addition to the initial state Gaussian perturbations to initialize the ensemble. Since no velocity observation is available at the initial time, the parameter initial observation is given by the initial ensemble variance:

a = Cδt N -1 N j=1 [w j (t) -w(t)][w j (t) -w(t)] T . (42)

Comparison criterion

With a reference state available in both synthetic cases, we used the Root Mean Square Error (RMSE) to measure the differences between the values predicted by a model or an estimator and the reference value. Such RMSE is defined through L2 norm:

RMSE = ||x target -x ref ||2, (43)
In the 2D synthetic case, we check the lengthwise velocity besides the surface height. The lengthwise velocity is the primary velocity component along Lx. In the real image case, since there is no ground truth, we compare the background/analysis at discrete time levels against the observations as well as the surface height profile obtained by different methods.

Results and discussions

1D case results, synthetic observations

As for the 1D experiment, the spatial mean of the correlation coefficient calculated from equation (39) are cor(h, axx) = 0.1649 and cor(u, axx) = 0.1804. Such an order of magnitude of the correlation coefficient ensures the identifiability of a through the state variables.

A 30 pure assimilation cycles' test has been performed for the 4DEnVar schemes with perturbed observations and localized covariance. In order to evaluate the robustness of the proposed scheme and the effectiveness of the stochastic model, we carried out one test with the model under uncertainty (with suffix 'stochastic-model') and two tests on the "standard-model" in which no subgrid parameterization is introduced and the ensemble is generated by perturbing the initial condition state with the uncertainties computed a few time steps after the starting of the experiment. The two tests on the "standard-model" differs on whether or not the ensemble multiplicative inflation technique discussed in section 3.2.3. is used or not. The figure (1) shows the RMSE evolution in the course of the 30 cycles (here we use PO prefix denoting perturbed observations and LC prefix denoting localized covariance). The results are clearly in favor of our proposed inflation-free uncertainty model. The ensemble inflation technique improved to a certain extent the result but further improvement can hardly be obtained by simply enlarging inflation coefficient γ since this causes numerical instability for the standard model.

The RMSE curves plotted in figure (2) represent the results obtained by the different assimilation configurations mentioned in section 3.2. After a 6 windows cycling of assimilation procedure, we run a forecast stage, whose duration equals the length of 3 assimilation windows. We have discussed in section 3.2.4. about the two different sources of the parameter estimation effect. The purpose of a prevision stage is to study only the effect upon the numerical integration using the prescribed stochastic model. In terms of the performance with the model under uncertainty, we have found that with no damping in parameter evolution model, insufficient number of ensemble members (< 50) lead to model blow-up characterized by unrealistic model variables' values. This is not shown in figure 2 since it makes the figure scale unreadable. By introducing adequate damping factor β ∈ [0.2, 0.5] which attenuates the change rate of variance tensor a, the black curve displayed in figure 2 only experiences a strong divergence in the first cycle, while this behavior is dampened at later cycles. The analysis trajectories across the assimilation stage show that increasing the ensemble size allows achieving a considerable improvement in the analysis quality in terms of the initial condition and indirectly of the parameters. We also find that the localization technique with optimally tuned localization scales enables improving the result for ensemble of small dimension. Note that compared to a localization procedure, increasing the ensemble size do not improve significantly the analysis performance. In fact the RMSE curve using 512 sampling numbers and no localization scheme was very close to the one with only 32 members. Therefore it has not been included in the figure for sake of readability.

In terms of the forecast stage (starting from 15000 in figure 2), the forecast trajectories associated with better analysis states from the last cycle yield less error and distinct error growth schematic compared to the background. Conversely, the error evolutions with the standard model (with or without inflation) coincide closely with the one of the background. This reflected the power of the subgrid modeling in producing more relevant analysis state. Nevertheless, a similar divergence can be found in all the forecast trajectories. This behavior indicates that the effect of the parameters here is only slightly significant with respect to forecast. This seems to indicate that the strategy consisting to keep the variance tensor constant along the assimilation window (and therefore for the whole forecast) does not provide a slower slope of the error growth at long terms. Its assimilation performance enables nevertheless to get forecast results of much better quality. Choosing a time-dependent variance tensor may lead to better result for long-terms forecasting.

Eventually we conjectured that the effectiveness of the corresponding methods is actually related to the ensemble spread. The ensemble spread curves shown in figure (3) correspond to three identical experiment setup (No localization scheme with same sampling numbers) but with the state driven by the stochastic model and the standard model with or without inflation scheme respectively. Note that the first cycle is not shown here since during this cycle the ensemble spread is dominated by arbitrary noise. We observed that the ensemble spread for stochastic model is significantly larger than the spreads associated with the standard model (with or without inflation). The spread curve of the standard model without inflation not only decreases faster but also collapses of nearly 100% after the fourth cycle. The inflation technique do increase the ensemble spread but its power is intrinsically limited by the model stability. On the other hand, we observed that the ensemble spread of the stochastic model also decreases but is still maintained at a rather high level across the sliding windows. A relatively high-level ensemble spread explains why the forecast state is corrected effectively at the beginning of each cycle in terms of the RMSE curve (e.g. the black dotted line in figure 2); similarly, a collapsed ensemble related to the standard model loses gradually the ability of improving the forecast state.

We have also conducted a study on the different terms composing the subgrid model in Eqs.( 17). The configuration of different cases are shown in table (1). The RMSE evolution along 3 cycles are illustrated in figure (4). We can conclude that the stochastic shallow water model works equally well with the case aI (complete terms presented) and the case aII (no term associated with the gradient of a in the momentum equation). Case aIII does not take account of the terms associated with the gradient of a in both mass and momentum conservation equations. The initial analysis resulting from this case is relatively bad Table 1. Different combination of the terms composing the subgrid model.

Momentum equation Continuity Equation

a I -1 2 axx∂x(hux) -1 2 axx∂x(uhx) -1 2 axx∂xxh -∂xaxx∂xh -∂xaxx∂x(hu) a II -1 2 axx∂x(hux) -1 2 axx∂x(uhx) -1 2 axx∂xxh -∂xaxx∂xh a III -1 2 axx∂x(hux) -1 2 axx∂x(uhx) -1 2 axx∂xxh a IV -1 2 axx∂x(hux) none
compared to the previous two cases, which indicates the relevance of the missing modeling part of the continuity equation in generating larger ensemble spread. Case aIV assumes a similar form as the Smagorinsky-like eddy viscosity model described in the previous section without any subgrid term in the continuity equation. Besides the significantly lower quality of initial analysis state, the evolution of RMSE diverges faster and to a greater extent compared to other 3 counterparts. This test suggests that the subgrid modeling term presented in the continuity equation is crucial to keep an accurate representation of the small-scale uncertainty constituting a more balanced ensemble spread between the height and the velocity fields. Note that case aI is the default case for all the other tests conducted in this article.

2D case results, synthetic observations and real data

Under the 2D experimental setup, due to the sizable degrees of freedom of the system, the localization-free paradigm basically does not work anymore. So here we focus on the comparison of two localization techniques: localized covariance (LC) and local analysis (LA). [START_REF] Sakov | Relation between two common localisation methods for the enkf[END_REF] have shown the similarities of the two types of localization scheme within the framework of EnKF under some restricted conditions. In our practice, other factors interfere so that the similarity of the two different localization scheme may break down. One of factors is the effectiveness of the minimization scheme since our localization techniques are constructed on a variational basis. The background error statistics as well as the relevance of the initial ensemble also contribute to the inequality. Another important factor is the ensemble update schemes used in combination with LC or LA technique. We use perturbed observation scheme (PO) with LC and ensemble transform scheme (ET) in conjunction with LA. In order to separate the two factors (the ensemble update scheme and the localization scheme), we employed another scheme (section 3.2.1.3.) using global covariance localization and a deterministic ensemble update approach (denoted as DE-LC) without perturbing the observations. All those have been run for the 2D stochastic large-scale dynamics described eqs. (15a,15b).

To be able to diagnose the performances, the first step is to study the sensitivities of localization scales. The localization scale is termed as cut-off distance in terms of LC and local domain size in terms of LA respectively. Figure 5 shows the mean of RMSE of height per cycle against different localization scales. We observed that the LC approach tolerates a larger scale while the LA approach clearly is more in favor of a smaller scale. Although the average RMSE monotonically increased in function of the local domain size when using LA, the local domain size smaller than 10%Lx contains only several grid points which results in numerical instability. Eventually, we chose 20%Lx as the cutoff distance for LC and 10%Lx as the local domain size for LA.

The RMSE results regarding the 2D Synthetic Results are shown in figure (6). In terms of the group of three LC scheme, the error evolution of the localized covariance (LC) cases with stochastic model indicated consistent sensible results. The result of a PO-LC test using the standard shallow water model and multiplicative inflation technique shows relatively higher error and a diverged velocity RMSE from the 2nd cycle compared to those associated to the stochastic model using LC scheme. Between two ensemble update schemes using LC and the stochas- 1 of stochastic shallow water model ( 17). The meaning of different a is shown in table 1. tic subgrid modeling, the DE scheme yield comparable result in height field but slightly worse result in velocity field. In terms of ET-LA scheme, as shown by the black solid line in figure (6), a similar catastrophic filter divergence to the 1D case can be found in the first cycle with underdamped parameter model. We can conclude that the damping parameter model is essential to the first cycle for the parameters estimated in the first cycle are badly constrained. After optimally tuning the value of β, we can see that the analysis state for LA approach (the black dashed line in figure ( 6)) is corrected at the 1st cycle and its trajectory does not diverge anymore. However, the analysis across later cycles is characterized by a monotonically increasing error caused by the collapsing of ensemble spread. Because when the ensemble spread is too small, the state variables will conform to the background and the observation will barely exert any effect. Parameter ensemble collapsing is a common problem encountered by many studies of parameter estimation in an augmented state vector framework [START_REF] Aksoy | Ensemblebased simultaneous state and parameter estimation in a twodimensional sea-breeze model[END_REF]Simon et al., 2015). Regarding our test, firstly, after the 1st cycle, the estimated parameters have to be severely tampered before they can be used to evolve the analysis ensemble to the beginning of the 2nd cycle (therefore formulating the forecast ensemble). This damping process is a source of loss of signal in the ensemble parameter fields. Secondly, during the assimilation process, the variances of both the state variables and the parameters are naturally decreased. We concluded that the reason why the ensemble fields deteriorates in LA is linked to the form of calculating the ensemble anomaly fields in local domains. Since the noise that we are trying to characterize in this paper is mainly due to the non-consistent spatial and temporal sampling numerical steps. Consequently, using a transformation matrix of small dimension can lead to excessive small-scale variations and eliminate the longer correlations since it restricts the directions of the variations of increments in the ensemble space. To verify this, we plot the ensemble anomalies' PSD (power spectral density) in terms of height with respect to the forecast ensemble, the analysis ensemble of LC and LA respectively for the 1st cycle in figure 7. We observed that after integrating the stochastic model for some time, the forecast ensemble is clearly dominated by few leading modes. The analysis ensemble yielded by LC preserves considerably the spectral form while the magnitude associated with each mode decreases. The result clearly reveals why the LC approach is more robust even with a relative small ensemble size. However, the computational cost associated with LC approach of 64 members ensemble is already considerably high. The analysis ensemble produced by LA, on the contrary, spreads the energy over much more modes. These small-scales energy variations are then very likely smoothed out by the diffusive subgrid model.

In order to make LA work better, we adopt the relaxation methods discussed in section 3.2.3. in which α is tuned as 0.9 for state ensemble. As shown by the black dotted line in figure 6, the LA approach with relatively large ensemble generates sensible analysis comparable with the LC experiments.

One possible fix besides this relaxation method is to use leftmultiplication ensemble transformation schemes which is in the state vector space. Note that this scheme is both compatible with LC and LA. However, the left-multiplied form will require the minimization process to be done in state vector space which is more computational demanding. Another possible approach would be to introduce a static full-rank B and to apply a 3DVar-FGAT in order to add more directions in the incremental field which is not solely covered by the ensemble B as suggested in [START_REF] Gustafsson | Four-dimensional ensemble variational (4D-En-Var) data assimilation for the HIgh Resolution Limited Area Model (HIRLAM)[END_REF]. This will in fact add an hybrid taste in our EnVar scheme, which constitutes an interesting perspective for advanced study. From the perspective of image data, it is also beneficial to introduce correlation into the observation error covariance matrix [START_REF] Chabot | Accounting for observation errors in image data assimilation[END_REF]. This remains another interesting outlook for future work.

In terms of the real image data depicting the free surface height of a laboratory tank, figure (8a) shows the evolution of the mean surface height of the wave crest region using different methods. We particularly chose to focus on this region as the flow is dominated by a single wave. Figure (8b) shows the RMSE between background/analysis state and the observation at discrete observation times. Note the analysis trajectories denoted by 4DVar, 4DEnVar-PO-LC and 4DEnVar-ET-LA are the same as in figure 8 in [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF], which are produced by the standard model without any subgrid parameterizations. No ensemble inflation technique is used in this test. Basically, the group of 4DEnVar outperforms the 4DVar in terms of the surface height at the beginning of the assimilation window. We observed an underestimation of the surface height for 4DVar in figure (8a), while the group of 4DEnVar follows quite well the first 2 observations. The analysis trajectories generated by EnVar schemes with stochastic models (with Sto suffix) show a slight improvement between the 3rd and the 6th observations. Their trajectories started to attenuate in the last cycle compare to their counterpart using standard model. This local phenomenon is again validated in the global RMSE results of figure (8b). We can observe that the RMSE values regarding the stochastic models between the 3rd and the 6th observations are consistent with the 4DVar case. Lower RMSE associated with stochastic models are also found in the last cycle (black symbols versus cyan symbols). We recalled in this experiment the initial background state is unknown, which makes the numerical performance of the model highly unstable. Hence, the parameterization procedure can only tolerate a relatively low magnitude of the variance tensor, which cannot induce significant ensemble perturbations. All in all, the stochastic version has been found more stable than the other schemes.

We have also plotted the ensemble spread versus time in figure (9). As in the synthetic case, the ensemble spread associated with subgrid model is larger than those of the standard model for the whole assimilation windows. However, the difference of ensemble spread between the two models is less evident compared to the synthetic case, which indicate the perturbed initial states are mainly responsible for generating relevant background error estimations especially for the later cycles. This behavior also sheds light on the divergence of the analysis trajectory.

Concluding remarks

In this paper, we started by presenting shortly the principle of the stochastic fluid dynamics proposed in [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] offers us a way to reinterpret the so called subgrid stress models used in CFD for many years. In this framework, the subgrid models are theoretically sound as they are directly related to the small-scale uncertainty components and inferred from physical conservation laws. They are therefore not prone to the assumptions underlying the eddy viscosity models (to list few of them: filter permutation with differentiation, specific shape of the dissipation, isotropic turbulence, diffusion alignment with the rate of strain proper directions).

Then we have focused on the numerical evaluation of the stochastic model with a coherent scheme 4DEnVar proposed in [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF]. The idea is this: since we argue that the subgrid-scale stress term attached to the uncertainties plays an essential role in formulating the effects of the uncertainties on the resolved scales, we can indeed try to estimate these effects through high-resolution data by a parametrization of the associated terms. In our 4DEnVar scheme, the unknown augmented initial vector is adjusted to fit the 4D observations simultaneously. Besides, due to its cycling and iterative ensemble update strategy, we are able to fully utilize the 4D Ensemble in order to meet the needs of dealing with bad background and incomplete observations.

The numerical results have been encouraging. Using the stochastic shallow water model with synthetic data sets, the 4DEnVar scheme is tested to compensate the gap between the coarse resolution background state and fine resolution observations. The results clearly indicated the pertinence of the quantification of the unknown small scale physical processes or numerical artifacts. We have a much better analysis state in global quality in terms of spatial errors; and the error growth tendency of the analysis trajectory with the estimated subgrid-scale stress parameters is distinct from other cases driven by the standard model (with/without inflation), which suggests that more relevant information form the observation can be extracted and contribute to the formulation of a more correct and coherent analysis state. Such an improvement is largely due to two related subjects. The stochastic model introduced more relevant ensemble spread growth and avoided the common catastrophe of ensemble collapse due to the insufficient ensemble spread. The ensemble-based DA scheme, therefore benefited from these 4D ensemble in formulating the necessary components (innovation vector, ensemble anomaly matrix, etc.), and consequently gen-erated more accurate analysis states and parameter fields for the stochastic model. It is important to note that the subgrid-scale stress term appeared in the stochastic model, although only labeled as resolution-related, can actually embrace broader types of errors as long as these errors can be inferred from the discrepancies between the dynamics and the high-resolution data.

The numerical results are also discussed in detail with two localization techniques. In ensemble-based methods, localization is useful not only in removing the unrealistic correlation, but also in enlarging the control space to a space with greater degree of freedom. With 2D synthetic or real case, using localization is obligatory. The localized covariance approach is proved to be more robust in the sense that the posterior ensemble produced by such scheme not only kept its magnitude, but also maintained to a large extent its energy spectral form. The spectral form is important in terms of our experimental setups since it prevented the small scale energy variations from diffusing by the stochastic model, which is exactly the case with the local analysis scheme. The results revealed by the real data case is less convincing due to the numerical instabilities associated with the variance tensor. The analysis trajectory only showed a slight improvement with a more stable behavior in the analysis. More work needs to be done in order to provide a better initial background, which eventually allows generating a more pertinent prior ensemble. Employing the full stochastic model with the noise term constitutes an interesting perspective for this path.

Our results suggested that the two localization techniques both suffered some limitations in their current forms. In the localized covariance approach, the posterior ensemble can only be sought through an ensemble analysis, which is computational demanding. The local analysis, although much cheaper and easy to implement, has shown some difficulties in bridging the gap between the local and the global representations. Following works will investigate new localization scheme that do not have the above restrictions. One way to improve the posterior ensemble efficiency and accuracy is to construct the estimation scheme on adjoint models. By introducing the adjoint models, the optimization is done in the state space rather than the ensemble space. Although maintaining and computing an adjoint model is cumbersome and time-consuming, using a state space formulation allows to express the left transformation of the prior ensemble to the posterior ensemble [START_REF] Bocquet | Localization and the iterative ensemble kalman smoother[END_REF]. In this form, the ensemble transformation takes the form of EAKF [START_REF] Anderson | An ensemble adjustment Kalman filter for data assimilation[END_REF][START_REF] Sakov | Relation between two common localisation methods for the enkf[END_REF]. In the future we are planning to adopt similar strategy with complex data sets. Since the scheme is formulated on variational base, the transition is straightforward. Adding a static full-rank B and applying a 3DVar-FGAT [START_REF] Gustafsson | Four-dimensional ensemble variational (4D-En-Var) data assimilation for the HIgh Resolution Limited Area Model (HIRLAM)[END_REF] can also be beneficial to the posterior ensemble quality as more directions are permitted to grow in the incremental field. This path constitutes another promising direction to follow. figures the small-scale velocity component associated to a much thinner time-scale. This component is modeled as a deltacorrelated random field in time as it represents a highly non regular process at the resolved time scale. In this work, we assume that this small-scale random component is incompressible and therefore associated to a divergence-free diffusion tensor: ∇ • σ = 0. This assumption, which obviously does not prevent the drift component (and therefore the whole field) to be compressible, leads to much simpler modeling and remains realistic for the models considered in this study.

Let us consider now a spatially regular enough scalar function φ of compact support, transported by the stochastic flow (A1) and that vanishes outside volume V(t) and on its boundary ∂V(t). As this function is assumed to be transported by the stochastic flow, it constitutes a stochastic process defined from its initial time value g: φ(Xt, t) = g(x0).

We will assume that both functions have bounded spatial gradients. Besides, the initial function g : Ω → R vanishes outside the initial volume V(t0) and on its boundary. Let us point out that in this construction, function φ cannot be a deterministic function. As a matter of fact, if it was the case, its differential would be given by a standard Ito formula ∂tφ + ∇φ • w + 1 2 i,j d X i t , X j t ∂ 2 φ ∂xi∂xj dt+ ∇φ • σd Bt = 0, (A5) which here cancels as φ is transported by the flow. A separation between the slow deterministic terms and the fast Brownian term would yield a null uncertainty or a specific uncertainty orthogonal to ∇φ. This would lead us back to the deterministic case, which is not the general goal followed here. As φ is a random function, the differential of φ(Xt, t) involves the composition of two stochastic processes. Its evaluation requires the use of a generalized Ito formula usually referred in the literature to as the Ito-Wentzell formula (see theorem 3.3.1, [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF]. This extended Ito formula incorporates in the same way as the classical Ito formula a quadratic variation terms related to process Xt but also co-variation terms between Xt and the gradient of the random function φt. Its expression is in our case given by dφ(t, Xt) = dtφ + ∇φ In these expressions f (resp. σ) designates the kernel associated to operator f (resp. σ). Now, identifying first the Brownian term and next the deterministic terms of equations ( A6) and (A7), we infer f (x, y, t) T = -∇φ(x, t) T σ(x, y, t), v(x, t) = -∇φ • w + With the complete expression of L * (remarking that the second right-hand term of A12 is self-adjoint) and if φ(x, t) → 1I V(t)/∂V(t) , where 1I stands for the characteristic function, we get the sought form of this differential:

V(t) dtq + ∇• qw - 1 2 ∇• (qa) T dt + ∇qσdBt dx,

APPENDIX B: Quadratic variation and covariation

We recall here the notion of quadratic variation and covariation, which play a central role in stochastic calculus. We will here restrict ourselves to standard Ito processes. Quadratic variation and co-variation correspond respectively to the variance and covariance of the process increments along time. The quadratic co-variation process denoted as X, Y t, (respectively the quadratic variation for Y = X) is defined as the limit in probability over a partition {t1, . . . , tn} of [0, t] with t1 < t2 < • • • < tn, and a partition spacing δti = ti -ti-1, noted as |δt|n = max For Brownian motion these covariations can be easily computed and are given by the following rules B, B t = t, B, h t = h, B t = h, h t = 0, where h is a deterministic function and B a scalar Brownian motion.

Fig. 1 .

 1 Fig. 1. RMSE comparison in terms of free surface height (a) and velocity (b) between various configurations. PO: Perturbed observation. LC: Localized covariance approach. Standard model: no subgrid tensor parameterization. N: number of ensemble member.

Fig. 2 .

 2 Fig. 2. RMSE comparison in terms of free surface height (a) and velocity (b) between various configurations of 4DEnVar. PO: Perturbed observation. LC: Localized covariance approach. Standard model: no subgrid tensor parameterization. N: number of ensemble member. 6 DA cycle phase: t ∈ [6000, 16800)s. Forecast phase: t ∈ [15000, 22200]s

Fig. 3 .

 3 Fig. 3. Ensemble spread in terms of free surface height (a) and velocity (b): standard model (circle) with ensemble in function of initial states; standard model (triangle) with ensemble in function of initial states but with multiplicative inflation technique; subgrid model (square) with ensemble in function of initial parameters.

Fig. 4 .

 4 Fig. 4. RMSE comparison in terms of free surface height (a) and velocity (b) between various subgrid model configurations in table1of stochastic shallow water model (17). The meaning of different a is shown in table 1.

Fig. 5 .

 5 Fig. 5. The RMSE of height comparison in the 1st and 2nd cycle as a function of Cut-off distance (LC, localized covariance) or local domain size (LA, local analysis).

Fig. 6 .

 6 Fig. 6. RMSE comparison in terms of free surface height (a) and lengthwise velocity (b) between various configurations of 4DEnVar. ET, ensemble transform. PO, perturbed observation. DE, deterministic ensemble. LC, localized covariance. LA, local analysis. ET-LA-N128-Stochastic (solid line) case with β = 1 shows ensemble 'blow up'; ET-LA-N128-Stochastic (dashed line) with β = 0.5 shows ensemble collapse; ET-LA-N512-Stochastic (dotted line) used ensemble relaxation inflation α = 0.9 technique for state ensemble; PO-LC-N64-Standard used standard model and multiplicative inflation technique for state ensemble.

Fig. 7 .

 7 Fig. 7. Power spectral density of ensemble anomalies for the forecast ensemble (a), the analysis ensemble by LC (b) and LA (c).

Fig. 8 .

 8 Fig. 8. Comparison of various configurations of 4DEnVar w.r.t. (a) Mean surface height of the wave crest region as a function of time and (b) The RMSE against observation at observation time. 4DVar, 4DEnVar-PO-LC and 4DEnVar-ET-LA are tested with the standard model; 4DEnVar-PO-LC-Sto and 4DEnVar-ET-LA-Sto are tested with the stochastic model. LC, localized covariance. PO, perturbed observation. LA, local analysis. ET, ensemble transformation.

Fig. 9 .

 9 Fig. 9. Ensemble Spread in terms of free surface height (a) and lengthwise velocity (b) for: standard model cycle 1 and cycle 5 (circle with solid/dashed line); stochastic model cycle 1 and cycle 5 (square with solid/dashed line).

=V

  -∇φ • w + ∇φ • (∇• a) a∇φ). (A12)This differential at a fixed point, x, defines the equivalent in the deterministic case of the material derivative of a function transported by the flow. The differential of the integral over a material volume of the product qφ is given by d + dt q, φ dx, where the first line comes from φ(t, x) = 0 if x ∈ Ω\V(t) and the second one from the integration by part formula of two Ito processes. Hence, from (A11) this differential isΩ dtqφ + qLφ + ∇φ • a∇q dt dx -Ω q∇φ • σdBt.Introducing L * the (formal) adjoint of the operator L in the space L 2 (Ω) with Dirichlet boundary conditions, this expression can be written as Ω dtq + L * q -∇• (a∇q) dt +∇• qσdBt φ dx.

i

  δti and such that |δt|n → 0 when n → ∞: )-X(ti) Y (ti+1)-Y (ti)T .

  , which

		0.022							Background
										4DVar
		0.0215							4DEnVar¡PO¡LC 4DEnVar¡ET¡LA
										4DEnVar¡PO¡LC¡Sto
		0.021							4DEnVar¡ET¡LA¡Sto Observations
	(m)	0.0205
	h
		0.02
		0.0195
						Cycle 1
		0.019				Cycle 2	Cycle 3
									Cycle 4
										Cycle 5
		0.03 0.0185	0.06	0.09	0.12	0.15	0.18	0.21	0.24	0.27
								t(s)
								(a)
			2.8	x 10 -3
			2.6
		Surface height error(m)	1.4 1.6 1.8 2 2.2 2.4
			1.2
			1
			0.03 0.8	0.06	0.09	0.12	0.15	0.18	0.21	0.24	0.27
								t(s)

  • dXt +It can be immediately checked that for a deterministic function, the standard Ito formula is recovered since the co-variations terms between Xt and ∇φt cancel.It follows from A6 that for a fixed grid point, function φ(x, t) is solution of a stochastic differential equation of the form dtφ(x, t) = v(x, t)dt + f (x, t) • dBt,

									(A7)
	where the Brownian term must compensate the Brownian term
	of (A6). The quadratic variation term involved in (A6) is given
	through (appendix B) as
		d X i t , X j t = a ij (x, t)		(A8)
				=		σ ik (x, t)σ kj (x, t)
					k
				=			σik (x, y, t)σ kj (y, x, t),	(A9)
					k	Ω
	and the covariation term reads
	d	∂φt ∂xi	, X i t	=	j	Ω	∂ f j ∂xi	(x, y, t)σ
									i	d	∂φ ∂xi	, X i t dt+
									1 2 i,j	d X i t , X j t	∂ 2 φ ∂xi∂xj	dt
									= 0.	(A6)

ij (y, x, t). (A10)

Note that the relationship between the Hessian and the covariance matrix holds rigorously in a linear sense; in a nonlinear scenario, this relationship is only an approximation. As the minimization algorithm LBFGS relies on an approximation of the inverse Hessian matrix H -1 , we can use this byproduct to evaluate equation (30).

APPENDIX A: Stochastic Reynolds transport theorem

We derive in this appendix the expression of the Reynolds transport theorem for the drift plus noise decomposition: In a Lagrangian stochastic picture, the infinitesimal displacement associated to a trajectory Xt of a particle is noted:

The random field involved in this equation is defined for all points of the fluid domain, Ω through the kernel σ(., ., t) associated to the diffusion operator σ(., t)

This operator is assumed to be of finite norm for any orthonormal basis of the associated Hilbert space H d and to have a null boundary condition on the domain frontier:

Function Bt denotes a d-dimensional Brownian function 2 and the resulting d-dimensional random field, σ(x, t)dBt ∈ H d , is a centered vectorial Gaussian function correlated in space and uncorrelated in time with covariance tensor:

We observe that those random fields have a (mean) bounded norm: E σdBt 2 = tr Q < ∞, where the trace of the covariance tensor is given for any basis {e k | k ∈ N} of H d as tr Q = k (e k , Qe k )H. In the following we will note the diagonal of the covariance tensor as: a(x)dt = Q(x, x). Tensor, a, that will be referred to as the variance tensor

is by definition a symmetric positive definite matrix at all spatial points, x. This quantity corresponds to the time derivative of the so-called quadratic variation process:

The notation f, g stands for the quadratic cross-variation process of f and g. This central object of Stochastic Calculus, can be interpreted as the covariance along time of the increments of f and g. The quadratic variation process is briefly presented in Appendix B. The drift term, w, of Lagrangian expression (A1), represents

2 Formally it is a cylindrical I d -Wiener process (see Da Prato and Zabczyk (1992) for more information on infinite dimensional Wiener process and cylindrical I d -Wiener process).

the "smooth" -differentiable -part of the flow whereas the random part, Ẇ = σ dBt dt , (A4)