G Moldovan

G Lehnasch

L Cordier

M Meldi

A multigrid/ensemble Kalman filter strategy for assimilation of unsteady flows

Keywords: Kalman Filter, Data Assimilation, compressible flows

A sequential estimator based on the Ensemble Kalman Filter for Data Assimilation of fluid flows is presented in this research work. The main feature of this estimator is that the Kalman filter update, which relies on the determination of the Kalman gain, is performed exploiting the algorithmic features of the numerical solver employed as a model. More precisely, the multilevel resolution associated with the multigrid iterative approach for time advancement is used to generate several low-resolution numerical simulations. These results are used as ensemble members to determine the correction via Kalman filter, which is then projected on the high-resolution grid to correct a single simulation which corresponds to the numerical model. The assessment of the method is performed via the analysis of one-dimensional and two-dimensional test cases, using different dynamic equations. The results show an efficient trade-off in terms of accuracy and computational costs required. In addition, a physical regularization of the flow, which is not granted by classical KF approaches, is naturally obtained owing to the multigrid iterative calculations. The algorithm is also well suited for the analysis of unsteady phenomena and, in particular, for potential application to in-streaming Data Assimilation techniques.

Introduction

In Computational Fluid Dynamics (CFD), newly developed numerical methods are generally assessed in terms of accuracy via comparison with experimental data [START_REF] Oberkampf | Verification and validation in computational fluid dynamics[END_REF]. In practice, this validation step is far from being trivial since many sources of error are inevitably introduced in the simulations. First, the partial differential equations used to derive the numerical scheme may be restricted to oversimplified physical models, such as the Boussinesq approximation applied in thermal convection or the incompressibility condition. Second, the discretization process and the use of iterative numerical methods introduce computational errors in the representation of the flow features [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. Third, boundary and initial conditions are usually very sophisticated in complex applications but detailed a priori knowledge is insufficiently available. Last, for very high Reynolds number configurations, turbulence/subgrid-scale modelling must be included in order to reduce the required computational costs [START_REF] Pope | Turbulent flows[END_REF]. All of these sources of error exhibit complex interactions owing to the non-linear nature of the dynamical models used in numerical application, such as the Navier-Stokes equations.

The experimental results are also affected by uncertainties and biases. In many cases, the set-up of the problem can be controlled up to a finite precision (alignment between the flow and the wind tunnel or immersed bodies, mass flow rate, . . . ). This kind of uncertainty, which is clearly affecting every physical system but cannot be exactly quantified, is usually referred to as epistemic uncertainty. In addition, experimental results are also affected by the precision (or bias) of the acquisition and measurement system. Thus, the main difficulty in the comparison between numerical and experimental results is understanding how much of the differences observed is due to an actual lack of precision of the numerical model, and how much is instead associated to differences in the set-up of investigation.

One possible strategy to obtain an optimal combination of the knowledge coming from simulations and experiments is to derive a state estimation which complies with both sources of information. The degree of precision of such es-timation is connected with the confidence in the sources of information. This has the advantage of naturally incorporating uncertainty and bias present in the sources of information in the analysis. Tools providing such state estimation are usually included in different disciplines, control theory for state observers [START_REF] Zhou | Robust and optimal control[END_REF], Data Assimilation (DA) [START_REF] Asch | Data Assimilation: methods, algorithms, and applications[END_REF][START_REF] Evensen | Data Assimilation: The Ensemble Kalman Filter[END_REF][START_REF] Daley | Atmospheric Data Analysis[END_REF] for weather prediction, ocean modelling and more recently mechanical engineering problems. Essentially, DA methods combine information issued from two sources: i) a model, which provides a dynamical description of the phenomenon in the physical domain, ii) a set of observations, which are usually sparse and/or local in time and space. These methods are classified in different families according to the way the state estimation is performed. One of the classical criterion of classification deals with the operative strategy used to derive the state estimation. Variational approaches resolve a constrained optimization problem over the parametric space characterizing the model (usually coefficients defining boundary conditions or physical models).

The solution of the variational problem minimizes prescribed error norms so that the assimilated model complies with the observation provided over a specified time window. Methods from this family, which include 3D-Var and 4D-Var, usually exhibit very high accuracy [START_REF] Onder | Optimal control of a transitional jet using a continuous adjoint method[END_REF][START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean flow reconstruction[END_REF][START_REF] Mons | Kriging-enhanced ensemble variational data assimilation for scalar-source identification in turbulent environments[END_REF][START_REF] Chandramouli | 4d large scale variational data assimilation of a turbulent flow with a dynamics error model[END_REF]. However, they are also affected by several drawbacks. First, the formulation of the adjoint problem that is introduced to perform the parametric optimization can be difficult, if not impossible, when automatic differentiation is not employed. Second, the adjoint problem is defined backward in time which may lead to numerical difficulties of resolution related to the unaffordable data storage that is needed, and the amplification of the adjoint solution that frequently happens when multi-scale interactions are dominant [START_REF] Asch | Data Assimilation: methods, algorithms, and applications[END_REF][START_REF] Onder | Optimal control of a transitional jet using a continuous adjoint method[END_REF][START_REF] Sirkes | Finite Difference of Adjoint or Adjoint of Finite Difference?[END_REF]. Third, standard variational data assimilation methods are intrusive from a computational point of view, requiring the development of an adjoint calculation code, or the availability of the source code when the use of automatic differentiation is planned. In most cases, modifications are required in order to maintain the code or extend the applications (change of the definitions of errors, for instance). For commercial software without opensource licence, these modifications are expensive. Non-intrusive methods that require no modification of original calculation code are therefore preferable in data assimilation.

Another family of DA methods is represented by the sequential approaches.

These methods, mostly based on Bayes' theorem, provide a probabilistic description of the state estimation. A well known approach is the Kalman Filter (KF) [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]. Extensions and generalizations of this method have also been developed, such that the Extended Kalman Filter (EKF) [START_REF]An Introduction to the Kalman Filter[END_REF] which is tailored for nonlinear systems, and Ensemble Kalman Filter (EnKF) [START_REF] Evensen | Data Assimilation: The Ensemble Kalman Filter[END_REF]. This class of methods solves the state estimation problem by transporting error covariance matrices of the model and observations. These methods are usually more flexible than variational approaches (no required computation of first order sensitivities), but the time advancement and update of the covariance matrices are prohibitively expensive for large scale problems, encountered in practical applications [START_REF] Rozier | A Reduced-Order Kalman Filter for Data Assimilation in Physical Oceanography[END_REF].

One possible strategy consists in reducing the order of the Kalman filter [START_REF] Suzuki | Reduced-order Kalman-filtered hybrid simulation combining particle tracking velocimetry and direct numerical simulation[END_REF] or filtering the error covariance matrix. Inspired by a domain localization procedure, Meldi & Poux [START_REF] Meldi | A reduced order model based on Kalman Filtering for sequential Data Assimilation of turbulent flows[END_REF][START_REF] Meldi | Augmented Prediction of Turbulent Flows via Sequential Estimators: Sensitivity of State Estimation to Density of Time Sampling for Available Observation[END_REF] proposed a strategy based on an explicit filter of the error covariance matrix. The application of this estimator to different turbulent flows exhibited encouraging results considering the relatively small increase in computational resources. A more popular strategy for data assimilation of engineering applications is the Ensemble Kalman Filter [START_REF] Asch | Data Assimilation: methods, algorithms, and applications[END_REF][START_REF] Evensen | Data Assimilation: The Ensemble Kalman Filter[END_REF][START_REF] Evensen | The ensemble Kalman Filter for combined state and parameter estimation -Monte Carlo techniques for data assimilation in large systems[END_REF], which relies on a Monte Carlo implementation of the Bayesian update problem. The EnKF (and follow-up models) was introduced as an extension of the original Kalman filter made for high-dimensional systems for which transporting the covariance matrix is not computationally feasible. EnKF replaces the covariance matrix by the sample covariance matrix computed from an ensemble of state vectors. The main advantage of EnKF is that advancing a high-dimensional covariance matrix is achieved by simply advancing each member of the ensemble. Several research works have been reported in the literature in recent years for application in fluid mechanics [START_REF] Mons | Reconstruction of unsteady viscous flows using data assimilation schemes[END_REF][START_REF] Xiao | Quantifying and reducing model-form uncertainties in reynolds-averaged navier-stokes simulations: A data-driven, physics informed bayesian approach[END_REF][START_REF] Rochoux | Towards predictive data-driven simulations of wildfire spread -Part I: Reduced-cost Ensemble Kalman Filter based on a Polynomial Chaos surrogate model for parameter estimation[END_REF]. Despite the interest of this non-intrusive technique, and the possibility to perform efficient parametric inference, the computational costs can still be prohibitive for realistic applications. Statistical convergence is usually obtained for a typical ensemble size going from 60 to 100 ensemble members [START_REF] Asch | Data Assimilation: methods, algorithms, and applications[END_REF].

In addition, the state estimation obtained via sequential tools does not necessarily comply with a model solution i.e. the conservativity of the dynamic equations of the model is violated. This aspect is a potentially critical issue in fluid mechanics studies. Violation of conservativity may result in loss of conservation of some physical properties of the flow (such as mass conservation or momentum conservation) as well as in the emergence of non-physical discontinuities in the flow quantities. The aforementioned issues significantly affect the precision of the prediction of the flow and may eventually produce irreversible instabilities in the time advancement of the dynamical model. A number of works in the literature have provided advancement in the form of additional constraints to be included in the state estimation process. Meldi & Poux [START_REF] Meldi | A reduced order model based on Kalman Filtering for sequential Data Assimilation of turbulent flows[END_REF] used a recursive procedure and a Lagrangian multiplier (the pressure field) to impose the zero-divergence condition of the velocity field for incompressible flows. Other proposals deal with imposing hard constraints in the framework of an optimization problem [START_REF] Simon | Kalman filtering with state equality constraints[END_REF], ad-hoc augmented observation [START_REF] Nachi | Kalman filtering in the presence of state space equality constraints[END_REF] and generalized regularization [START_REF] Zhang | Regularized ensemble kalman methods for inverse problems[END_REF]. These approaches are responsible for a significant increase in the computational resources required, which is due to augmentation in size of the state estimation problem or to the optimization process, which usually needs the calculation of gradients of a cost function.

The investigation of physically constrained sequential state estimation is here performed using an advanced estimator strategy, which combines an EnKF approach and a multigrid method. For this reason, we refer to our algorithm as Multigrid Ensemble Kalman Filter (MEnKF). Multigrid methods [START_REF] Brandt | Multi-level adaptive solutions to boundary-value problems[END_REF][START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF] are a family of tools which employ multi-level techniques to obtain the timeadvancement of the flow. In particular, the geometric multigrid [START_REF] Wesseling | Geometric multigrid with applications to computational fluid dynamics[END_REF] uses different levels of the resolution in the computational grid to obtain the final state.

The method here proposed exploits algorithmic features of iterative solvers used in practical CFD applications. The EnKF error covariance matrix reconstruction is performed using information from a number of ensemble members which are generated over a coarse level mesh of a multigrid approach. This procedure is reminiscent of reduced order / multilevel applications of EnKF strategy reported in the literature [START_REF] Hoel | Multilevel ensemble kalman filtering[END_REF][START_REF] Siripatana | Combining ensemble kalman filter and multiresolution analysis for efficient assimilation into adaptive mesh models[END_REF][START_REF] Fossum | Kalman filtering with state equality constraints, Assessment of multilevel ensemble-based data assimilation for reservoir history matching[END_REF][START_REF] Brajard | Combining data assimilation and machine learning to infer unresolved scale parametrisation[END_REF]. However, the state estimation obtained at the coarse level is used to obtain a single solution calculated on a high resolution mesh grid, similarly to the work by Debreu et al. [START_REF] Debreu | Multigrid solvers and multigrid preconditioners for the solution of variational data assimilation problems[END_REF] for variational DA. Because of the algorithmic structure of the problem, all of the simulations on the fine and coarse level can be run simultaneously in parallel calculations, providing a tool able to perform in-streaming DA for unsteady flow problems.

The strategy is tested on different configurations, using several physical models represented by the Burger's equation and the compressible Euler and Navier-Stokes equations for one-dimensional and two-dimensional test cases.

The article is structured as follows. In Sec. 2, the sequential DA procedure is detailed including descriptions of the classical KF and EnKF methods. The numerical discretization and the multigrid strategy are also presented. In Sec 3, the algorithm called Multigrid Ensemble Kalman Filter (MEnKF) is discussed.

In Sec. 4, the MEnKF is used to investigate a one-dimensional case using the Burgers' equation. In Sec. 5, a second one-dimensional case is investigated, but in this case the dynamical model is represented by a Euler equation. In Sec. 6, we investigate the case of the two-dimensional compressible Navier-Stokes equations, namely the spatially evolving mixing layer. Finally, in Sec. 7 concluding remarks are drawn.

Sequential data assimilation in fluid dynamics

In Sec. 2.1, we introduce sequential data assimilation methods starting with the Kalman filter. The objective is to gradually arrive at the Dual ensemble Kalman Filter methodology (Dual EnKF), which is an essential ingredient of the estimator proposed in this work. In Sec. 2.2, the steps that are necessary to transform a general transport equation into a discretized model usable in sequential DA are described. A brief description of the multigrid approach employed is also provided.

Sequential data assimilation

Kalman filter

The role of the Kalman filter (KF) is to provide an estimate of the state of a physical system at time k (x k ), given the initial estimate x 0 , a set of measurements or observations, and the information of a dynamical model (e.g., first principle equations):

x k = M k:k-1 (x k-1 , θ k ) + η k (1)
where M k:k-1 is a non linear function acting as state-transition model and θ k contains the parameters that affect the state-transition. The term η k is associated with uncertainties in the model prediction which, as discussed before, could emerge for example from incomplete knowledge of initial / boundary conditions. In the framework of KF applications, these uncertainties are usually modelled as a zero-mean Gaussian distribution characterized by a variance

Q k , i.e. η k ∼ N (0, Q k ).
Indirect observations of x k are available in the components of the observation vector y o k . These two variables are related by:

y o k = H k (x k ) + o k (2)
where H k is the non linear observation operator which maps the model state space to the observed space. The available measurements are also affected by uncertainties which are assumed to follow a zero-mean Gaussian distribution characterized by a variance R k , i.e. o k ∼ N (0, R k ). The model and observation errors being Gaussian, we can show that the solution is characterized entirely by the first two moments of the state. Following the notation generally used in DA literature, the forecast/analysis states and error covariances are indicated as x f/a k and P f/a k , respectively. The error covariance matrix is defined as

P f/a k = E x f/a k -E(x f/a k ) x f/a k -E(x f/a k ) .
In the case of a linear dynamical model (M k:k-1 ≡ M k:k-1 ) and a linear ob-servation model (H k ≡ H k ), the estimated state is obtained via the following recursive procedure:

1. A predictor (forecast) phase, where the analysed state of the system at a previous time-step is used to obtain an a priori estimation of the state at the current instant. This prediction, which is obtained relying on the model only, is not conditioned by observation at time k:

x f k = M k:k-1 x a k-1 (3)
P f k = M k:k-1 P a k-1 M k:k-1 + Q k (4)
2. An update (analysis) step, where the state estimation is updated accounting for observation at the time k:

K k = P f k H k H k P f k H k + R k -1 (5)
x

a k = x f k + K k y o k -H k x f k (6)
P a k = (I -K k H k ) P f k (7)
The optimal prediction of the state (x a k ) is obtained via the addition to the predictor estimation (x f k ) of a correction term determined via the so called Kalman gain K k . The classical KF algorithm is not suited for direct application to the analysis of complex flows. First of all, KF classical formulation is developed for linear systems. Applications to non-linear systems can be performed using more advanced techniques such as the extended Kalman filter [START_REF]An Introduction to the Kalman Filter[END_REF] or exploiting features of the numerical algorithms used for numerical discretization [START_REF] Meldi | A reduced order model based on Kalman Filtering for sequential Data Assimilation of turbulent flows[END_REF].

The canonical Kalman filter is difficult to implement with realistic engineering models because of i) nonlinearities in the models and observation operators, ii) poorly known error statistics and iii) a prohibitive computational cost. Considering point ii), the matrices Q k and R k are usually unknown and their behaviour must be modelled. One simple, classical simplification is to consider that errors for each component are completely uncorrelated in space and from other components i.e. Q k and R k are considered to be diagonal [START_REF] Stonebridge | Diagonal Approximations to the ObservationError Covariance Matrix in Sea Ice ThicknessData Assimilation[END_REF][START_REF] Brunton | Closed-Loop Turbulence Control: Progress and Challenges[END_REF]. Point iii) is easily explained by the fact that in engineering applications, the dimension of the physical state x k is huge (N ∈ [10 6 , 10 9 ]). By examining the equations of the canonical Kalman filter, it is immediate to see that the procedure relies on the transport of a very large error covariance matrix P k . It is therefore necessary to store it but also to invert very large matrices (see [START_REF] Asch | Data Assimilation: methods, algorithms, and applications[END_REF]). To bypass this computational cost, one popular strategy is to rely on a statistical description obtained via an ensemble of trajectories of the model dynamics.

Ensemble Kalman filter

The Ensemble Kalman Filter (EnKF) [START_REF] Evensen | Data Assimilation: The Ensemble Kalman Filter[END_REF], introduced by Evensen in 1994 [START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using monte-carlo methods to forecast error statistics[END_REF], relies on the estimation of P k by means of an ensemble. More precisely, the error covariance matrix is approximated using a finite ensemble of model states of size N e . If the ensemble members are generated using stochastic Monte-Carlo sampling, the error in the approximation decreases with a rate of 1 √ N e .

In the following, we are describing the stochastic EnKF flavour of EnKF for which random sampling noise is introduced in the analysis [START_REF] Asch | Data Assimilation: methods, algorithms, and applications[END_REF].

Given an ensemble of forecast/analysed states at a certain instant k, the ensemble matrix is defined as:

E E E f/a k = x f/a,(1) k , • • • , x f/a,(Ne) k ∈ R Nx×Ne (8)
To reduce the numerical cost of implementation, the normalized ensemble anomaly matrix is then specified as:

X f/a k = x f/a,(1) k -x f/a k , • • • , x f/a,(Ne) k -x f/a k √ N e -1 ∈ R Nx×Ne , (9)
where the ensemble mean x f/a k is obtained as:

x f/a k = 1 N e Ne i=1 x f/a,(i) k (10)
The error covariance matrix P f/a k can thus be estimated via the informa-tion derived from the ensemble. This estimation, hereafter denoted with the superscript e, can be factorized into:

P f/a,e k = X f/a k X f/a k ∈ R Nx×Nx (11)
The goal of the EnKF is to mimic the BLUE (Best Linear Unbiased Estimator) analysis of the Kalman filter. For this, Burgers et al. [START_REF] Burgers | On the analysis scheme in the ensemble kalman filter[END_REF] showed that the observation must be considered as a random variable with an average corresponding to the observed value and a covariance R k (the so-called data randomization trick). Therefore, given the discrete observation vector y o k ∈ R Ny at an instant k, the ensemble of perturbed observations is defined as:

y o,(i) k = y o k + o,(i) k , with i = 1, • • • , N e and o,(i) k ∼ N (0, R k ). (12)
In a similar way to what we did for the ensemble matrices of the forecast/analysed states, we define the normalized anomaly matrix of the observations errors as

E o k = 1 √ N e -1 o,(1) k -o k , o,(2) k -o k , • • • , o,(Ne) k -o k , ∈ R Ny×Ne ( 13
)
where

o k = 1 N e Ne i=1
o,(i) k

.

The covariance matrix of the measurement error can then be estimated as

R e k = E o k (E o k ) ∈ R Ny×Ny . ( 14
)
By combining the previous results, we obtain (see [START_REF] Asch | Data Assimilation: methods, algorithms, and applications[END_REF]) the standard stochastic EnKF algorithm. The corresponding analysis step consists of updates performed on each of the ensemble members, as given by

x a,(i) k = x f,(i) k + K e k y o,(i) k -H k x f,(i) k (15)
The expression of the Kalman gain is

K e k = X f k Y f k Y f k Y f k + E o k (E o k ) -1 (16)
where

Y f k = H k X f k .
Interestingly, it can be shown that the computation of the tangent linear operator H k can be avoided. The details can be found in [START_REF] Asch | Data Assimilation: methods, algorithms, and applications[END_REF]Sec. 6.3.3]. A version of the Ensemble Kalman filter algorithm using the previously defined anomaly matrices is given in Appendix A.2. This is the version we use in our applications.

State-of-the-art approaches based on the EnKF are arguably the most advanced forms of state estimation available in the field of DA methods. These techniques have been extensively applied in the last decade in meteorology and geoscience [START_REF] Asch | Data Assimilation: methods, algorithms, and applications[END_REF]. Applications in mechanics and engineering are much more recent, despite a rapid increase in the number of applications in the literature. Among those, studies dealing with wildfire propagation [START_REF] Rochoux | Towards predictive data-driven simulations of wildfire spread -Part I: Reduced-cost Ensemble Kalman Filter based on a Polynomial Chaos surrogate model for parameter estimation[END_REF], combustion [START_REF] Labahn | Data assimilation using high-speed measurements and LES to examine local extinction events in turbulent flames[END_REF], turbulence modeling [START_REF] Xiao | Quantifying and reducing model-form uncertainties in reynolds-averaged navier-stokes simulations: A data-driven, physics informed bayesian approach[END_REF] and hybrid variational-EnKF methods [START_REF] Mons | Kriging-enhanced ensemble variational data assimilation for scalar-source identification in turbulent environments[END_REF] have been reported. These applications reinforce the idea that approaches based on EnKF have a high investigative potential despite the highly non-linear, multiscale features of the flows studied by the fluid mechanics community.

Dual Ensemble Kalman filter

In this section, we extend the classical EnKF framework presented in Sec. 2.1.2 by considering the case of a parameterized model such as [START_REF] Oberkampf | Verification and validation in computational fluid dynamics[END_REF]. The objective is to enable the model to generate accurate forecasts. For this, we need to determine good estimates of both model state variables x k and parameters θ k given erroneous observations y o k . One approach is provided by joint estimation where state and parameter vectors are concatenated into a single joint state vector (state augmentation). After [START_REF] Hamid | Dual state-parameter estimation of hydrological models using ensemble kalman filter[END_REF], the drawback of such strategy is that, by increasing the number of unknown model states and parameters, the degree of freedom in the system increases and makes the estimation unstable and in-tractable, especially in the non linear dynamical model. For this reason, we follow the procedure developed by [START_REF] Hamid | Dual state-parameter estimation of hydrological models using ensemble kalman filter[END_REF], called dual estimation. The principle is to apply successively two interactive filters, one for the estimation of the parameters from a guessed state solution, the other for the updating of the state variables from the estimated previous parameters.

In the first step of the algorithm, the ensemble of the analysed parameters is updated following the classical KF equation:

θ a,(i) k = θ f,(i) k + K θ,e k y o,(i) k -y f,(i) k with i = 1, • • • , N e (17)
where

y f,(i) k = H k x f,(i) k .
The Kalman gain responsible for correcting the parameter trajectories in the ensemble is obtained as follows:

K θ,e k = Θ f k Y f k Y f k Y f k + E o k (E o k ) -1 , (18)
where the variable Θ f k plays the same role for the parameters as the variable X f k defined in (9) for the states. We then have:

Θ f/a k = θ f/a,(1) k -θ f/a k , • • • , θ f/a,(Ne) k -θ f/a k √ N e -1 ∈ R N θ ×Ne (19)
with

θ f/a k = 1 N e Ne i=1 θ f/a,(i) k (20)
Knowing a better approximation of the model's parameters, we can update the state by EnKF (see Sec. 2.1.2). The Dual Ensemble Kalman filter allows to perform a recursive parametric inference / state estimation using the information from the ensemble members. The algorithm that we use is given in Appendix A.3.

From transport equation to multigrid resolution

The general expression for a conservation equation in local formulation over a continuous physical domain reads as:

Dx Dt = 1 ρ ∇ • σ + f (21)
where D/Dt is the total (or material) derivative of the physical quantity of investigation x and ρ is the flow density. The divergence operator is indicated as ∇• while σ is the stress tensor. Finally, f represents the effects of volume forces.

The evolution of the flow is obtained via time advancement of the discretized solution, which is performed in a physical domain where initial and boundary conditions are provided. A general expression of the discretized form of ( 21) for the time advancement from the step k -1 to k is given by:

x k = Φ k x k-1 + B k b k (22)
where Φ k is the state transition model which includes the discretized information of [START_REF] Xiao | Quantifying and reducing model-form uncertainties in reynolds-averaged navier-stokes simulations: A data-driven, physics informed bayesian approach[END_REF]. In case of non-linear dynamics described by [START_REF] Xiao | Quantifying and reducing model-form uncertainties in reynolds-averaged navier-stokes simulations: A data-driven, physics informed bayesian approach[END_REF], the state-of-the-art algorithms used for the discretization process are able to preserve the non-linear information in the product Φ k x k-1 , up to a discretization error which is usually proportional to the size of the time step. The term b k represents the control vector reflecting, for instance, the effect of the boundary conditions. B k is the control input model which is applied to the control vector b k . Equation ( 22) is consistent with a time explicit discretization of [START_REF] Xiao | Quantifying and reducing model-form uncertainties in reynolds-averaged navier-stokes simulations: A data-driven, physics informed bayesian approach[END_REF]. It is well known that this class of methods, despite the very high accuracy, may exhibit some unfavorable characteristics for the simulation of complex flows, such as limitations to the time step according to the Courant-Friedrichs-Lewy (CFL) condition [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. To bypass this limitation, one possible alternative consists in using implicit schemes for time discretization. In this case, the general structure of the discretized problem is usually cast in the form:

Ψ k x k = Ψ k x k-1 + B k b k = c k (23)
where Ψ k , Ψ k and B k are matrices obtained via the discretization process.

Obviously, considering Φ k = Ψ -1 k Ψ k , we retrieve [START_REF] Rochoux | Towards predictive data-driven simulations of wildfire spread -Part I: Reduced-cost Ensemble Kalman Filter based on a Polynomial Chaos surrogate model for parameter estimation[END_REF]. However, this manipulation is in practice not performed due to the prohibitive costs associated to large scale matrices inversions at each time step. Instead, an iterative procedure is used until the residual δ n , determined at the n-th iteration, falls below a pre-selected threshold value ε. In other words, the procedure is stopped when

δ n = Ψ k x n k -c k < ε.
Among the various iterative methods proposed in the literature, multigrid approaches are extensively used in CFD applications [START_REF] Hackbusch | Multi-grid methods and applications[END_REF][START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. The solution is found on the computational grid by updating an initial guess via multiple estimations obtained on a hierarchy of discretizations. Two well-known families of multigrid approaches exist, namely the algebraic multigrid method and the geometric multigrid method. With algebraic multigrid methods, a hierarchy of operators is directly constructed from the state transition model Ψ. On the other hand, the geometric multigrid obtains the solution via a set of operations performed in two (or more) meshes. In this paper, we consider the simplified case of two grids. Thereafter, the variables defined on the fine grid will be denoted with the superscript F (x F for instance), those defined on the coarse grid will be denoted with the superscript C (x C for instance).

The coarse-level representation x C is usually obtained suppressing multiple mesh elements from the initial fine-level one x F . This operation may be defined by a coarsening ratio parameter r C , which indicates the total number of elements on the fine grid over the number of elements conserved in the coarse grid. Among the numerous algorithms proposed for geometric multigrid, we use the Full Approximation Scheme (FAS), which is a well documented strategy [START_REF] Brandt | Multi-level adaptive solutions to boundary-value problems[END_REF][START_REF] Wesseling | Geometric multigrid with applications to computational fluid dynamics[END_REF]. A general formulation for a two-grid algorithm is now provided. The time subscript k is excluded for clarity. The superscript n represents the iteration step of the procedure.

1. Starting from an initial solution on the fine grid x 0 F (which is usually equal to x at the previous time step k -1), an iterative procedure is applied to obtain a first solution

x 1 F . A residual δ 1 F = c F -Ψ F x 1 F is calculated. 2.
x 1 F and δ 1 F are projected from the fine grid to the coarse grid space via a projection operator Π C , so that x 1 C and δ 1 C are obtained. Similarly, the state transition model Ψ F is projected on the coarse grid (that is reestimated based on the projection of the solution of the fine grid onto the coarse grid) to obtain Ψ C . Finally, we evaluate

c C = Ψ C x 1 C + δ 1 C .
3. An iterative procedure is employed to obtain x 2 C on the coarse grid using as initial solution x 1 C . 4. The updated variable on the fine grid is obtained as

x 2 F = x 1 F + Π F x 2 C -x 1 C
where Π F is a projection operator from the coarse grid to the fine grid. 5. At last, the final solution x 3 F is obtained via a second iterative procedure on the fine grid starting from the intermediate solution x 2 F . This procedure can be repeated multiple times imposing x 0 F = x 3 F at the beginning of each cycle. When the convergence is reached, the fine grid solution at time instant k is equal to x 3 F . Performing part of the calculations on a coarse grid level provides two main advantages [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. First, a significant reduction in the computational resources is obtained since the calculations performed over the coarse grid are usually much less expensive than a full set of iterations over the fine grid. Second, spurious high-frequency numerical oscillations, due to the discretization error in x F , are naturally filtered via the operators Π C and Π F , which globally improves the quality of the numerical solution. In this paper, the two projections (Π F and Π C ) are defined as 4-th order Lagrange interpolators.

Multigrid Ensemble KF method (MEnKF)

Despite the game-changing advantage that EnKF offers for the analysis of large-scale dynamical systems, the use of a sufficiently large ensemble (usually 60 to 100 members are required for convergence [START_REF] Asch | Data Assimilation: methods, algorithms, and applications[END_REF]) may still be too prohibitive for advanced applications. In the following, we present a Kalman filter strategy which relies on the use of a coarse mesh for the ensemble Kalman filter step. These computations on the coarse mesh are jointly run with a single high-refinement simulation, which is updated using the coarse mesh assimilation results. For this reason, the clock time required for the time advancement of the ensemble members and the memory storage of the physical variables are dramatically reduced. The multigrid-ensemble algorithm works through the steps described below. In the following description, the notation Ψ might hold for both Ψ and Ψ introduced in ( 23), depending on the choice of the time integration strategy. For its part, the state (x C k ) * represents the projection on the coarse-grid of (x F k ) f , the fine-grid solution. Finally, (x C k ) denotes the coarse-grid state obtained by an ensemble Kalman filtering process applied on the coarse grid. In the FAS multigrid method described in Sec. 2. 2,(x C k ) * and (x C k ) correspond to x1 C and x 2 C , respectively.

1. First iteration on the fine grid. Starting from an initial solution on the fine grid x F k-1 a , a forecasted state (x F k ) f is obtained by using θ a k as parameter for the model Ψ F .

Projection on the coarse grid. (x F

k ) f is projected on the coarse grid space via a projection operator Π C , so that (x C k ) * is obtained. Similarly, the state matrix Ψ F is projected on the coarse grid to obtain Ψ C . This critical passage will be discussed in detail in the following.

Time advancement of the ensemble members used in the Dual

EnKF. For each member i of the ensemble, the state matrix (Ψ C ) (i) used for the advancement in time on the coarse grid is determined f,e is used to determine the coarse grid solution (x C k ) through a Kalman filter estimation. 5. Final iteration on the fine grid. A first estimation (x F k ) of the fine grid state is determined using the results obtained on the coarse space:

(x F k ) = (x F k ) f + Π F (x C k ) -(x C k ) * . The state (x F k )
a is obtained from a final iterative procedure starting from (x F k ) . An overview of the assimilation cycle described in the previous algorithm is presented in Fig. 1. Two important aspects need to be discussed:

FINE GRID COARSE GRID

-As previously stated, Ψ C and its role in the determination of the matrices

(Ψ C ) (i)
is an essential step in the MEnKF strategy. In non-linear problems of interest in fluid mechanics, the state transition matrix Ψ includes information of the multi-scale interactions that are specific for every case investigated. The simplest possible choice, which is the one adopted in this work, is to calculate the coefficients of the matrices Ψ C and (Ψ C )

(i)
separately for each simulated state. Thus, the similarities between the employed state matrices are limited to the use of the same discretization schemes / structure of Ψ. However, one can envision to use the non-linear information conserved in Ψ C , which is supposedly accurate, to improve the accuracy of the prediction of the ensemble members. This aspect is discussed in the perspectives included in Sec. 7.

-The recursive structure of the algorithm allows for integration of iterative corrections for non-linear systems [START_REF] Sakov | An iterative enkf for strongly nonlinear systems[END_REF] as well as hard constraints (see the discussion in the introduction of [START_REF] Simon | Kalman filtering with state equality constraints[END_REF][START_REF] Nachi | Kalman filtering in the presence of state space equality constraints[END_REF][START_REF] Zhang | Regularized ensemble kalman methods for inverse problems[END_REF]) to respect the conservativity of the model equations. However, these corrections may result in an increase of the computational resources required. Here, the multigrid algorithm itself is used for regularization (i.e. for smoothing the discontinuities in the physical variables produced by the update via Kalman Filter) of the flow. In fact, if an intentionally reduced tolerance is imposed in the iterative steps 4 and 5, the final solution will keep memory of the features of the state estimation produced in step 3. However, the iterative resolution will smooth the estimation via the state transition model Ψ, which will perform a natural regularization of the flow. Clearly, if a reduced tolerance is imposed, the final solution will not necessarily respect the conservativity constraints of the model equations. However, one can argue that complete conservativity is not an optimal objective in this case, if the model state at the beginning of the time step is not accurate.

The advantages of our strategy with respect to classical approaches based on EnKF may be summarized in the following points:

-The RAM requirement necessary to store the N e ensemble members during the assimilation is usually moderate. The reduction in computational costs is driven by N e and by the size of the coarse variables. To illustrate, let us consider the case of a simple two-level geometric multigrid approach for a 3D test case with a constant coarsening ratio r C = 4 and a size of ensembles N e = 100. Each ensemble member is then described by -Considering that the ensemble members in the coarse grid and the simulation over the fine grid are running simultaneously, communication times are optimized.

-Owing to the iterative procedures of steps 4 and 5, regularization of the final solution is naturally obtained.

-The algorithm is here described and tested in the framework of geometric multigrid, but it can actually be integrated within other algorithmic structures. For iterative methods, the only essential operation to be performed is the determination of the state transition model Ψ C and of the projections Π C and Π F . This implies that the method can be easily extended to other popular procedures, such as the algebraic multigrid. If the multigrid operations are removed, the methods becomes a classical multilevel EnKF method, which relies on just one level of resolution for the ensemble members. However, in this case, no regularization is obtained unless specific corrections are included.

This general algorithm may be easily tailored accounting for the complexity of the test case investigated, in particular for the requirements of iterative loops on both the coarse grid level and the fine grid level. The algorithm that we used to validate our approach is described in Appendix A.4.

Application: one-dimensional Burgers' equation

The MEnKF method introduced in Sec. 3 is now applied to the analysis of different test cases. Several dynamical systems of increasing complexity were chosen in order to highlight different properties of the algorithm. Also, a set of different tests is performed in order to obtain a comprehensive validation of the method. At first, let us consider a 1D Burgers' equation:

∂u ∂t + u ∂u ∂x = 1 Re ∂ 2 u ∂x 2 ( 24
)
where x is the spatial coordinate, u the velocity and Re is the Reynolds number.

Equation ( 24) is non-dimensionalized with a reference velocity u 0 and a reference length λ. This equation is solved with a second-order centered finite difference scheme for the space derivatives and a first-order scheme for the time integration to obtain the general form of discretized representation as given by [START_REF] Simon | Kalman filtering with state equality constraints[END_REF]. A Dirichlet time-varying condition is imposed at the inlet:

u(x = 0, t) = u 0 (1 + θ 1 sin(ωt + θ 2 )) (25)
where u 0 = 1 and ω = 2πf c . The characteristic frequency f c is set to 1 and the characteristic length λ is defined as λ = u0 fc . The characteristic time is defined as t c = λ/u 0 . This implies that, over a complete period of oscillation 1/f c , the solution is advected of a length λ by the velocity u 0 . θ 1 and θ 2 represent the amplitude and phase of the sinusoidal signal, respectively. Note that from now on, the units are the characteristic magnitudes. The outlet boundary condition is extrapolated from the nearest points to the outlet. The initial condition imposed for t = 0 is u(x, t = 0) = u 0 everywhere in the physical domain (θ 2 = 0 for the reference simulation). The value of the Reynolds number is Re = 200.

The time advancement step is chosen as ∆t = 0.0002. It is kept constant throughout the simulation. The analysis is performed over a physical domain of size [0,[START_REF] Mons | Kriging-enhanced ensemble variational data assimilation for scalar-source identification in turbulent environments[END_REF]. The distance between the computational nodes in the fine mesh is constant and set to ∆x = 0.0125. This choice has been performed to discretize the characteristic length λ using 80 mesh elements. This also implies that the total number of nodes employed to perform the calculation is N x = 800. A simulation of reference is run on the fine grid with values θ 1 = 0.2 for the amplitude and θ 2 = 0 for the phase. Thereafter, the solution obtained by this reference simulation is called the true state or truth. A flow visualization at t = 10 is shown in Fig. 2. For the investigated value of Reynolds number, the non-linear effects and viscous mechanisms can be clearly identified. Data assimilation is performed in the following conditions:

-The observations are sampled each 30 time steps of the reference simulation on the space domain [0, 1] (80 sensors) and on the time window [START_REF] Mons | Kriging-enhanced ensemble variational data assimilation for scalar-source identification in turbulent environments[END_REF][START_REF] Siripatana | Combining ensemble kalman filter and multiresolution analysis for efficient assimilation into adaptive mesh models[END_REF]. The observations are sampled as soon as the flow is fully developed. These observations are artificially perturbed using a Gaussian noise of variance R = 0.0025.

-The model is chosen to be the discretized version of [START_REF] Nachi | Kalman filtering in the presence of state space equality constraints[END_REF]. The numerical test consists of one main simulation, which is run on the fine grid previously introduced, and an ensemble of N e = 100 coarse simulations used for assimilation purposes. The coefficients θ 1 and θ 2 are initially assumed to be described by Gaussian distributions, so that θ 1 ∼ N (0, Q θ1 ) and The data estimation is run for a time window of T DA = 19 characteristic times, which encompass roughly 3000 DA analysis phases. The sensitivity of the parametric inference procedure to the resolution of the coarse simulations is investigated considering several coarsening ratios r C = 1, 2, 4, 8, 16. The fine grid is unchanged, so that the MEnKF is performed using information of progressively coarser grids as r C increases. The interest of this test is to analyse the loss of accuracy of the estimator as r C increases and to ascertain the potential for efficient trade-off between accuracy and computational resources required for the estimation process.

θ 2 ∼ N (0.3, Q θ2 ).
The time-evolution of the estimation of θ 1 is shown in Fig. 3. Very rapid convergence (less than 2t c ) is observed for r C ≤ 4. In addition, the parameter estimation is extremely precise (discrepancy lower than 0.01% for r C = 1, lower than 2% for r C = 4). For higher values of the parameter r C , the estimation of θ 1 becomes progressively more degraded. For the case r C = 8, θ 1 is initially overestimated and it finally converges to a value of θ 1 = 0.194, 3% smaller that the true value. Larger errors in the optimization of θ 1 are observed for r C = 16.

In this case, the optimized amplitude parameter is θ 1 = 0.27, which is 35% larger than the true value. the parameter θ 2 , which is shown in Fig. 4. For r C = 1, 2, 4, we obtain an accurate prediction of the parameter, while a loss in accuracy is observed for the cases r C = 8, 16. This observation can be justified considering the number of mesh elements representing one characteristic length λ for these cases, which are 10 and 5 for r C = 8 and 16, respectively. If one considers that the frequency ω is set so that one complete oscillatory cycle is performed on average over a characteristic length λ, this means that the average phase angle between mesh elements is equal to 0.63 radians for r C = 8 and 1.26 radians for r C = 16.

Thus, the values observed for the optimization of θ 2 for these two cases, which are around θ 2 ≈ 0.15, are significantly lower than the uncertainty due to the coarse-level resolution. State estimation results for the case r C = 1 are shown in Fig. 5. This case, which is performed using the same grid for the coarse and fine mesh level, is equivalent to a standard Dual EnKF. However, owing to the final multigrid iterative loop, the final solution is naturally regularized. The results, which are shown for t = 1, 3.88, 10.60, show that the estimator successfully represents the behaviour of the dynamical system. A full domain advective time (i.e. 10 characteristic time units) must be simulated in order to observe the effect of the MEnKF in the whole domain. In fact, the parametric information imposed at the inlet for the ensemble members must affect the whole physical domain before a reliable correlation between the state variables can be established. However, once this initial transient is faded, the state estimation almost perfectly captures the behaviour of the true state.

Results are now investigated for increasing values of r C . Results for r C = 8 are shown in Fig. 6. Minor differences between the state estimation and the true state can be observed in this case. This discrepancy is due to the lack of resolution of the ensemble members. In fact, the resolution in this case is of 10 mesh elements per characteristic length. This number of points is arguably not enough to provide an accurate representation of the sinusoidal waves which are imposed at the inlet. However, one can see that no spurious numerical effects are observed as the estimator provides a smooth, continuous prediction of the velocity. The discrepancy between the true state and the state estimation is mainly associated with an erroneous calculation of the Kalman gain due to the under-resolution of the ensemble, which also affects the parameter estimation.

A combined analysis of Fig. 3 and6 shows that, due to the lack of accuracy in the estimation of θ 1 , the variable u is over-predicted for t < 2t c while it is slightly under-predicted for t > 4t c . At last, the case for r C = 16 is shown in Fig. 7. In this case, the mesh elements are only 5 times smaller than the characteristic length λ. Despite the important under-resolution of the ensemble members, which severely affects the estimation of θ 1 and θ 2 , the state estimation still adequately represents the main features of the dynamical system. The discrepancy between the truth and the state estimation is measured via the time-dependent relative Root Mean Square Error (RMSE), i.e.

RMSE(k) = x (u F k ) a (x) -(u F k ) True (x) 2 dx x (u F k ) True (x) 2 dx ( 26
)
The results are shown in Fig. 8 for different values of the coarsening ratio r C .

One can see that the error decreases slowly for t < 10t c . This threshold time corresponds to a complete advective cycle in the physical domain. After this transient, the error may rapidly decrease before reaching a quasi-asymptotic behaviour. One can also see that, once convergence is reached, the asymptotic RMSE value decreases with lower r C values, as expected. However, lower r C values are as well associated with larger computational costs, so that a trade-off between accuracy and required resources must be found. This aspect is further investigated considering the computational resources required to perform a full assimilation window for each r C value. In Fig. 9, results are shown and normalized over the case r C = 1. One can see that the computational resources required rapidly decrease with increasing r C , even for this simple one-dimensional test case. For large values of r C = 8, 16, one can see that the computational resources reach a plateau. Here the computational time to perform the DA procedures, which is the same for every r C , is of similar order of magnitude of the calculation for the time advancement of the ensemble members. In summary, the present analysis assesses the performance of the MEnKF tool for varying mesh resolution of the ensemble members. As expected, the accuracy of the state and parameter estimations diminishes for increasing r C , but so does the computational cost. In addition, it was observed that the accuracy significantly drops when the mesh resolution is not able to provide a suitable description of the main scales characterizing the flow. Such a significant discrepancy for a relatively simple test case stresses how a minimal resolution threshold must be achieved in order to capture the essential physical features and to obtain a successful state estimation. In conclusion, the results of the parametric optimization may become unreliable with extreme under-resolution of the ensemble members. The dispersive and diffusive effects of the coarse grid on the state prediction deserve an in-depth analysis, which is out of the scope of this work. Here, the discretization schemes used are the same for the fine and coarse grids. However, one could argue that choosing adapted schemes for progressively coarser grids could improve the performance of the estimator.

Acoustic propagation of sinusoidal wave

The MEnKF strategy is now applied to a more complex physical system, namely the inviscid one-dimensional Euler equations:

∂ρ ∂t + ∂(ρu) ∂x = 0 (27) ∂(ρu) ∂t + ∂((ρu)u) ∂x + ∂p ∂x = 0 (28) ∂ (ρE) ∂t + ∂((ρE) u) ∂x + ∂(pu) ∂x = 0 ( 29
)
where ρ is the density, u is the velocity, p is the pressure and E is the total energy per unit mass. In this case, viscous effects are absent, but acoustic propagation affects the evolution of the flow. The equations are discretized using the finite difference method. Second-order centered schemes are used for the derivatives in space and a first-order scheme for the time integration to obtain the general form of discretized representation as given by [START_REF] Simon | Kalman filtering with state equality constraints[END_REF]. A centered sixth-order numerical filter is included to damp numerical spurious oscillations [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for flow and noise computations[END_REF]. The numerical algorithm is used to analyse the acoustic propagation of a sinusoidal wave with a time-varying amplitude. To do so, a Dirichlet time-varying condition is imposed at the inlet for the velocity:

u(x = 0, t) = u 0 (1 + θ(t) sin(ωt)) ( 30
)
u 0 is set in order to impose an inlet Mach number M = u0 a = 0.4, where a is the speed of the sound. The amplitude of variation in θ is set at sufficiently low value to allow a flow evolution mainly driven by acoustic phenomena. The inlet velocity perturbation creates an acoustic wave that is transported along the domain with a speed of u 0 + a. The characteristic velocity and length of the test case are u c = u 0 + a and λ. Similarly to the analysis in Sec. 4, λ represents the wavelength of the signal imposed at the inlet. The characteristic time of the system is defined as t c = λ/u c .

The sinusoidal behaviour of the velocity at the inlet is characterized by a constant frequency f c = 1/t c with ω = 2πf c . However, the time-varying amplitude of the sinusoidal wave is driven by the parameter θ(t) = θ 0 (1 + sin(ω θ t)), where θ 0 , ω θ = 2πf c /b and b are constants. The density at the inlet is fixed ρ(x = 0, t) = ρ 0 as well as the total energy per unit mass E(x = 0, t) = E 0 , which is calculated as E 0 = e + 0.5u 2 0 . The outlet boundary condition is extrapolated from the nearest points to the outlet. The initial condition imposed at t = is u(x, t = 0) = u 0 , ρ(x, t = 0) = ρ 0 and E(x, t = 0) = E 0 everywhere in the physical domain. The fluid is an ideal gas with γ = 1.4, ρ 0 = 1.17 and T 0 = 300 in S.I. units. Note that from now on, the units of measure are the characteristic magnitudes.

The computational domain has been set to a size of L x = 10. A uniform mesh distribution is used for every calculation. Similarly to the analysis in Sec. 4, 80 mesh elements are used to discretize a characteristic length λ for a total of N x = 800 elements in the domain. Finally, the normalized value of ∆t is set to ∆t = 0.0006. The DA procedure is identical to the analysis presented in Sec. 4. The assimilation is composed by the base simulation, which is run on the fine grid, and N e = 100 simulations on a coarse-grid level. In this section, only one coarsening ratio is considered (r C = 4). The estimator is used to dynamically track the value of the parameter θ, which evolves in time. No a priori knowledge about the behaviour of the parameter is used. This time optimization is much more challenging when compared with the inference of a constant parameter as done in Sec. 4. A similar analysis using a classical Kalman smoother was recently proposed by Mons et al. [START_REF] Mons | Reconstruction of unsteady viscous flows using data assimilation schemes[END_REF].

For each coarse grid simulation of the estimator, θ is initially assumed to be a random Gaussian phenomenon θ ∼ N (0, Q θ ). The initial value of the covariance is Q θ (t = 0) = 6.4 × 10 -5 . As shown in the previous case, the value initially imposed at the inlet of the fine-grid simulation is θ = 0, while random values are selected for each ensemble member on the coarse grid level following the normal distribution introduced above. The variance of the parameter θ for the ensemble members is artificially increased, as in the classical Dual EnKF algorithm. As described in Appendix A.4, we add to the estimated parameter of each member of the ensemble a Gaussian noise of diagonal covariance Σ θ k = 10 -10 . Extensive numerical tests have been performed and the results show the importance of Σ θ k . The value chosen for Σ θ k is arbitrary and works only for this test-case. In general, the estimation of Σ θ k is challenging, since a priori information about the dynamic behaviour of θ is required. Nonetheless, θ can be determined heuristically: the artificial variance added to the parameters, Σ θ k , should be of the same order of magnitude as the true rate of change of θ in a single assimilation cycle. When this parameter is underpredicted, the variance in θ is too small for the Dual EnKF to perform, while showing a oscillatory behaviour in the parameter estimation if Σ θ k = 10 -10 is overpredicted.

This study is performed for three different values of f a , the normalized period between successive assimilations defined as f a = t c /t a . This normalized analysis frequency corresponds to the number of analysis phases per characteristic time of simulation t c . Three different values of f a are investigated: f a = 2, 10, 55.

The estimator is run for a total simulation time T ref = 110, which encompasses 220 to 6000 DA analysis phases, depending on the value of f a . At the end of each analysis, the mean value and the variance of the amplitude θ are updated following the Dual EnKF technique [START_REF] Hamid | Dual state-parameter estimation of hydrological models using ensemble kalman filter[END_REF], similarly to what was done in Sec. 4.

The results for the estimation of the time-varying parameter θ are reported in Fig. 11. The time evolution of θ are correctly estimated for the three values of f a . This is an important result, considering that no a priori information was provided for the evolution of this parameter. A more detailed analysis reveals a lag in the parameter estimation. The application of a simple Kalman filter seems to be responsible for this result, while a Kalman Smoother (KS) should have been used to obtain a better synchronization. However, considering that the implementation of a KS is straightforward in this case and that observation is always provided close to the inlet, we considered that the increase in computational resources required by the KS were not deserved. We find that the lag increases when a relatively small number of DA analyses is done. One can see that the prediction is significantly degraded for f a = 2, while similar results are obtained for f a = 10, 55. In addition, θ tends to be underestimated (around 10%) when it reaches its maximum value. This result is arguably associated with the under-resolution of the coarse level of the grid, where the gradients of physical variables are calculated with lower accuracy. Now, results dealing with the state estimation are discussed. The predicted physical variable ρu, normalized over the initial value ρ 0 u 0 , is shown in Fig. 12, 13 and 14 for f a = 2, 10 and 55, respectively. For f a = 2, the state estimation is significantly distant from the truth. It appears that the field correction applied via the Kalman gain is not able to compensate the poor estimation of θ.

However, accurate results are observed for f a = 10 and 55. Even though the value of the parameter θ is not exact, the total state estimation including the correction via Kalman gain is very precise. For the case f a = 55, almost no discernible difference is observed between the state estimation and the truth.

At last, the relative Root Mean Square Error (RMSE) defined as

RMSE(k) = x (ρu) F k a (x) -(ρu) F k True (x) 2 dx x (ρu) F k True (x) 2 dx (31)
is shown in Fig. 15. The error achieves a quasi-constant asymptotic behaviour after a complete propagation of the signal in the physical domain (t ≈ 10t c ).

As expected, a low global error is obtained for the cases f a = 10 and f a = 55.

On the other hand, the error for f a = 2 case is around 2 -3 times larger. The very small difference in performance between the cases f a = 10 and f a = 55 indicates that, once a minimal threshold in the assimilation period t a is reached, the prediction of the estimator exhibits a robust convergence. Similar results were previously observed in three-dimensional simulations [START_REF] Meldi | Augmented Prediction of Turbulent Flows via Sequential Estimators: Sensitivity of State Estimation to Density of Time Sampling for Available Observation[END_REF].

Spatially evolving compressible mixing layer

In this section, we consider the compressible Navier-Stokes equations in a two-dimensional physical domain:

∂ρ ∂t + div(ρu) = 0 (32)
∂ (ρu) ∂t + div(ρu ⊗ u) = -grad p + div τ (33)
∂ (ρE) ∂t + div(ρEu) = -div(pu) + div(τ u) + div (λ(T )grad T ) ( 34
)
where ρ is the density, u is the velocity (components u in the streamwise direction and v in the normal direction), p is the pressure, E is the total energy per unit of mass, τ is the tensor of the viscous constraints and T is the temperature. To obtain the representation given by ( 23), the equations are discretized using the finite difference method. Second-order centered schemes are used for the derivatives in space and a first-order scheme for the time integration. A centered sixth-order numerical filter is included to damp numerical spurious oscillations [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for flow and noise computations[END_REF].

As flow configuration, we consider the two-dimensional spatially evolving mixing layer at Re = 100. For this value of Reynolds number, the flow exhibits unsteady features. It can be shown [START_REF] Ko | Sensitivity of two-dimensional spatially developing mixing layers with respect to uncertain inflow conditions[END_REF][START_REF] Meldi | Quantification of errors in large-eddy simulations of a spatially evolving mixing layer using polynomial chaos[END_REF] that the characteristics of the mixing layer are strongly affected by the inlet and, in particular, by imposed ad-hoc time perturbations. The computational domain has been set to a size of 14Λ × 6Λ in the streamwise direction x and normal direction y, respectively.

The characteristic length Λ, which is taken as reference length from now on, is given by Λ = Aδ 0 , where δ 0 is the initial vorticity thickness imposed at the inlet. The value of the parameter A is determined from the most unstable wavelength determined by Linear Stability Theory (LST). At Re = 100, we have A = 14.132. The mesh resolution in the horizontal direction is constant for x ≤ 10. The size of the elements is ∆x = δ0 8 . For x ≥ 10, a sponge zone is established with a coarsening ratio between successive elements which increases from 1.025 to 1.04. The resolution in the normal direction is constant and equal to ∆y = δ0 20 for -0.18 ≤ y ≤ 0.18. Outside this zone, the mesh elements increase in size moving away from the centerline with a constant coarsening ratio of 1.01. using a constant single mode for the inlet perturbation 1 = = cte. In the second case, 1 varies in time following a sinusoidal form: 1 = (1 + sin(ω t)). The values of the numerical parameters characterizing the perturbation are = 0.15, n 1 = 0.4π, ω 1 = 1/t c and ω = 0.62ω 1 , where t c = 2Λ/(U 1 + U 2 ) is the average advection time. A flow visualization of ρv at t = 10 is shown in Fig. 16 for the two cases. In both cases, one can clearly observe the emergence of coherent structures. When a constant value of 1 is used, the coherent structures are virtually aligned. When 1 is time-varying, one can observe the emergence of more complex pairing patterns. The DA procedure is performed using the following elements:

-The model is the discretized version of the system given by ( 32) - [START_REF] Brunton | Closed-Loop Turbulence Control: Progress and Challenges[END_REF].

The features of the fine mesh level were previously introduced. For the coarse grid level, a homogeneous coarsening ratio r C = 4 is employed. The simulation is initialized in the physical domain using a simple hyperbolic tangent profile with no perturbation (see [START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using monte-carlo methods to forecast error statistics[END_REF]). We consider that no prior information is available on the time evolution of 1 . At t = 0, this coefficient is fixed to be a random Gaussian phenomenon 1 ∼ N (0, Q a ) where the initial value of Q a (t = 0) = 0.0625. Similarly to the cases analyzed in Sec. 4, the value imposed on the main fine-grid simulation at t = 0 is 1 = 0, while random values are imposed for each ensemble member on the coarse grid level. The size of the ensemble is N e = 100.

-The observation is sampled from the reference simulations shown in Fig. 16 The value of ρv that is predicted at the centerline of the mixing layer is shown in Fig. 18. At t = 1, the flow is strongly affected by the KF corrections. After five characteristic times, the flow is perfectly matching the true state upstream. Finally, after the initial transient is dissipated, we observe that the state estimation almost perfectly matches the true state in the whole physical domain.

The results dealing with the time-varying parameter 1 are now discussed.

The time evolution of the estimated value of 1 is reported in Fig. 19. The overall sinusoidal trend is generally respected, although a relatively small phase lag is visible. This lag does not appear to be larger than the one previously observed for the one-dimensional case based on the Euler equation. The presence of this delay has probably the same reasons as in Sec. 5. However, in this case, some over prediction of the parameter is locally observed in time, which was not obtained for the wave propagation test case.

The results obtained for the prediction of the normal momentum ρv are shown in Fig. 20. One can see that the combination of parameter and state estimations produces an accurate prediction of the flow. Minor differences are observed with the true state. In particular, the momentum ρv does not exhibit spurious oscillations which could stem from the field correction determined via the Kalman gain. In order to evaluate the respective influence of the parameter estimation step and state estimation phase, a test case is run in which only the parameter estimation is performed. That is, the state estimation obtained on the coarse-grid level is not included in the steps 4 and 5 of the algorithm presented in Sec. 3. While the results of the parameter estimation are the same for the two cases, one can see in Fig. 21 that the prediction is sensibly deteriorated.

This observation is quantified by the evaluation of the relative Root Mean Square Error (RMSE), defined as:

RMSE(k) = x (ρv) F k a (x) -(ρv) F k True (x) 2 dx x (ρv) F k True (x) 2 dx (38)
The results, which are shown in Fig. 22, indicate that the accuracy of the complete algorithm is higher when compared to the case in which only the parameter estimation is performed. Therefore, the two operations concurrently provide an improvement in the prediction of the flow.

At last, an analysis of the conservativity of the algorithm is performed. As previously discussed, the state estimation obtained via EnKF does not necessarily comply with the dynamical equations of the model. This drawback can be responsible for discontinuities in the physical field, which can significantly affect the accuracy and stability of the global algorithm. The analysis is performed considering an indicator Γ F k which measures the conservation of the transversal momentum equation [START_REF] Stonebridge | Diagonal Approximations to the ObservationError Covariance Matrix in Sea Ice ThicknessData Assimilation[END_REF] in discretized form:

(ρv) F k a -(ρv) F k-1 a ∆t -F ρv ρ k , u k , p k , τ k = Γ F k , (39)
where F ρv represents the spatial discretization terms in the transversal momentum equation. In the forecast step performed via the model, Γ F k = 0 down to a convergence rate δ which is prescribed. However, the value of Γ F k , at the end of a time step where a forecast-analysis is performed, is strictly connected with the computational strategy employed. Here, three cases are considered for -A classical Dual EnKF is performed on the coarse grid and a fine-grid correction is obtained through the ensemble statistics. In this scenario, the state estimation obtained in the step 4 of the MEnKF algorithm presented in Sec. 3 is directly projected in the fine mesh space and used as final solution. The step 5 of the algorithm is not performed.

-A standard MEnKF algorithm, as described in Sec. 3.

-A MEnKF algorithm where the ensemble prediction is just used to estimate the unknown parameter of the system. No update of the physical solution is performed using the correction via Kalman gain.

The results are shown in Fig. 23 after the first forecast / analysis step. For clarity, we introduce a normalized criterion (Γ

F k ) * = Γ F k C where C is defined as max k ((ρv) F k ) a -((ρv) F k-1 ) a ∆t
. As expected, Γ F k = 0 everywhere when MEnKF is only used for the parameter estimation. Here, the time advancement of the solution is performed using the model only, which exactly complies with the discretized equation and respects conservativity (up to a convergence error which is negligible). On the other hand, results in Fig. 23 (a) show some lack of conservativity in the physical domain for the first scenario. This is also expected, since no constraint is imposed to force the Kalman gain correction to comply with the dynamical equations. Finally, results for the MEnKF are shown in Fig. 23 (b). The evolution of (Γ F k ) * is very similar to the results observed for the first scenario. However, one can clearly see that this field appears to be sensibly smoothed out by the multigrid iterative procedures in step 4 and 5 of the MEnKF algorithm. As previously discussed, complete conservativity starting from an erroneous state at k -1 is possibly not an optimal objective, while one wants a regularized solution to avoid affecting the precision of the global calculation. On this last objective, the MEnKF appears to provide a better result when compared with the classical Dual EnKF, described in the first scenario.

Considering also that the MEnKF showed better accuracy than the algorithm relying on parameter estimation only, one can conclude that the MEnKF provides an efficient compromise between global accuracy and regularization of the solution. In order to draw more information about this important aspect, the MEnKF algorithm needs to be tested for the simulation of three-dimensional compressible flows, where the Kalman gain correction may be responsible for important acoustic phenomena which are not observable in 2D and 1D dynamical systems.

Conclusions

A sequential estimator based on a Kalman filter approach for Data Assimilation of fluid flows is presented in this research work. This estimator exploits iterative features which are employed in several CFD codes for the resolution of complex applications in fluid mechanics. More precisely, the multilevel resolution associated with the multigrid iterative approach for time advancement is used to generate several low-resolution numerical simulations. These results are then employed as ensemble members to determine i) the correction via Kalman filter, which is then projected on the high-resolution grid to correct a single simulation which corresponds to the numerical model and ii) an optimization of the free parameters driving the simulation. One of the main advantages of the model is that, owing to the iterative procedure for the calculation of the flow level to improve the predictive capabilities of the ensemble members. This implies a more accurate prediction of the state estimation and of the parameter optimization on the coarse level, from which the main refined simulation will benefit. This loop has the potential to improve even more the performance of the MEnKF model, and strategies for efficient application are currently under investigation. Moreover, we plan to test the MEnKF algorithm to more complex configurations involving complex geometries and more challenging parametric Compute the Kalman gain

K k = P f k H k H k P f k H k + R k -1 x a k = x f k + K k y o k -H k x f k P a k = (I -K k H k )P f k Analysis/Correction
K e k = X f k Y f k Y f k Y f k + E o k (E o k ) - 1 7
k-1 + τ (i) k with τ (i) k ∼ N (0, Σ θ k ) x f,(i) k = M k:k-1 (x a,(i) k-1 , θ f,(i) k ) y f,(i) k = H k x f,(i) k 4:
Compute the normalized anomalies Parameter update: Compute the normalized anomalies

K θ,e k = Θ f k Y f k Y f k Y f k + E o k (E o k ) -1 θ a,(i) k = θ f,(i) k + K θ,e
X f k :,i = x f,(i) k -x f k √ N e -1 ; Y f k :,i = y f,(i) k -y f k √ N e -1 ; [E o k ] :,i = o,(i) k -o k √ N e - 1 8:
State update:

K x,e k = X f k Y f k Y f k Y f k + E o k (E o k ) -1
x a,(i) k Parameter forecast:

θ f,(i) k = θ a,(i) k-1 + τ (i) k with τ (i) k ∼ N (0, Σ θ k ) (x C k ) f,(i) = M C k:k-1 x C k-1
a,(i) , θ f,(i) k

if Observation available then

(y C k ) f,(i) = H C k (x C k )
f,(i)

3:

Observation ensemble:

(y C k ) o,(i) = (y C k ) o + ( C k ) o,(i) with ( C k ) o,(i) ∼ N (0, R C k ) (R C k ) e = 1 N e -1 Ne i=1 ( C k ) o,(i) ( C k ) o, (i) 4:
Compute the normalized anomalies Parameter update:

(K C k ) θ,e = Θ f k Y f k Y f k Y f k + E o k (E o k ) -1 θ a,(i) k = θ f,(i) k + (K C k ) θ,e (y C k ) o,(i) -(y C k ) f, (i) 6:
State forecast:

(x C k ) f,(i) = M C k:k-1 x C k-1 a,(i) , θ a,(i) k (y C k ) f,(i) = H C k (x C k )
f,(i)

7:

Compute the normalized anomalies

X f k :,i = (x C k ) f,(i) -(x C k ) f √ N e -1 ; Y f k :,i = (y C k ) f,(i) -(y C k ) f √ N e -1 ; [E o k ] :,i = ( C k ) o,(i) -( C k ) o √ N e - 1 8:
State update:

(K C k ) x,e = X f k Y f k Y f k Y f k + E o k (E o k ) -1 (x C k ) a,(i) = (x C k ) f,(i) + (K C k )
x,e (y C k ) o,(i) -(y C k )

f,(i)

Algorithm 4: Multigrid EnKF algorithm. We have i = 1, • • • , N e . Initialize (x F 0 ) a , θ a 0 , θ a,(i) 0

, (x C 0 )

a,(i)

for k = 1, . . . , K do 2:

Fine grid forecast:

(x F k ) f = M F k:k-1 x F k-1
a , θ a Fine grid state correction using the ensemble statistics:

(x C k ) = (x C k ) * + (K C k ) x,e (y C k ) o -H C k (x C k ) * (x F k ) a = (x F k ) f + Π F (x C k ) -(x C k ) * 6:
Matrix-Splitting iterative procedure on the final solution starting from (x F k ) a .

Figure 1 :

 1 Figure 1: Schematic representation of the Multigrid Ensemble Kalman Filter (MEnKF). Two different levels of representation (fine and coarse grids) are used to obtain a data-driven fine grid estimation. The Dual Ensemble Kalman filter procedure is solved in the coarse grid. The full algorithm is given in Appendix A.

Figure 2 :

 2 Figure 2: Solution of the 1D Burgers' equation at t = 10 for θ 1 = 0.2 and θ 2 = 0 (true state).

  The initial value of the covariance of the parameters is chosen equal to Q θ1 (t = 0) = Q θ2 (t = 0) = 0.0025. The values prescribed on the fine grid simulation are the mean values of the Gaussian distribution i.e. θ 1 = 0 and θ 2 = 0.3. Random values for the parameters are imposed at the inlet for each ensemble member on the coarse grid level.The initial mean values for the parameters are significantly different when compared with the values prescribed in the reference simulation, which are θ 1 = 0.2 and θ 2 = 0. This choice allows to analyze the rate of convergence of the optimization procedure. Also, the initial condition u(x, 0) = u 0 is imposed for every simulation.

  Similar considerations can be drawn by the analysis of the optimization of

Figure 3 :

 3 Figure 3: Values of the parameter θ 1 for different coarsening ratios rC = 1, 2, 4, 8, 16. In the zoomed region, the shaded area represents the 95% credible interval for the shown cases.

Figure 4 :

 4 Figure 4: Values of the parameter θ 2 for different coarsening ratios rC = 1, 2, 4, 8, 16. In the zoomed region, the shaded area represents the 95% credible interval for the shown cases.

  (a) t = 1 (b) t = 3.88 (c) t = 10.60

Figure 5 :

 5 Figure 5: Estimations obtained by MEnKF for rC = 1 at t = 1 (a), t = 3.88 (b) and t = 10.60 (c). Times are given in tc units. The grey shaded area corresponds to the observation window.

  (a) t = 1 (b) t = 3.88 (c) t = 10.60

Figure 6 :

 6 Figure 6: Estimations obtained by MEnKF for rC = 8 at t = 1 (a), t = 3.88 (b) and t = 10.60 (c). Times are given in tc units. The grey shaded area corresponds to the observation window.

  (a) t = 1 (b) t = 3.88 (c) t = 10.60

Figure 7 :

 7 Figure 7: Estimations obtained by MEnKF for rC = 16 at t = 1 (a), t = 3.88 (b) and t = 10.60 (c). Times are given in tc units. The grey shaded area corresponds to the observation window.

Figure 8 :

 8 Figure 8: Time evolution of the RMS error of u for rC = 1, 2, 4, 8, 16.

Figure 9 :

 9 Figure 9: Computational time required to perform a full assimilation cycle for rC = 1, 2, 4, 8, 16. Results are normalized over the computational time required for rC = 1.

A

  preliminary simulation is performed using an explicit Euler scheme for time derivative. The parameters which describe the amplitude of the velocity oscillation imposed at the inlet are set to θ 0 = 0.015 (true value) and b = 10. A flow visualization of the wave patterns is shown in Fig. 10 at t = 17.3. The fully developed state obtained for t = 10 is used to initialize a new simulation from t = 0. This simulation is run for a total time of T ref = 110. As in Sec. 4, this simulation is used as true state and it is sampled every 30 time steps to obtain observations in the space region x ∈ [0, 1]. The sampled data, which is used as observation, is artificially perturbed using a Gaussian noise of variance R = 0.09.

Figure 10 :

 10 Figure 10: Field ρu of the one-dimensional inviscid Euler equations at t = 17.3 for θ 0 = 0.015 and ω θ = ω/10 (true state).

Figure 11 :

 11 Figure 11: Time estimation of the parameter θ driving the amplitude of the sinusoidal acoustic wave. Results are shown for fa = 2, 10, 55 and compared to the true value of θ. In the top image, θ is shown in the whole assimilation window t ∈ [0, 110]. In the bottom image, a zoom for t ∈ [0, 15] is shown. The shaded area represents the 95% credible interval for the case fa = 10.

  (a) t = 1.23 (b) t = 8.32 (c) t = 16.30

Figure 12 :

 12 Figure 12: Estimations by MEnKF of the momentum ρu normalized by ρ 0 u 0 for fa = 2 at t = 1.23 (a), t = 8.32 (b) and t = 16.30 (c). Times are given in tc units. The grey shaded area corresponds to the observation window.

  (a) t = 1.23 (b) t = 8.32 (c) t = 16.30

Figure 13 :

 13 Figure 13: Estimations by MEnKF of the momentum ρu normalized by ρ 0 u 0 for fa = 10 at t = 1.23 (a), t = 8.32 (b) and t = 16.30 (c). Times are given in tc units. The grey shaded area corresponds to the observation window.

  (a) t = 1.23 (b) t = 8.32 (c) t = 16.30

Figure 14 :

 14 Figure 14: Estimations by MEnKF of the momentum ρu normalized by ρ 0 u 0 for fa = 55 at t = 1.23 (a), t = 8.32 (b) and t = 16.30 (c). Times are given in tc units. The grey shaded area corresponds to the observation window.

Figure 15 :

 15 Figure 15: Time evolution of the RMS error of ρu for fa = 2, 10, 55.

(a) 1

 1 : constant (b) 1: time-varying

Figure 16 :

 16 Figure 16: Visualization of the normal momentum ρv (S.I. units) for the 2D compressible Navier-Stokes equation. Reference simulation at t = 10 for a constant value of 1 (a), and a time-varying value (b).

  , which are run for a total simulation time of T ref = 40 in t c units. A fully developed state obtained from a prior simulation at t = 10 is used to initialize the simulations at t = 0. Data are sampled every 30 time steps in the region x ∈ [0, 0.55] and y ∈ [-0.16, 0.16]. Considering the results obtained in Sec. 5, the update frequency is chosen sufficiently high to assure a good estimation. The observations are made from the instantaneous fields ρu and ρv. The data used as observation are artificially perturbed using a Gaussian noise of variance R = 1.The estimation algorithm is run over a time window equal to T DA = 40 which encompass roughly 3200 DA analysis phases. This value corresponds to four complete advections in the whole physical domain. At the end of each analysis, the mean value and the variance of the coefficient 1 are updated following the Dual EnKF technique[START_REF] Hamid | Dual state-parameter estimation of hydrological models using ensemble kalman filter[END_REF], similarly to what was done in Sec. 4.First, results dealing with the constant value of 1 are discussed. The parameter estimation for this case is almost exact, as observed in Fig.17

  . The convergence towards the exact value of 1 is obtained in roughly 0.5 units of time, after a first transient where a slight over-prediction is observed. When combined with the correction obtained via the Kalman gain, the accuracy of the state estimation is excellent, as shown in Fig.18.

Figure 17 :

 17 Figure 17: Time evolution of the inferred values of 1 . The truth corresponds to 1 = cte = 0.15. In the zoomed region, the shaded area represents the 95% credible interval for the estimated parameter.

Figure 18 :

 18 Figure 18: Estimations obtained by MEnKF of the momentum ρv (S.I. units) at the centerline y = 0 of the mixing layer. Results at t = 1 (a), t = 5 (b) and t = 30 (c) for the case 1 = cte = 0.15.

Figure 19 :

 19 Figure 19: Time evolution of the inferred values of 1 for the time-varying reference case. (a) Large time window. (b) Zoomed region. The shaded area represents the 95% credible interval for the estimated parameter.

Figure 20 :

 20 Figure 20: Estimations obtained by MEnKF of the momentum ρv (S.I. units) at the centerline y = 0 of the mixing layer. Results at t = 1 (a), t = 5 (b) and t = 30 (c) for the time-varying 1 .

Figure 21 :

 21 Figure 21: Estimations obtained by MEnKF of the momentum ρv (S.I. units) at the centerline y = 0 of the mixing layer. Here, MEnKF is only used to provide the estimation of 1 . Results at t = 1 (a), t = 5 (b) and t = 30 (c) for the time-varying 1 .

Figure 22 :

 22 Figure 22: Time evolution of the RMS error of ρv for the case of a time-varying inlet parameter 1 . The symbol P.E. corresponds to the case where MEnKF is only used for the estimation of 1 . The notation MEnKF corresponds to the standard version of the algorithm, including parameter estimation and physical state correction via Kalman gain.

  optimization problems. The compressible effects are relatively low for the tests performed so far. Further research is planned on test-cases where the compressible effects are more accentuated. = M k+1:k x a k P f k+1 = M k+1:k P a k M k+1:k + Q k+1 Forecast/Prediction k ← k+1

Figure A. 24 :Algorithm 1 : 1 : 2 : 3 : 5 :

 2411235 Figure A.24: Kalman Filter algorithm. The initialization is made with the analysed state.

: 8 :)Algorithm 2 :

 82 Update the ensemble; i = 1, • • • , N e Compute the ensemble forecast ; i = 1, • • • , N e x f,(i) k+1 = M k+1:k (x a,(i)k Dual Ensemble Kalman Filter (slightly adapted from[START_REF] Hamid | Dual state-parameter estimation of hydrological models using ensemble kalman filter[END_REF]). Use of anomaly matrices withY f k = H k X f k . We have i = 1, • • • , N e .Input: For k = 1, . . . , K: the forward models M k:k-1 , the observation models H k , the observation error covariance matrices R k Output: {θ a,(i) k } and {x a,(i) k } ; k = 0, • • • , K begin

Algorithm 3 :}

 3 Dual Ensemble Kalman filter of Algo. 2 applied on the coarse mesh. We have i = 1, • • • , N e . Input: For k = 1, . . . , K: the forward models M C k:k-1 , the observation models H C k , the observation error covariance matrices R C k Output: {θ a,(i) k } and {(x C k ) a,(i) } ; k = 0, • • • , K begin and {(x C 0 ) a,(i) }; for k = 1, . . . , K do 2:

k 3 :

 3 Dual EnKF on coarse mesh: Apply Algo. 3 if Observation available then

  43 = 64 times less mesh elements than the single simulation on the fine grid. If one considers that one main simulation and 100 ensemble members are run simultaneously, and if the RAM requirement is normalized over the main simulation, this implies that R RAM , the non-dimensional RAM requirement, is equal to 1 + 100/64 = 2.56. In other words, the total cost in RAM is increased to just 2.56 times the cost of the simulation without EnKF. For r C = 8, the normalized RAM requirement is R RAM = 1 + 100/8 3 = 1.195, thus just a 20% increase in RAM requirements. This is clearly orders of magnitude more advantageous than a fine-grid classical EnKF application with N e = 100, since in this case R RAM = N e = 100.

While the advancement model is unique (the Navier-Stokes equations, for example), the discretization process contained in Ψ is unique for each member of the ensemble. To distinguish them, it is therefore necessary to introduce an exponent i in the notations.

= cte = 0.15:

Acknowledgements: Our research activities are supported by the Direction Générale de l'Armement (DGA) and the Région Nouvelle Aquitaine. Prof. Heng Xiao is warmly acknowledged for valuable discussion on the subject.

The Reynolds number of the flow is calculated as Re = (U 1 -U 2 )δ 0 /ν = 100 with asymptotic velocities set to U 1 = 173.61 and U 2 = 104.17. These values correspond to a Mach number Ma = 0.5 and Ma = 0.3, for each stream, respectively. The kinematic viscosity and thermal diffusivity of the flow are considered to be constant and their value is fixed to ν = 1.568 × 10 -5 and α = 22.07 × 10 -7 , respectively. All these quantities are expressed in S.I. units.

The inlet boundary condition is taken from [START_REF] Ko | Sensitivity of two-dimensional spatially developing mixing layers with respect to uncertain inflow conditions[END_REF]. For the velocity field, one has:

V in = 0 [START_REF] Burgers | On the analysis scheme in the ensemble kalman filter[END_REF] where U in is the streamwise velocity at the inlet and V in is the normal velocity.

U in is estimated as a classical hyperbolic tangent profile plus a time-varying perturbation component:

where N in is the total number of perturbation modes and i quantifies the magnitude of each mode. The function f i (y) = cos(4n i y δ0 )h(y) controls the shape of the perturbation of the inlet velocity profile in the normal direction.

The role of h(y) = 1-tanh 2 ( 2y δ0 ) is to damp the perturbation component moving away from the centerline. The wavelength parameters n i are tuned according to the LST results. In the following, we consider N in = 1 i.e. the inlet perturbation consists of a single mode. The inlet density is set to be constant so that ρ in = 1.177, as well as the temperature T = 300 in S.I. units.

In this section, the parameters θ of the model correspond to i , the variable governing the amplitude of the perturbation. Two different scenarios are studied. In the first one, the parameter i is constant. In the second one, it is a time-dependent coefficient. In the first case, the reference simulation is done variables, the final solution is regularized.

The method, which is referred to as Multigrid Ensemble Kalman Filter (MEnKF), is assessed via the analysis of one-dimensional and two-dimensional test cases and using different dynamic equations. First, the one-dimensional Burgers equation for Re = 200 is analyzed. Here, the performance of MEnKF is assessed considering several coarsening ratios r C , which determines the difference in resolution between the main simulation and the ensemble members. The results for r C = 1 (i.e. method equivalent to a EnKF) indicate that the State Estimation and the parametric optimization of the inlet provide very high accuracy in the results. With increasing coarsening ratios the quality of the results is progressively degraded, but the main features of the flow are obtained even for very under-resolved ensemble members. In addition, higher r C values are associated with significantly decreased computational costs, so that this method exhibit a potential to be explored for efficient trade off between accuracy and resources required.

Then, MEnKF is used to track the time evolution of a free parameter for the case of a wave propagation, using a one-dimensional Euler model. Three cases are here investigated, varying the time window between successive assimilations.

The estimator can efficiently represent the evolution in time of the parameter, as well as to provide an accurate state estimation. However, the global prediction is significantly degraded if the assimilation window is larger than a threshold value, which is arguably connected to the physical features of the flow.

At last, the analysis of the two-dimensional spatially evolving mixing layer for Re = 100 is performed. The algorithm appears to be well suited for the analysis of unsteady phenomena, in particular for the analysis of time varying free parameters of the simulation. These features are promising for potential application to in-streaming Data Assimilation techniques.

Future research over this method will target improvement of the projector operators between coarse and fine grid level, in particular for the state matrix Ψ. Preliminary tests have shown that, in the case of non-linear models, information from the main, refined simulation can be projected on the coarse mesh Conflict of interest: The Authors have no conflict of interests.

Appendix A.2. Ensemble Kalman filter algorithm

An efficient implementation of the EnKF relying on anomaly matrices is given in Algo. 1. We have used the secant method described in [START_REF] Asch | Data Assimilation: methods, algorithms, and applications[END_REF] to change the definition of the variable Y f k .

Appendix A.3. Dual Ensemble Kalman filter algorithm

An efficient implementation of the Dual EnKF relying on anomaly matrices is given in Algo. 2. We have slightly adapted this algorithm from [START_REF] Hamid | Dual state-parameter estimation of hydrological models using ensemble kalman filter[END_REF].

Appendix A.4. Multigrid Ensemble Kalman filter algorithm

The algorithm 4 represents a simplified, ready-to-use application of the conceptual methodology presented in Sec. 3. This algorithm was tailored for the relatively simple physical models used in this work. While it may not be suited for complex three-dimensional applications, it proved an optimum trade-off in accuracy and computational resources for the present analysis.

First of all, when observation is not available, the two main forecast operations (fine grid forecast and ensemble coarse forecast) are performed using explicit time advancement schemes. This choice allows to reduce the computational costs. However, when observation is available, the following strategies are employed:

1. The two forecast operations (main simulation and ensemble members) are performed using an implicit matrix-splitting iterative procedure, using a single iteration. As stated in Sec. 3, the state transition model for each ensemble member is determined independently, but using the same structure and discretization schemes of the main simulation.

2. The number of iterative solutions for the main simulation on the coarsegrid level is equal to zero. That is, the solution from the first forecast is projected on the coarse grid, and the difference between the KF state estimation and this forecast is re-projected over the fine grid.

3. In the final iteration on the fine grid, an implicit matrix-splitting iterative procedure is employed, using a single iteration and a relaxation coefficient α = 0.5. This choice, which provides the best compromise between accuracy and regularization, has been identified after extensive tests for the configurations investigated.