G Artana

A Cammilleri

J Carlier

E Mémin

Strong and weak constraint variational assimilations for reduced order fluid flow modeling

Keywords: POD, variational assimilation, reduced order dynamical systems, PIV, wake Flow

In this work we propose and evaluate two variational data assimilation techniques for the estimation of low order surrogate experimental dynamical models for fluid flows. Both methods are built from optimal control recipes and rely on proper orthogonal decomposition and a Galerkin projection of the Navier Stokes equation. The techniques proposed differ in the control variables they involve. The first one introduces a weak dynamical model defined only up to an additional uncertainty time-dependent function whereas the second one, handles a strong dynamical constraint in which the dynamical system's coefficients constitute the control variables. Both choices correspond to different approximations of the relation between the reduced basis on which is expressed the motion field and the basis components that have been neglected in the reduced order model construction. The techniques have been assessed on numerical data and for real experimental conditions with noisy Image Velocimetry data.

Introduction

The reduction to a set of coupled ordinary differential equations of high dimensional dynamical systems associated to the Navier Stokes equations has been an intensive subject of studies since a long time. Despite the extreme temporal and spatial complexity of fluid flows, model reduction principles can often be applied to capture and represent the main characteristics of the flow dynamics. The construction of so-called low order dynamical systems (LODS) or reduced order models can be viewed as an approach in which one sacrifices the ambitious goal of having a full representation of the flow to a more modest one that consists in capturing and representing just the most representative components of the flow (with respect to a given criteria) through a small number of coupled ordinary differential equations. Such a representation is of major interest when one aims at characterizing the current dynamical and kinematical states of a given flow with a low computational load. This is crucial for instance in flow control applications.

This system of ordinary differential equations can be determined based on experimental results, so these kinds of models are also named empirical models. The snapshots of flow fields used as experimental data to generate the models may be issued either from numerical or physical experiments. Models produced from numerical simulations have been used in general as a benchmark to assess the model performances, since in that case all the data are available to validate the results obtained. On the other hand models provided by physical experiments carried out in laboratories environments have been settled in order to get closer to real world experiments and applications. Such models have however to face some additional complexities due to inherent measurements uncertainties, to the absence of data in some regions or to poor spatial and temporal resolution of the images. These reduced dynamical systems have been used for instance to improve numerical simulations in tuning the boundary conditions to more realistic situations [START_REF] Perret | Polynomial identification of POD based low-order dynamical system[END_REF].

One of the possibilities to neglect the influence of some degrees of freedom over the main characteristics of the flow consists in representing the flow as a sum of northonormal spatial basis functions weighted with temporal coefficients. To determine the basis function of this expansion, one of the usual approaches relies on the Karhunen Loeve decomposition (referred as proper orthogonal decomposition -POD -in the fluid mechanics domain). In this case the spatial basis functions, also called modes, and the associated temporal coefficients are computed from eigenvalues problems associated either to the two points correlation tensor or to the temporal correlation tensor (see for instance Holmes et al. [START_REF] Holmes | Turbulence, coherence structures, dynamical systems and symetry[END_REF] for a review). Beyond its easy implementation, this decomposition offers several advantages. It is optimal with respect to the L 2 -norm and the recovered modes are organized with respect to a decreasing mean kinetic energy order. Thus a truncation of the number of modes based on energetic criteria is straightforward. Also, as the spatial modes are orthonormal, a Galerkin projection of the Navier Stokes equations, under the hypothesis of incompressible flow, reduces the system of partial differential equations to a system of ordinary differential equations with respect to the temporal coefficients of the decomposition.

Unfortunately, these reduced order models determine in general a dynamical system that in many cases may converge to erroneous states after a relative short time of integration. Holmes et al. [START_REF] Holmes | Turbulence, coherence structures, dynamical systems and symetry[END_REF] have mentioned as potential source of this problem the following issues:

• Neglect of the incidence of boundary or pressure terms in the computation of the dynamical system coefficients,

• Inaccurate estimation of derivatives of spatial modes that determine the dynamical system coefficients,

• Problems related to low-dimensional truncation associated either to the dissipation loss induced by the truncated modes or to a generation of a non-realistic incomplete system.

Concerning the first point several authors have disregarded in their analysis the influence of the pressure term in wake flows arguing that this effect could be ignored when large wake domains were considered. However, depending on the problem or on the truncation level, neglecting this term may add some uncertainties in the ability of the computed dynamical system to faithfully represent the actual flow dynamics (Kalb and Deane [START_REF] Kalb | An intrinsic stabilization scheme for proper orthogonal decomposition based low-dimensional models[END_REF], Noack et al. [START_REF] Noack | The need for a pressure-term representation in empirical Galerkin models of incompressible shear flows[END_REF]). As a remedy, the introduction of correcting terms to take into account the effect in the outflow boundary through modifications of the linear [START_REF] Galetti | Accurate model reduction of transient and forced wakes[END_REF] or quadratic terms have been proposed [START_REF] Noack | The need for a pressure-term representation in empirical Galerkin models of incompressible shear flows[END_REF].

The second point refers to the fact that the coefficients of the reduced dynamical system involve spatial derivatives that can only be poorly determined on experimental data. To avoid the direct explicit computation of these coefficients it is possible to rely on regression techniques such as the one proposed in [START_REF] Perret | Polynomial identification of POD based low-order dynamical system[END_REF]. In this approach the system's coefficients are estimated within a least square goodness of fit from the sequence of temporal coefficients observed at the snapshot instants and a finite difference approximation of their temporal derivatives (assuming a very small sampling time to get an accurate derivative).

As for the third point, the effect of dissipation of the unresolved modes on the dynamics of the system has been proposed to be included through the addition of an artificial viscosity. The dissipation coefficients could be determined either by a proper modeling or through an appropriate tuning allowing obtaining qualitatively correct results. One of the classical solution consists in adding a constant viscosity acting in the same way on all the POD modes [START_REF] Delville | Examination of large scale structures in a turbulent mixing layer part.1[END_REF][START_REF] Podvin | A low-dimensional approach for the minimal flow unit[END_REF]. This artificial viscosity is added to the flow kinematic viscosity and enables after a proper tuning of its value to improve the system's numerical stability. Rempfer and Fasel [START_REF] Rempfer | Evolution of three-dimensional coherent structures in a flat-plate boundary layer[END_REF] and Rempfer [START_REF] Rempfer | Investigation of of boundary layer transition via Galerkin projection on empirirical eigenfunctions[END_REF] have proposed modal viscosities that affect differently the modes of the decomposition. More recently, Karamanos and Karniadakis [START_REF] Karamanos | A spectral vanishing viscosity method for large-eddy simulations[END_REF] employed a dissipative model called spectral vanishing viscosity model to formulate an alternative approach. Cazemier et al. [START_REF] Cazemier | Proper orthogonal decomposition and lowdimensional models for driven cavity flows[END_REF] used also modal kinetic equations to determine viscosities to be added on each POD mode. The main drawback of all these methods is that they require fixing a priori a great number of parameters.

Some authors have also proposed to enhance the system by considering larger basis than the one provided by a given snapshots sequence associated to a single case of the flow characteristic parameters. A better response of the system is expected by increasing the subspace spanned by the truncated dynamical system and the possibility to include a statistical variance attached to the attractor [START_REF] Bangia | Unsteady two-dimensional flows in complex geometries: Comparative bifurcation studies with global eigenfunction expansions[END_REF][START_REF] Jorgensen | Low-dimensional modeling of a driven cavity flow with two free parameters[END_REF]. However, while this may stabilize the scheme, the limit cycles amplitude obtained with such extended basis may be different from those obtained for the original basis [START_REF] Deane | Low-dimensional models for complex geometry fows: Application to grooved channels and circular cylinders[END_REF].

Besides, for the case of wake flows as pointed out by [START_REF] Noack | A hierarchy of low-dimensional models for the transient and post-transient cylinder wake[END_REF], the error sources are inherent to the reduced system as defined from a Galerkin projection. The authors highlighted that it is not possible to guarantee, a priori, the system structural stability as little perturbations on the system's coefficients can provide qualitatively very different solutions. In order to restore dissipation or nonlinear large-scales/small scales interactions to the reduced system some authors have proposed to use non-linear Galerkin techniques in which an explicit relation between the resolved and non-resolved modes together with a clear time scale separation principle between the associated temporal coefficients are hypothesized.

Nonlinear Galerkin projection techniques, which are based on the approximate inertial manifold concept [START_REF] Constantin | Integral Manifolds and Inertial. Manifolds for Dissipative Partial Differential Equations[END_REF][START_REF] Foias | Exponential tracking and approximation of inertial manifolds for dissipative equation[END_REF][START_REF] Marion | Nonlinear Galerkin methods[END_REF][START_REF] Margolin | The postprocessing Galerkin and nonlinear Galerkin methods: A truncation analysis point of view[END_REF] have been successfully used for the reconstruction of low order dynamical systems from experimental particle images measurements [START_REF] Ma | Dpiv-driven simulation: a new computational paradigm[END_REF]. Despite some sensibility to the initial condition such approaches have shown to provide long term stability for discrete systems [START_REF] Shen | Long time stability and convergence for fully discrete nonlinear Galerkin methods[END_REF] even if some authors report situations for which such frameworks have revealed to be inadequate [START_REF] Sirisup | A spectral viscosity method for correcting the long-term behavior of POD models[END_REF]. It must be also outlined that, in the same way as for the classical Galerkin projection, non-linear Galerkin methods require to have an accurate estimation of the coefficients of the original ODE system.

It is also possible to try to improve the accuracy and stability of the dynamical system itself without directly addressing the error sources. In this spirit, Kalb and Deane [START_REF] Kalb | An intrinsic stabilization scheme for proper orthogonal decomposition based low-dimensional models[END_REF] proposed an intrinsic stabilization scheme that enables to make corrections on the constant and linear terms of the LODS Galerkin system. With this formulation the dynamical system is reinitialized at each snapshot instant to the temporal eigenfunctions value, and time-stepped to the next snapshot time so that time integration of the dynamical system does not diverge. The model seems more adapted to numerical data since for experimental data the "true value" of the temporal eigenfunctions are unavailable due to the presence of noise caused by experimental conditions.

The method considered in this paper makes use of an alternative solution for building a reduced order flow representation from the data and to cope with an inaccurate specification of the system's coefficients or of the modes truncation closure. The approach proposed relies on data assimilation principles built from optimal control theory [START_REF] Lions | Optimal Control of Systems Governed by Partial Differential Equations[END_REF] that are intensively used in meteorology and oceanography sciences for the analysis of geophysical flows [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF]. In the domain of wake flows other authors have used approaches in the same spirit to improve the accuracy of the POD-Galerkin models [START_REF] Bergmann | Optimal control of the cylinder wake in the laminar regime by trustregion methods and POD reduced-order models[END_REF][START_REF] Couplet | Calibrated reduced-orer POD-Galerkin system for fluid flow modeling[END_REF][START_REF] Galetti | Low-order modelling of laminar flow regimes past a confined square cylinder[END_REF][START_REF] Galetti | Accurate model reduction of transient and forced wakes[END_REF]. These methods introduce some calibration terms, which may include an eddy viscosity term or additional constant or linear terms. Those terms are computed through the minimization of a goodness of fit cost function between the POD temporal coefficients and those predicted by the calibrated model under the constraint imposed by the calibrated dynamics. The method proposed here can be regarded as a calibration technique on the system's coefficients as well. It is nevertheless formulated slightly differently as an optimal control problem with an explicit constraint on the initial condition and on the system's coefficients and relies on general norms weighted by appropriate covariance matrices. The method extends a similar technique we proposed recently [START_REF] Adamo | Variational assimilation of POD low-order dynamical systems[END_REF]. This technique was formulated also as a weakly constrained optimal control problem incorporating an additive uncertainty function [START_REF] Adamo | Variational assimilation of POD low-order dynamical systems[END_REF]. We will show that considering for the control variables an additional uncertainty function or the dynamical system's coefficients corresponds indeed to different assumptions on the reduced dynamical model. Both methods have their own advantages and we will show how they can be nicely combined. These techniques will be compared on the basis of experimental and numerical data. Let us mention that POD reduced order models have been used to build efficient variational assimilation strategies of oceanic models [START_REF] Cao | A reduced order approach to four-dimensional variational data assimilation using proper orthogonal decomposition[END_REF][START_REF] Fang | Reduced order modelling of an adaptive mesh ocean model[END_REF]. However, in those approaches, only the initial condition is considered as a control variable. In these studies, the authors aimed at lowering the computational cost of the forward and backward integration of the full dynamical model required within a variational assimilation process. The coefficients of the reduced model are thus directly set from samples of the dynamics forward integration; they are not estimated from noisy experimental data as we propose to do it in this study.

The paper is structured as follow. In section 2, we briefly recall the principles of the Proper Orthogonal Decomposition and those governing the obtention of a reduced order modeling through a linear Galerkin projection. After pointing out the different limitations associated either to a direct modeling of the reduced system or to a polynomial identification of the dynamics from experimental data, we present in section 3 the main concepts of the variational data assimilation techniques we used in this work. Section 4 describes the successive application of two variational assimilation strategies to the definition of reduced order dynamical models from noisy velocity measurements. The last section is dedicated to numerical and experimental results of a turbulent wake of a 4 cylinder at Reynolds 3 900.

POD-Galerkin model

The POD-Galerkin method has been extensively described in many articles. The method foundations are for instance thoroughly described in [START_REF] Holmes | Turbulence, coherence structures, dynamical systems and symetry[END_REF]. Let us review very briefly its main characteristics. As previously presented, a flow field u(x, t) ∈ L 2 (Ω), where Ω is a physical domain with boundary ∂Ω, is decomposed on a finite set of a n-orthonormal spatial basis functions {φ i (x) : Ω ⊂ IR 2 → IR 2 , i = 1, . . . , n}, with temporal coefficients {a i (t) : IR + → IR, i = 1, . . . , n} such that:

u(x, t) = n i=1 a i (t)φ i (x). ( 1
)
The spatial modes φ(x) are global and mutually orthogonal. They are usually extracted through an eigenvalue problem from a set of observations u(x, t i ) obtained experimentally or numerically. When working on finite dimension, their identification reduces to a SVD decomposition of a matrix that gathers in its columns the different snapshots of the velocity fields u(x, t i ). The obtained modes are optimal with respect to the L 2 (Ω) scalar product:

f, g = Ω f (x)g(x)dx,
whose associated norm corresponds to the kinetic energy of the flow contained in Ω. The spatial modes also inherit by construction of properties of the flow that can be expressed as linear homogeneous equations of the velocity fields (incompressibility condition for instance). The Navier Stokes equations, and the incompressibility condition for the fluctuating components are then reduced by means of a Galerkin projection (with respect to the L 2 (Ω) inner product) onto the subspace spanned by the first n vectors of the POD basis:

∂u ∂t + (u • ∇)u + ∇p - 1 Re ∆u, φ j = 0. (2)
As outlined previously the incompressibility condition is implicitly imposed to the basis functions through the data and hence can be drop. Homogeneous Dirichlet or Neumann boundary conditions are also implicitly satisfied by the basis functions. Let us also recall that the POD representation is only defined for zero mean signals. The temporal average must be thus either incorporated as the first mode of the decomposition (associated to a unity coefficient) or considered as an additive component (that is previously removed from the data). For sake of simplicity, in this study we will not distinguish between these two cases that have obviously direct implications on the form of the Navier-Sokes equations considered (i.e. primitive form vs. Reynolds average form). Let us note that removing the mean field from the data has the advantage to enable coping with spatially non-homogeneous stationary boundary conditions and provides a simple so-called lifting procedure. For a non-stationary boundary condition a similar approach can be envisaged by subtracting a non-stationary velocity field respecting this condition.

Identifying u(x, t) to its modal decomposition in equation ( 2), leads to a quadratic system of ordinary differential equations describing the evolution of the temporal coefficients:

da k dt = f k + n i=1 l ik a i + n i=1 n j=1 a i b ijk a j , k = 1, • • • , n, (3)
with :

f k = Ω ∇p(x)φ k (x)dx l ik = - 1 Re Ω ∆φ i φ k dx + ϕ ik , (4)
b ijk = Ω φ j ∇φ i φ k dx,
and where the terms ϕ ij and p introduced in the constant and linear coefficient respectively depends on the pressure function. In order to define further their expressions, we proceed to a decomposition of the pressure term as proposed in [START_REF] Boree | Exended proper orthogonal decomposition: a tool to analyse correlated events in turbulent flows[END_REF]:

p(x, t) = p(x) + n k=1 a k (t)ψ k (x) + p d (x, t). (5)
In this decomposition the first term represents the pressure temporal average. The second expression represents the pressure term belonging to the subspace spanned by the truncated temporal modes basis. The last term consists of the remaining pressure. It depends on the POD basis components but also on the modes that have been neglected in the basis truncation stage. This decomposition introduces an extended POD basis with new basis functions defined as:

ψ i (x) = 1 σ i T to a i (t)p(x, t)dt, (6)
σ i = T to a 2 i (t)dt. (7)
From this decomposition, the projection of the pressure gradient on the POD basis components may be expressed as:

Ω ∇p(x, t)φ i (x)dx = Ω ∇pφ i dx + n k=1 a k (t)ϕ ki + ξ i (t), (8)
ϕ ki = Ω ∇ψ k φ i dx, (9)
ξ i (t) = Ω ∇p d (x, t)φ i dx. (10)
The second term of this expression corresponds to the influence of the pressure gradient projected on the subspace spanned by the POD temporal modes basis components. The third term represents the contribution of the rest of the fluctuating pressure components.

For incompressible flows the pressure terms can be further evaluated on the boundary of the domain after integration by part. In some cases such as wake flows for instance, considering a uniform velocity on the inflow boundary and the far downstream outflow boundaries, the contribution of these terms is negligible or small. In the general case however, they do not cancel and are completely unknown. A direct calculation of the dynamic system's coefficients requires a modeling of those terms and to restrict strictly their computation to numerical data. The estimation of the system's coefficients from the data enables to integrate implicitly the contribution of those terms in the estimated coefficients given of course that their influences are observable on the available measurements. The last pressure term, ξ i (t), cannot be included in the constant, linear or quadratic coefficients and should be added to the reduced dynamics as an additional time varying component. However, as this term is at least partly caused by the unresolved modes it is quite difficult to model. In this work we will show how a particular form of an optimal control strategy allows us specifying from the data the contribution of such a function. Before describing the data assimilation framework on which we rely to estimate the system coefficients, we describe first the effects of the low-dimensional truncation procedure through a nonlinear Galerkin technique.

Unresolved modes modeling with nonlinear Galerkin methods

Nonlinear Galerkin methods have been proposed by several authors in order to take into account the non-resolved components in flow numerical simulations. Those techniques separates the solution of a dynamical system on two subspaces related respectively to slow and fast modes or to large and small scales. The small scales dynamics are simplified and assumed to follow linear or stationary pde's. This is referred as a slaving principle in which the small scales/fast modes are slaved to the evolution of the slow/large scales components of the velocity. An exact (when possible) or approximate algebraic relation modeling this slaving is injected in the slow modes components in order to close the large scales dynamics. It is expected that such modeling bring a better approximation of the unresolved components than the linear Galerkin projection, which consists to neglect such a relation and to cancel the influence of fast modes.

In order to describe further those methods, let us consider an infinite expansion of the velocities on a L 2 (Ω) POD basis, where we identify the velocity u(x, t) with an infinite dimension vector of the temporal coefficient a = (a 1 , a 2 , . . . , a ∞ ). Decomposing this vector by defining two projections P : L 2 (Ω) → H n and Q = I -P : L 2 (Ω) → H r respectively defined as the projection on the space spanned by the truncated POD basis {φ 1 (x), . . . , φ n (x)} (i.e. the resolved modes) and the complement space spanned by {φ n+1 (x), . . . , φ ∞ (x)} (i.e. the unresolved modes) , we get

P (u) = p = (a 1 , a 2 , . . . , a n ), ( 11
) Q(u) = q = (a n+1 , a n+2 , . . . , a ∞ ). ( 12
)
When these projections are associated to a spectral decomposition of a linear dynamical operator (i.e the dynamics linear tangent operator for instance) those methods are referred as centered manifold (if there exists a set of pure imaginary eigenvalues associated to an eigenspace called the center eigenspace) or otherwise as (approximate) inertial manifold [START_REF] Carr | Applications of center manifold theory[END_REF][START_REF] Wiggins | Introduction to applied nonlinear dynamical systems and chaos[END_REF]. Let us note that inertial manifolds have been proved to exist for many dissipative systems but their existence remains open for the Navier-Stokes equations.

Considering the evolution model ( 3) for a = p + q written in compact form as

da dt = La + B(a, a) + f, (13)
where L, B and f denote respectively linear, bilinear and constant operators, the projections P and Q lead to the system:

dp dt = L n (p + q) + B n (p + q, p + q) + f n , (14)
dq dt = L ∞ (p + q) + B ∞ (p + q, p + q) + f ∞ , (15)
where the superscript n and ∞ refer to projections P and Q respectively. Assuming that the unresolved modes evolves on a much smaller time scale compared to the resolved modes, and considering that the resolved modes concentrate a higher energy than the unresolved ones (or equivalently the temporal p-modes exhibit much higher amplitudes than the q ones), the following approximate resolution of the second equation can be proposed [START_REF] Aling | Nonlinear model reduction for simulation and control of rapid thermal processing[END_REF][START_REF] Deane | Low-dimensional models for complex geometry fows: Application to grooved channels and circular cylinders[END_REF]:

q = ϕ(p). (16)
This expression relates, in a generic way, univocally the non-resolved modes to the sole resolved modes. Injecting ϕ : H n → H r into the resolved modes dynamics ( 14) enables to close the problem:

dp dt = L n (p + ϕ(p)) + B n (p + ϕ(p), p + ϕ(p)) + f n . ( 17
)
As a simple example considering a linear (eventually time dependent) relation q = a + M (t)p leads for the resolved modes to a dynamical system:

dp dt = L n (I n + M (t))p + B n (a, p) + B n (p, a) + B n (a, M (t)p) + B n (M (t)p, a) Ln (t)p + B n (p + M (t)p, p + M (t)p) Bn t (p,p) + f n + L n a + B(a, a) f n , (18)
which remains of the same form but depends now on new coefficients involving eventually time-dependent terms.

Considering invariant manifold reductions through a linear stability analysis around a base flow state, it can be proposed that the resolved/unresolved modes relations can be expressed as a second order polynomial [START_REF] Noack | A hierarchy of low-dimensional models for the transient and post-transient cylinder wake[END_REF] . In this case the associated extended dynamical system includes third and forth order terms. However, these two additional terms correspond to unresolved modes interactions, which are associated to low amount of the kinetic energy. Compared to the other terms involving the most energetic modes, it is far likely that these terms will have only a small influence on the short-term dynamics of the resolved modes. Thus, it is customary to simplify the system by sticking to a quadratic dynamical system.

Let us point out that a linear Galerkin projection associated to the traditional POD reduced order modeling is obviously recovered through a null dependency (independence assumption of the two modes sets). Compared to this solution, we see that a simple linear relation between the resolved and unresolved modes involves a dynamical system in which the system's coefficients are modified and depend on this interaction model. It can be expected that a sound modeling of this interaction enables to reproduce more faithfully the system's dynamics. However, this gain in precision is general paid by an increase of the order of the dynamical system, which impacts generally its stability. Instead of an accurate modeling of this interaction, in this work we will rather rely on procedures allowing an estimation of the dynamics' coefficients from the data. The credo we will follow consists to hypothesize that the data reflect in a way or another the influence of the truncated unresolved modes. An accurate and robust estimation method should therefore enable to propagate implicitly the relation between the two modes sets from the data to the coefficient values. To that end, we propose in this work to explore the use of two different optimal control techniques for the assimilation of measurement data into the reduced order dynamical systems and the estimation of its initial condition and of its coefficients. The first technique will consider a dynamics in which the coefficients L and B constitutes stationary control variables of the assimilation process. The second strategy will incorporate an additive time dependent uncertainty variable, ϑ(t), to the classical POD reduced order dynamics. This corresponds to a case in which this additional additive control variable gathers the time-dependent interaction terms:

ϑ(t) = L n M (t)p + B n (a, M (t)p) + B n (M (t)p, M (t)p) + B n (p, M (t)p) + B n (M (t)p, p). ( 19
)
The first method will allow us eventually to forecast future states of the system whereas the second one is limited to analysis tasks on the time interval in which observation measurements are available. We will show in particular that the combination of both strategies brought some very good results for the constitution of a reduced order model within the time period corresponding to the observation time range. Before entering precisely into the details of both methods we describe the general context of the variational data assimilation techniques.

Data assimilation

The framework considered in this paper relies on optimal control theory introduced by Lions [START_REF] Lions | Optimal Control of Systems Governed by Partial Differential Equations[END_REF]. It is intensively used in meteorology and oceanography in order to perform data assimilation of geophysical flows [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF]. In the following we present the general principles governing such variational data assimilation schemes. Even though we will use such principles in a particular context, we wish to describe this methodology in a general setup so that interested reader may eventually adapt the same derivation to different objectives.

The problem we are dealing with consists in recovering a system's state X(x, t) obeying to a dynamical law given some noisy and possibly incomplete measurements of the state. The measurements (also called observations) are assumed to be available only at discrete time. This is formalized, for any location x ∈ Ω at time t ∈ [t 0 , t f ], by the system:

∂tX(x, t) + M(X(x, t), u(t)) = 0 (20) X(x, t0) = X0(x) + ǫ(x), ( 21
)
where M is a nonlinear dynamical operator depending on a control parameter u(t). We assume here that u(t) ∈ U and X(t) ∈ V are square integrable functions in Hilbert spaces identified to their dual. The term X0 is the initial vector at time t 0 and ǫ is an (unknown) zero mean additive control variable on the initial condition. Besides, we assume that measurements of the unknown state Y(t) ∈ O are available. These observations are measured through the nonlinear operator H and belong also to an Hilbert space:

Y(x, t) = HX(x, t) + η(x), ( 22
)
where η is a zero mean Gaussian random field.

Cost functional

From this system, data assimilation aims at recovering an initial condition and a value of the dynamics' control parameter that leads to the lowest discrepancy between the measurements and the state variable. This objective can be formalized as the minimization of a cost functional, J : U × V → R, defined as:

J (u, ǫ) = 1 2 t f t0 Y(t) -H(X(u(t), ǫ, t)) 2 R dt + 1 2 ǫ 2 Ic + 1 2 t f t0 u(t) -u 0 2 F dt. ( 23
)
This functional of the initial condition and the dynamics' control parameter gathers three terms. The first term comes directly from the measurement equation. It is a quadratic best fit term between the observation and the state variable provided by the dynamics integration for a given value of the functional parameters. The second term aims at specifying a low error on the initial condition whereas the third term enforces the control variable to be close to a given a priori value u 0 of the control parameter. For a null value of this a priori a control of lowest energy is sought. The norms corresponds to the Mahalanobis distance defined from the inner products < R -1 ., . >O, < I -1 c ., . >V and < F -1 ., . >U of the measurements, the state variable and the control variable spaces respectively. They involve covariance tensors R, I c and F related to the measurement error, the error on the initial condition and the deviation between the control and its a priori value. In our applications, these covariance tensors have been defined as diagonal tensor (i.e. the noise are assumed to be uncorrelated in time and space). For example the observation covariance tensor has been set to a covariance tensor of the form:

R(x, t, x ′ , t ′ ) = σδ(x -x ′ )δ(t -t ′ ), (24)
and similar expressions hold for F and I c . In order to compute the gradient of this functional we assume that X(u(t), ǫ; t) depends continuously on (u(t), ǫ) and is differentiable with respect to the control variables u(t) and ǫ, on the whole time range.

Differentiation

Noting first that dX = (∂X/∂u)δu(t) + (∂X/∂ǫ)δǫ, the differentiation of equations [START_REF] Kalb | An intrinsic stabilization scheme for proper orthogonal decomposition based low-dimensional models[END_REF][START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF] in the direction (δu, δǫ) reads:

∂ t dX + ∂ X M(X, u(t))dX + ∂ u M(X, u(t))δu(t) = 0, ( 25
) dX(x, t 0 ) = δǫ(x), (26)
where ∂ X M and ∂ u M denotes the linear tangent operators defined by:

lim β→0 M(X + βX, u(t)) -M(X, u(t)) β = ∂ X M(X)dX. ( 27
)
The differentiation of the cost function [START_REF] Ma | Dpiv-driven simulation: a new computational paradigm[END_REF] in the direction (δu, δǫ) (denoting U T as the space of square integrable function on a spatio-temporal domain) reads then:

∂J ∂u , δu U T = t f t0 u(t)-u 0 , δu(t) F dt- t f t0 Y(t)-H(X(t)), (∂ X H) ( ∂X ∂u δu(t)) O dt, (28)
∂J ∂ǫ , δǫ V = (X(x, t 0 ) -X 0 (x)), δǫ Ic - t f t0 Y(t) -H(X(t)), (∂ X H) ( ∂X ∂ǫ δǫ) O dt.
(29) Introducing the adjoint of the linear tangent operator (∂ X H) * , defined as:

∀(x, y) ∈ (V, O), < (∂ X H) x, y > O =< x, (∂ X H) * y > V , (30)
these two relations can be reformulated as:

∂J ∂u , δu U T = t f t0 F -1 (u(t) -u 0 ), δu(t) U dt- t f t0 (∂ X H) * R -1 (Y(t) -H(X(t)), ∂X ∂u δu(t) V dt, (31)
and ∂J ∂ǫ , δǫ

V = I -1 c (X(x, t 0 ) -X 0 (x)), δǫ V - t f t0 (∂ X H) * R -1 (Y(t) -H(X(t)), ∂X ∂ǫ δǫ V dt. ( 32
)
Expression [START_REF] Perret | Polynomial identification of POD based low-order dynamical system[END_REF][START_REF] Podvin | A low-dimensional approach for the minimal flow unit[END_REF] provides the functional gradients in the directions (δu, δǫ m ). We can remark from these expressions that a direct numerical evaluation of these gradients is in practice completely unfeasible. As a matter of fact, such an evaluation would require to compute perturbations of the state variable along all the components of the control variables (δu, δǫ)i.e. integrate the dynamical model for all perturbed components of the control variables, which is computationally completely unrealistic.

Adjoint model

An elegant solution of this problem consists in relying on an adjoint formulation [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF][START_REF] Lions | Optimal Control of Systems Governed by Partial Differential Equations[END_REF]. To that end, the integration over the range [t 0 , t f ] of the inner product between an adjoint variable λ ∈ V T and relation ( 25) is performed:

t f t0 ∂dX ∂t (t), λ(t) V dt + t f t0 (∂ X M) dX(t), λ(t) V dt + t f t0 (∂ u M) δu(t), λ(t) V dt = 0.
An integration by parts of the first term yields:

-

t f t0 - ∂λ ∂t (t)) + (∂ X M) * λ(t), dX(t) V dt = λ(t f ), dX(t f ) V -λ(t 0 ), dX(t 0 ) V + t f t0 δu(t), (∂ u M) * λ(t) U dt, ( 33
)
where the adjoint of the tangent linear operators (∂ X M) * : V → V and (∂ u M) * : V → U have been introduced. At this point no particular assumptions nor constraints have been imposed on the adjoint variable. However, we are free to particularize the set of adjoint variables of interest in setting a particular evolution equation or a given boundary conditions allowing simplifying the computation of the functional gradient. As we will see it, imposing that the adjoint variable λ is solution of the system:

-∂ t λ(t) + (∂ X M) * λ(t) = (∂ X H) * R -1 (Y -H(X(t))) λ(t f ) = 0, (34)
will provide us a simple and accessible solution for the functional gradient.

As a matter of fact, injecting this relation into equation [START_REF] Rempfer | Investigation of of boundary layer transition via Galerkin projection on empirirical eigenfunctions[END_REF] with dX(t0) = δǫ and dX = (∂X/∂u)δu(t) + (∂X/∂ǫ)δǫ allows identifying the right hand second terms of the functional gradients [START_REF] Perret | Polynomial identification of POD based low-order dynamical system[END_REF][START_REF] Podvin | A low-dimensional approach for the minimal flow unit[END_REF] and we get

∂J ∂ǫ , δǫ V = -λ(t 0 ), δǫ V + I -1 c (X(t 0 ) -X 0 ), δǫ V , ∂J ∂u , δu U T = t f t0 δu(t), F -1 (u(t) -u 0 ) + (∂ u M) * λ(t) U dt = F -1 (u -u 0 )) + (∂ u M) * λ, δu U T .
From these relations, one can now readily identify the two components of the cost function derivatives with respect to the control variables:

∂J ∂ǫ = -λ(t 0 ) + I -1 c (X(t 0 ) -X 0 ), ∂J ∂u = F -1 (u -u 0 ) + (∂ u M) * λ. (35)
The partial derivatives of J are now simple to compute when the adjoint variable λ is available. The knowledge of the functional gradient enables then to define updating 1. Set an initial condition: X(t 0 ) = X 0 2. From X(t 0 ), compute X(t) with the forward integration of relation (20) 3. Compute the adjoint variable λ(t) with the backward integration of relation (34) 4. Update the initial value X(t 0 ) and the parameter model u with (36) 5. Loop to step 2 until convergence rules for the control variables from iterative optimization procedures. A quasi-Newton minimization process consists for instance of:

X n+1 (t 0 ) = X n (t 0 ) -α n H-1 Xn(t0) (I -1 c (X n (t 0 ) -X 0 ) -λ(t 0 )), u n+1 = u n -α n H-1 un ((∂ un M) * λ + F -1 (u n -u 0 )), (36)
where H-1 xn denotes an approximation of the Hessian inverse computed from the functional gradient with respect to variable x n ; the constant α n is chosen so that to respect Wolfe conditions. The adjoint variable is accessible through a forward integration of the state dynamics [START_REF] Kalb | An intrinsic stabilization scheme for proper orthogonal decomposition based low-dimensional models[END_REF][START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF] and a backward integration of the adjoint variable dynamics [START_REF] Rempfer | Evolution of three-dimensional coherent structures in a flat-plate boundary layer[END_REF]. Let us point out that considering a final condition for the state variable (through a similar cost function term as for the initial condition) would change the null initial condition of the adjoint dynamics into a term similar to the one involved in the derivative with respect to the initial condition control variable. The overall optimal control process is schematically summarized in figure 1.

Such a formulation provides us a practical framework for the estimation of state variables trajectory from a sequence of noisy and possibly incomplete observation variables. The variables used to control the system may be either parameters or global uncertainty functions of the dynamical model. Both of these possibilities will be exploited in the following to settle two assimilation schemes for low order dynamical systems.

Application to Low order dynamical systems

We apply here the assimilation formalism for the estimation of low order dynamical systems. The first method, which corresponds to a refinement of a method detailed in [START_REF] Adamo | Variational assimilation of POD low-order dynamical systems[END_REF], consists in performing an assimilation of the temporal modes considering an additive uncertainty function on the reduced dynamical model. The coefficients of the dynamical system are here first assumed to be directly provided by the data. These coefficients are then updated, in a second time, from the assimilated temporal modes trajectory. In the second technique, these coefficients are directly assimilated jointly with the temporal modes. The uncertainty is in that case put on each parameter of the model and not anymore through an additional uncertainty function. We will see that this latter techniques, brings a significant improvement in the accuracy of the low order dynamical system.

Recovery of the temporal modes trajectory with a weak constraint assimilation

This first technique corresponds to a continuation of the work developed in [START_REF] Adamo | Variational assimilation of POD low-order dynamical systems[END_REF]. For this assimilation problem, we consider the dynamics of the temporal POD modes defined in equations [START_REF] Bergmann | Optimal control of the cylinder wake in the laminar regime by trustregion methods and POD reduced-order models[END_REF][START_REF] Boree | Exended proper orthogonal decomposition: a tool to analyse correlated events in turbulent flows[END_REF] and augmented by an additive uncertainty function, ω(t) ≃ ξ(t) + ϑ(t), encoding both the contribution of the component of the fluctuating pressure belonging to the complement space of the subspace spanned by the POD temporal modes, ξ(t) [START_REF] Adamo | Variational assimilation of POD low-order dynamical systems[END_REF] and the interaction between the modes of the reduced POD basis and the secondary modes that have been truncated, ϑ(t) [START_REF] Jorgensen | Low-dimensional modeling of a driven cavity flow with two free parameters[END_REF]:

∂ t a(t) = f + La(t) + B(a(t), a(t)) + ω(t), a(t 0 ) = a 0 + η. ( 37
)
In this system the constant, linear and bilinear operators of the dynamics, f , L and B are assumed to be known and fixed.They can be provided for instance by a regression technique from the data [START_REF] Perret | Polynomial identification of POD based low-order dynamical system[END_REF]. Functions ω(t) and η represent zero mean uncertainty functions with covariances Q and I c associated to the dynamical model and to the initial condition respectively. They are assumed to belong to Hilbert spaces W and I respectively. The aim of the assimilation consists to control these two functions in order to fit at best a sequence of data a obs ∈ O and to recover the temporal modes a ∈ A([t 0 , t f ]), where the space of the temporal modes trajectory is defined as

A([t 0 , t f ]) = g/g ∈ L 2 ([t 0 , t f ], ∂ t g ∈ L 2 ([t 0 , t f ]) .
The measurements a obs of the temporal modes are assumed to be provided through a standard POD analysis (Snapshots method [START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF] or direct method followed by a projection of the velocity measurements onto the spatial modes). Under Gaussian hypothesis of the model error (i.e.

[ da(t) -E da(t)|a(0)] ∼ N (0, Q))1 , the background error (i.e. η ∼ N (0, I c )) and the observation error (i.e. [a obs (t)a(t)|a(t)] ∼ N (0, R)), finding a control of lower mean square energy associated to the smallest discrepancy between the measurement and the unknown modes can be formalized as the minimization of the following functional:

J(ω, η) = 1 2 t f t0 a obs (t) -a(ω, η, t) 2 R dt + 1 2 t f t0 ω(t) 2 Q dt + 1 2 η 2 Ic . (38)
This functional corresponds to a so-called weak constraint assimilation problem as this minimization problem is equivalent to a weakly constrained optimization in which the dynamics augmented by an uncertainty function constitutes the optimization constraint. This uncertainty function constitutes an additive time dependant control variable on the dynamics. The third term of the functional consists to fix, up to a Gaussian noise, the initial condition to an a priori value. This value can be provided by the measurement at the initial time or fixed through a previous data assimilation on a time interval ending at t o . Let us now derive the expression of the associated functional gradient as given through the introduction of an adjoint variable driven by the adjoint of the tangent linear operator associated to the considered dynamical model.

Linear tangent operator

The linear tangent operator of the considered dynamical model is obtained computing (∂ a M) * for a small perturbation θ = (θ 1 , . . . , θ s ) T :

(∂ a M) (θ k ) = -   n i=1 l ik θ i + ij a i b ijk θ j + θ i b ijk a j   , k = 1, . . . , n. (39)
The adjoint of this operator is very similar and includes the adjoint L * of matrix L ℓk = l ℓk and the adjoint of matrices B(a, •) and B(•, a) where B(a, •) kℓ = i a i b iℓk and B(•, a) kℓ = i a i b ℓik respectively. The dynamics associated to the adjoint variables λ(t) derived in a general case in the previous section [START_REF] Rempfer | Evolution of three-dimensional coherent structures in a flat-plate boundary layer[END_REF] reads here as:

-∂ t λ(t) + L * λ(t) + B * (a(t), λ(t)) + B * (λ(t), a(t)) = R -1 (a obs (t) -a(t)), λ(t f ) = 0. ( 40
)
The knowledge of the adjoint variable provides us the functional gradients with respect to the control variable, which are given by the system:

∂ η J = -λ(t 0 ) + I -1 c (a(t 0 ) -a 0 ), ∂ ω J = Q -1 (∂ t a(t) -f -La(t) -B(a(t), a(t))) -λ, (41)
where the control function ω has been replaced by its expression: ω = ∂ t a(t)-f -La(t)-B(a(t), a(t)). The whole assimilation process consists schematically to integrate forward the system (37) for a null initial uncertainty value. The adjoint system is then integrated backward from a null terminal condition. The trajectory of this adjoint variable enables to get the gradient of the functional with respect to the two control variables (41). These gradients can then be used as descent directions of a minimization procedure. In this work we used the limited storage variant of the BFGS quasi-Newton optimization technique [START_REF] Nocedal | Numerical optimization[END_REF]. A complete synopsis of the whole algorithm is summarized in figure 2. This synopsis corresponds to an improved version of the technique proposed in [START_REF] Adamo | Variational assimilation of POD low-order dynamical systems[END_REF], where a fixed step gradient descent was originally proposed. The supplementary degree of freedom brought by the uncertainty variable has shown to provide very good results [START_REF] Adamo | Variational assimilation of POD low-order dynamical systems[END_REF]. It allows us adapting the modes trajectory to noisy data. Nevertheless as indicated previously, in this approach, the dynamical system coefficients are considered as known and given through a least square regression technique. These coefficients can be subsequently updated from the estimated assimilated trajectory. This constitutes thus finally a two steps method that is not completely satisfactory, as those coefficients have not been directly considered as variables of the problem. As a matter of fact, it must be outlined that the coefficients first obtained through a regression technique on noisy data do not allow usually to reconstruct a complete trajectory of the dynamical system on the whole sequence of data. To that end, an adhoc damping of the dynamics was considered to stabilize the system and to get a first run of the assimilation process [START_REF] Adamo | Variational assimilation of POD low-order dynamical systems[END_REF]. This first run is all the more important that if it is too far away from the observation, the backward adjoint integration may fail to successfully correct the modes dynamics, which yields a local convergence of the method toward a meaningless local minima.

In order to improve the accuracy and the robustness of the dynamical model reconstruction, we propose in the following a modified assimilation scheme in which the POD-ROM Weak Constraint Assimilation (WCA) Algorithm

Initializations :

• Compute the POD of a set of velocity snapshots {u(•, ti), ti = 0, . . . , f };

• Fix the POD-ROM dynamical coefficients (through a least square estimation);

• Fix a null uncertainty ω n = 0;

• Fix an initial condition a 0 (to) = a obs (to) = ( u(to), φi , i = 1, . . . , k) T ;

• n = 0; 2. Minimization procedure (a) Integrate forward the dynamics from a n (to):

∂ta n (t) = f + La n (t) + B(a n (t), a n (t)) + ω n (t) (b) Integrate backward the adjoint dynamics from λ(t f ) = 0 -∂tλ(t) + L * λ(t) + B * (a n (t), λ(t)) + B * (λ(t), a n (t)) = R -1 (a obs (t) -a n (t))
(c) Quasi-Newton descent iteration (LBFGS [START_REF] Nocedal | Numerical optimization[END_REF])

a n+1 (to) = a n (to) -α n a H-1 a n (to) I -1 c (a n (to) -a obs (to)) -λ(to) , ω n+1 = ω n -α n ω H-1 ω n Q -1 ω n -λ , (42)
where H-1

x n denotes the LBFGS inverse Hessian approximation computed from the functional gradient with respect to variable x n ; the constant α n

x is chosen with a line search procedure so that to respect Wolfe conditions. 3. Loop to step 2 until convergence.

results: initial condition a(t o ) and uncertainty function ω(t)

Figure 2: Synopsis of the Weak constraint assimilation algorithm system coefficients are now defined as control variables of the system. This new assimilation process will allow us to properly estimate the dynamics' coefficients and to initialize efficiently a weak constraint assimilation.

Assimilation of the low order dynamical system with a strong constraint assimilation

We are now considering an assimilation strategy defined from a control on the initial condition and on the coefficient f , L and B(, ) of the dynamical system. Apart from the initial condition, the dynamical model we are considering here does not include anymore any uncertainty terms:

∂ t a(t) + f + La(t) + B(a(t), a(t)) = 0, a(t 0 ) = a 0 + η. ( 43
)
Opposite to the previous assimilation this dynamical model will constitute a strong constraint on the associated optimization problem. More precisely, this scheme consists to seek a control that minimizes a cost functional encoding the variances of the measurement discrepancy and the unknown coefficients subject to the constraint given by the dynamics (43). The measurement and the unknown coefficients are here assumed to be distributed according to Gaussian laws centered on a(t), f 0 , L 0 and b 0 = (b 0 ijk , 0 ≤ i, j, k ≤ n) respectively. The corresponding cost function reads hence:

H(f, L, b, η) = t f t0 1 2 a obs (t) -a(f, L, B, η, t) 2 R dt+ 1 2σ 2 f f -f 0 2 + 1 2σ 2 L L -L 0 2 F + 1 2σ 2 b b -b 0 2 F + 1 2 η 2 Ic . (44)
The Frobenius norm is here considered for the matrix L and the tensor b. This model differs from the previous one, as it introduces three Gaussian uncertainties (of constant variances) on the parameter values instead of a global additive uncertainty model on the temporal modes. The expected value of the parameters, (or the mean of the uncertainties) f o , L o , b o can be fixed either through a standard polynomial identification technique from the POD modes [START_REF] Perret | Polynomial identification of POD based low-order dynamical system[END_REF] or fixed to a priori values. As for the different variances, in the case of a polynomial identification they can be set to the least squares estimation residues.

The adjoint system corresponding to this new problem is unchanged [START_REF] Wiggins | Introduction to applied nonlinear dynamical systems and chaos[END_REF]. The functional gradient with respect to the control variables, I, L and b, reads now:

         ∂ η H = -λ(t o ) + I -1 c (a(t 0 ) -a obs (0)), ∂ f H = σ 2 f (f -f o ) - t f to λ(t)dt, ∂ L H = σ 2 L (L -L o ) - t f to a(t)λ(t)dt, ∂ b H = σ 2 b (b -b o ) - t f to a(t)a(t) T ⊗ λ(t)dt. ( 45
)
As previously, the assimilation of the low order dynamical system consists to integrate the dynamical system (43) starting from the initial measurements, and initial values of the coefficients f , L and b. The initial value of the temporal coefficients are assumed to be given by the measured coefficients, a obs (0). As for the system coefficient's a priori values, in our experiments, we systematically set them to a zero value as we did not want to rely on an external estimation process. A backward integration of the adjoint system Type enables then to compute the functional gradient, which provides a descent direction for the cost functional minimization. The updated values of the control variables are then used to integrate again the direct dynamical model ( 43) and so on until convergence. The whole algorithm is summarized in figure 3. Opposite to the previous scheme, this new assimilation system has the advantage to authorize the direct estimation of the low order dynamical system and the recovering of the restored modes trajectory in a single step. We will see in the following experimental section that this unified estimation process provides good results on a difficult case. We will furthermore show that the combination of both assimilation processes allows us to estimate very precisely the trajectory of the POD-Galerkin system.

N ∆t T [vs] L x × L y n x × n y PIV 2D -2C

Data benchmarks for experimental validation

In order to assess the performances of the two data assimilation methods for the specification of low order dynamical models we have constituted a benchmark composed of large eddy simulation results and experimental PIV data. The numerical simulation and the experimental data concerns both a wake behind a circular cylinder at Reynolds 3 900. The experimental setup and the large eddies simulation technique used for this wake flow are described in Appendix A and Appendix B respectively. Below we detail the data benchmarks corresponding to these two cases.

Time sequences used

For the PIV, the whole sequence constituted of 3 072 successive 2D -2C velocity fields has been selected for the first benchmark. For the LES, a sequence of only 300 full 3D -3C velocity fields with a 20∆t time step between them has been kept for building a snapshots sequence of the flow motion. Velocity fields similar to the PIV motion fields have been extracted from this database. The time-resolved sequence is made up of 300 2D -2C velocity fields just behind the circular cylinder in the plane z = 0 with a square size of 3D × 3D and with a number of approximately 80 snapshots by vortex shedding.

Table 1 summarizes the main characteristics of the velocity fields estimated with the PIV technique and those obtained from LES. Few general remarks can be outlined. The spatial and temporal resolutions are much higher for LES than for PIV. The temporal length of the sequence is much greater for PIV than for LES. These characteristics reflect well the usual respective advantages of experimental measurements (free Reynolds number, ability to generate long time sequence at low cost) and numerical simulation approaches (no noise, full components in space with a high spatial density).

Figure 4 illustrates two snapshots of the spanwize vorticity ω y (adimensionned with U e and ν) in the plane z = 0 for LES and PIV. Both fields present similar vorticity magnitude and features. The spatial resolution of the PIV measurements (32 × 32 pixels final window size and 50% overlapping) acts as a low-pass filter but the noise level is POD-ROM Strong Constraint Assimilation (SCA) Algorithm • Fix an initial condition a 0 (to) = a obs (to) = ( u(to), φi , i = 1, . . . , k) T ;

• n = 0;

2. Minimization procedure (a) Integrate forward the dynamics from a n (to):

∂ta n (t) = f n + L n a n (t) + B n (a n (t), a n (t)) (b) Integrate backward the adjoint dynamics from λ(t f ) = 0 -∂tλ(t) + L n * λ(t) + B n * (a n (t), λ(t)) + B n * (λ(t), a n (t)) = R -1 (a obs (t) -a n (t))
(c) Quasi-Newton descent iteration (LBFGS [START_REF] Nocedal | Numerical optimization[END_REF])

a n+1 (to) = a n (to) -α n a H-1 a n (to) I -1 c (a n (to) -a obs (to)) -λ(to) , f n+1 = f n -α n f H-1 f n σ 2 f (f n -f o ) - t f to λ(t)dt , L n+1 = L n -α n L H-1 L n σ 2 L (L n -L o ) - t f to a(t)λ(t)dt , b n+1 = b n -α n b H-1 b n σ 2 b (b n -b o ) - t f to a(t)a(t) T ⊗ λ(t)dt , (46)
where H-1

x n denotes the LBFGS inverse Hessian approximation computed from the functional gradient with respect to variable x n ; the constant α n

x is fixed with a line search procedure so that to respect Wolfe conditions. 3. Loop to step 2 until convergence.

results: initial condition a(t o ) and dynamics coefficients f , L and b. emphasized by the derivative scheme (centered 2nd order). The spanwise vorticity ω z highlights the two shear layers on both sides of the cylinder. For this Reynolds number (Re D = 3 900), the boundary layer on the body is laminar and transition to turbulence takes place in the shear layers. While one of the two shear layers is almost straight and laminar up to x/D = 2, the second one becomes unstable and rolls up to form a primary vortex. Small vortices resulting from the instability of the shear layer seem to accumulate in this primary large scale structure. This scenario is observed alternatively in the two shear layers and at different stages of their evolution in the instantaneous fields. The Bénard-von Karman vortex street begins further downstream. The 3D instantaneous flow fields, which are not illustrated here, show clearly the fully three-dimensional character of the dynamics. As expected from the flow's Reynolds number these fields exhibits a wide range of scales. The vortex street, which is the main structure of the flow, is quasi-2D with shedding parallel to the cylinder. Numerous filament vortices can be also clearly observed.

Proper Orthogonal Decomposition

The velocity field sequences have been preprocessed removing the mean part (as the Reynolds decomposition) and considering an adimensionalization of the fluctuating part with respect to the external velocity U e . This preprocessing is useful to remove the mean kinetic energy part and to facilitate the comparison between LES and PIV (and for readers who might be eventually interested in using these data). A Proper Orthogonal Decomposition is performed with the snapshot method [START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF] using the full sequences of the PIV measurements (3 072 snapshots) and of the LES results (300 snapshots).

Figure 5 presents for the LES and PIV sequences the distribution of the turbulent kinetic energy with respect to the number of modes. Figure 6 shows the cumulative distribution of this energy. The slope of the spectra changes beyond 30 modes. This seems to indicate that the first spatial modes are probably converged. The PIV sequence and the LES sequences possess close energy levels for the first modes. From these figures we can remark that the first two modes are the principal modes of the flow. They represent 62% of the total turbulent kinetic energy for the PIV measurements and 60% for the LES sequence. These coupled modes correspond to the large vortices of the von Karman street and to their advection downstream. The following modes have a much lower contribution. The first 32 modes represent 84% (+22% compared to the two first modes) of the total turbulent kinetic energy for PIV and 94% (+34%) for LES. This relatively low decreasing can be explained by the relatively "high" Reynolds number of the flow in this study. The fastest energy decrease of the LES is due to a lower variability of the coherent structures present in a shorter number of snapshots. In this experimental study as we focus on the ability of the assimilation techniques to faithfully represent the flow dynamics along a given time range from the corresponding available data, questions of statistical convergence of the data are not primordial. In particular, we do not aim at comparing precisely the results obtained on these two sequences but rather to show that the assimilation works either for noise free numerical data or for noisy and sparse experimental data. 0.5 The two variational assimilation strategies (the strong and weak constraint assimilation) have been applied in order to assess their ability to accurately reproduce the trajectory of a reduced order dynamical systems with 32 temporal modes from discrete sequences of POD temporal modes provided by the snapshots method [START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF].

The Gaussian uncertainties σ f , σ L and σ b of the dynamical coefficients have been arbitrarily identically fixed to 10 3 . The covariance tensors R(x, x ′ ) and I c (x, x ′ ) have been defined as the diagonal tensors 1 2σ R δ(xx ′ ) and 1 2σ I δ(xx ′ ) with σ R = 1 and σ I = 10 3 .

As for the experiments, the strong constraint assimilation strategy has been first performed on three sequences of POD temporal modes. The first sequence is a LES sequence depicting 3 vortex shedding and corresponding to 256 successive snapshots. The other sequences are two PIV sequences showing 3 and 12 vortex shedding, corresponding respectively to 128 and 512 successive snapshots. The target values for the initial condition a 0 is given by the POD data whereas the target values for f 0 , L 0 and b 0 have been fixed to zero values. In this way, no previous calibration method (such as a polynomial identification associated to a damped dynamics) is used to set the values of these coefficients. At the very first iteration of the assimilation, the initial values of the dynamics' coefficients L, b are fixed to zero as well and the initial conditions of the unknown modes variables are initialized to the observed modes value at the initial time.

The weak constraint assimilation strategy is then run on the same sequence of temporal modes after convergence of the strong constraint assimilation. The dynamical coefficients and the initial conditions of the modes are set to those estimated from the previous assimilation. The target values for a 0 are not changed and are given by the POD data. At the first iteration, the uncertainty function is fixed to a null value.

In the following section we show first the results obtained from the strong constraint assimilation, noted as (SCA). The trajectories corresponding to the weak constraint assimilation, noted (WCA) are presented in a subsequent section.

Strong Constraint Assimilation

Figures 7, 8 and 9 present respectively the temporal modes trajectories of the first 32 modes recovered after the assimilation of the observed POD temporal modes.

As can be observed from this curves the strong constraint assimilation technique enables for the short time LES and PIV sequences to build a reduced order models that faithfully recover with a good accuracy the trajectories corresponding to the first modes. For the longer PIV sequence, the less energetic modes trajectories are less accurately recovered (Fig. 9). However the mean tendency of the modes are well represented. It can be outlined, that the corresponding reduced order models have been obtained considering a null initialization. In addition, the targeted values of the dynamical coefficients being fixed to zero values the assimilation consequently favors coefficients with a low magnitude. Despite this zero initialization and non informative targeted values of the dynamical coefficients, the assimilation succeeds to converge to a meaningful solution for the first modes of the dynamics. To illustrate the velocity fields reconstructed from the reduced dynamical model estimated through the strong constraint assimilation procedure, we plot in figure 10

Weak Constraint Assimilation

The coefficients of the dynamical system estimated from the previous strong constraint assimilation are now injected within the weak constraint assimilation technique. The initial condition of this new assimilation is initialized with the initial condition estimated from the previous assimilation scheme. It is however important to outline that this initial state remains a control variable of the weak constraint assimilation. The additional control variable consists in an uncertainty function representing the missing part of the dynamics. As previously stated from a simplified nonlinear Galerkin projection this additional variable can be interpreted as a term gathering the time dependent interaction terms between secondary modes and modes of the truncated basis (see section 2.1) and also as the eventual contribution of a component of the fluctuating pressure (see equation 10 and the text below this equation). The Figures 11,[START_REF] Delville | Examination of large scale structures in a turbulent mixing layer part.1[END_REF] and 13 present the results obtained with the weak constraint assimilation strategy. In these figures the trajectories of the dynamics complemented with an uncertainty variable are plotted. As can be observed from these results the trajectories fit perfectly the data in the three cases. The uncertainty control variable enables to complete the missing elements of the reduced dynamics to explain the data. This uncertainty function is however only defined for the time range on which the data are available. To define a reduced model enabling to forecast future states of the system, it could be interesting to rely on a decomposition of this uncertainty signal. Assuming the uncertainty function is stationary in a weak sense, the covariance signal matrix could be learned for each mode from this assimilation on the past data. The future states of this uncertainty function could be then adequately randomly sampled conditioned on its past trajectory. This study, which is beyond the scope of the present paper, would have the advantage to allow defining a reduced forecast model that could be useful for flow control prospect. This will be the route we intend to follow in a future work.

In order to emphasize the importance of the initial condition as a control variable, which is usually not considered in POD calibration techniques [START_REF] Bergmann | Optimal control of the cylinder wake in the laminar regime by trustregion methods and POD reduced-order models[END_REF][START_REF] Couplet | Calibrated reduced-orer POD-Galerkin system for fluid flow modeling[END_REF][START_REF] Galetti | Accurate model reduction of transient and forced wakes[END_REF][START_REF] Galetti | Low-order modelling of laminar flow regimes past a confined square cylinder[END_REF], we plot in figure 14 the modes trajectories obtained from the POD initial condition associated to the noise free LES sequence. We see clearly on the trajectory plots that the simulation diverges after the second vortex shedding. The initial states ensuing from the proper orthogonal decomposition are clearly not adapted to the estimated quadratic reduced dynamical model. The effect of phenomenon related to higher order dynamical terms or to 3D effects are probably too pregnant in this initial condition. As alternatives we can consider either to build a more accurate higher order polynomial dynamical model (of third or forth order) with eventual inherent numerical stability problems, or to adapt the initial condition to the kind of selected approximation of the Navier-Stokes POD-Galerkin projection. This latter solution was implemented in the data assimilation strategy adopted. The figure 15 illustrates however that the dependence on initial condition of the assimilated reduced dynamics can be slightly relaxed. In this figure we plot the modes trajectory obtained from an initial condition slightly perturbed in the direction of the POD data at the initial time. This initial condition is defined as: â(0) + 0.1(a obs oâ(0)), where â(0) refers to the initial condition estimated from the weak constraint data assimilation. The modes trajectories differ slightly from the trajectories starting from the assimilated initial condition but do not diverge. This shows a possible adaptation of the assimilated reduced order model to a small set of initial condition around the ideal assimilated initial modes values.

-0. 15: Weak constraint reduced order model forward integration from the assimilated initial conditions slightly perturbed in the direction of the POD initial condition associated to the LES sequence (dashed line). Reduced order model forward integration from the initial condition estimated with the weak constraint assimilation procedure (solid line).

CONCLUSION

From the benchmark of the flow around a cylinder problem at Reynolds 3 900, we proposed a method to improve the construction of a reduced order model from PIV measures, and numerical data (LES). As experimental data are prone to be affected by many errors (measurement noise due to 3D effects, measurements on a limited scales range, etc.) the estimation from regression techniques of reduced order models yields dynamical systems that often exhibit stability problems. Furthermore, those techniques do not allow providing an initial state of the different modes adapted to the chosen representation of the Navier-Stokes Galerkin projection. This is quite problematic when one aims, like in this study, at representing the quasi-2D component of the flow whereas the initial experimental condition includes strong 3D effects, important noise or inaccurate velocity measurements. To improve the deficiency of polynomial identification approaches we developed a variational data-assimilation scheme in which the initial state of the modes and the reduced order models coefficients constitute control variables of the system. The assimilation is conducted on the basis of POD modes estimated from a snapshots sequence of the flow velocity fields. This strategy enables the building of reduced order models initialized with a very uninformative prior (null initialization and null constant a priori values). This reconstructed dynamical model has shown to recover quite well the dynamics of a large amount of modes for short numerical and experimental sequences depicting three vortex shedding. For longer PIV sequences the last modes are less faithfully reconstructed but the mean tendency of the modes trajectories are still well described. The coupling of the strong constraint assimilation strategy with a weak constraint assimilation process enables to correct this deficiency by introducing an uncertainty control variable on the dynamics. Through a simplified nonlinear Galerkin projection it has been shown that such uncertainty variable can be interpreted as a time varying term of the dynamics representing both the coupling between the neglected secondary modes and the components of the POD basis and also the eventual contribution of the projection of the fluctuating pressure gradient component on the POD temporal modes. This assimilation procedure has enabled us building reduced order dynamical systems that fit quite perfectly the data for short or long temporal sequences of observed POD temporal modes. We believe this uncertainty function could be useful to design reduced order models for flows dynamics that are only partially known. The question of how to forecast such a reduced order dynamical system stays open, but we believe that the analysis of this uncertainty function could serve as a basis to generate for instance adequate random samples mimicking the behavior along time of the secondary modes.

Appendix A Wake behind a cylinder: experimental setup

The wake behind a circular cylinder at Reynolds 3 900 has been generated in the Cemagref wind tunnel. This fifteen meters long wind tunnel is usually used to study plane mixing layers where two parallel flows with different velocities and temperatures give rise to a turbulent mixing. For this study, this configuration has been slightly modified by an extension of the downstream separating plate. This modification provides two independent testing zones with identical dimension (a rectangular cross section of 100 × 50 cm 2 large and 300 cm long). The PIV experiments have been carried out in the upper testing zone for which the upper wall is slightly tilted to suppress the longitudinal 34 pressure gradient. The flow velocity can be chosen continuously between 1 and 5 m/s with a good stability. The uniformity of the velocity profile at the entrance of the test section has been checked by hot-wire anemometry. The free stream turbulence intensity was less than 0.5% and the temperature was kept within 0.2 • c by using an air-water heat-exchanger. The circular cylinder has been mounted horizontally at the entrance of the upper testing zone. It has been equipped with two thin rectangular end-plates to decrease the effects of the boundary layer which develops on the wall. The dimensions of the plates follow the specification recommended by Stansby [START_REF] Stansby | The effects of end plates on the base pressure coefficient of a circular cylinder[END_REF]. The distance between the end-plates is L = 89.6 cm and the diameter of the circular cylinder is D = 32 mm. The aspect ratio is L/D = 20 and the blockage ratio is 3.2%. The reference frame is direct and the origin is located at the center of the cylinder. The longitudinal x-axis is the testing zone axis of the wind tunnel, y-axis is normal to the wake central plane (symmetry plane) and z-axis is the cylinder axis. The velocity components u, v and w are related to x, y and z axis respectively.

The free-stream velocity has been adjusted at 1.82 m/s so that the Reynolds number be 3 900 (Reynolds number based on the free stream velocity Ue, the kinematic viscosity ν and the diameter of the circular cylinder D). For this Reynolds number, the flow is relatively well documented in the literature (Parnaudeau et al. [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3 900[END_REF]). The Strouhal number St is 0.208. The Strouhal frequency f s (vortex shedding frequency of the von Karman street) is then about 12 hz.

Time-resolved 2D-2C PIV experiments (two in-plane velocity components in a plane field) have been carried out with Nd-YAG laser (NewWave Pegasus with nominal energy by pulse of 10 mj at 3 000 hz) and a CMOS camera (Photron APX-RS with 50 -3 000 hz acquisition frequency range for 1 024 × 1 024 px image size). The laser sheet has been produced by means of a Rodenstock telescope and a cylindrical lens. The diameter of the particle seeding (diluted polyglycol in water) is about 1 µm.

The measurement area was located just behind the circular cylinder, in the plane z = 0. The camera was located perpendicularly to the laser sheet at a distance of 60 cm. For the mounted lens, the focal length was 50 mm with an aperture of 5.6. The resulting field of view is of dimension 3D × 3D. One sequence of 3 072 successive image pairs (equivalent to the 8 Go buffer memory of the fast camera) has been obtained with a 500 hz acquisition rate (≈ 40 snapshots by vortex shedding). The time interval between two pulses is 250 µs.

The PIV records have been analyzed through a cross-correlation technique implemented with a Fast-Fourier-Transform algorithm in a multi-grid process with 5 iterations and shift sub-pixel (1 × 128 × 128, 1 × 64 × 64 and 3 × 32 × 32 pixels) and 50% overlapping. The cross-correlation peaks were fit with a Gaussian function on 3 pixels. Erroneous velocity vectors have been identified by a median filter and replaced by their local mean value. To give an idea of the corresponding experimental accuracy, the particle displacement range is about -1 to +5 pixels (with an error less than 0.2 pixel on the instantaneous displacement) and less than 0.1% of the velocity vectors are erroneous.

Figure 1 :

 1 Figure 1: Schematic representation of the variational data-assimilation algorithm

•

  Compute the POD of a set of velocity snapshots {u(•, ti), ti = 0, . . . , f }; • Fix the POD-ROM dynamical coefficients to zero (or to a priori values f o , L o , b o if any);

Figure 3 :

 3 Figure 3: Synopsis of the Strong constraint assimilation (SCA) algorithm

Figure 4 :

 4 Figure 4: Snapshots of the vorticity field obtained from a Large Eddies Simulation (top) and estimated from a PIV technique (bottom).

Figure 5 :

 5 Figure 5: Turbulent kinetic energy distribution with respect to the number of modes for the 300-snapshots LES sequence (red curve) and the 3 072-snapshots PIV sequence (green curve).

Figure 6 :

 6 Figure 6: Cumulative turbulent kinetic energy distribution with respect to the number of modes for the 300-snapshots LES sequence (red curve) and the 3 072-snapshots PIV sequence (green curve).

32 Figure 7 :

 327 Figure 7: Reconstructed temporal modes trajectories for the LES sequence showing 3 vortex shedding (256 snapshots): dashed line, the data (every 2 snapshots instants); solid line, results of the strong constraint assimilation.

32 Figure 8 :

 328 Figure 8: Reconstructed temporal modes trajectories for the PIV sequence showing 3 vortex shedding (128 snapshots): dashed line, the data (every 2 snapshots instants); solid line, results of the strong constraint assimilation.

32 Figure 9 :Figure 10 :

 32910 Figure 9: Reconstructed temporal modes trajectories for the PIV sequence showing 12 vortex shedding (512 snapshots): dashed line, the data (every 4 snapshots instants); solid line, results of the strong constraint assimilation.

Table 1 :

 1 Characteristics of the available velocity fields from PIV and LES (vs for vortex shedding)

	3 072 0.11D/U e	72	3D × 3D 64 × 64
	LES 2D -2C 300	0.06D/U e	4	3D × 3D 144 × 144

  Reconstructed temporal modes trajectories for the LES sequence showing 3 vortex shedding (128 snapshots): dashed line, the data (every 2 snapshots instants); solid line, results of the weak constraint assimilation. Reconstructed temporal modes trajectories for the PIV sequence showing 3 vortex shedding (256 snapshots): dashed line, the data (every 2 snapshots instants); solid line, results of the weak constraint assimilation. Reconstructed temporal modes trajectories for the PIV sequence showing 12 vortex shedding (512 snapshots): dashed line, the data (every 4 snapshots instants); solid line, results of the weak constraint assimilation. Weak constraint reduced order model forward integration from the observed POD initial conditions of the LES sequence (dashed line). Reduced order model forward integration from the initial condition estimated with the weak constraint assimilation procedure (solid line)

	a 32 a 32 a 32 a 32 a 32	-0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.2 0.1
	a 30 a 30 a 30 a 30 a 30	-0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.2 0.1
	a 28 a 28 a 28 a 28 a 28	-0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.2 0.1
	a 12 a 12 a 12 a 12 a 12	-0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.2 0.1
	a 10 a 10 a 10 a 10 a 10	-0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.2 0.1
	a 8 a 8 a 8 a 8 a 8	-0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.1 0.2 -0.2 -0.1 0 0.2 0.1
	a 4 a 4 a 4 a 4 a 4	-0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.4 0.2
	a 3 a 3 a 3 a 3 a 3	-0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.4 0.2
	a 2 a 2 a 2 a 2 a 2	-0.8 -0.4 0 0.4 0.8 -0.8 -0.4 0 0.4 0.8 -0.8 -0.4 0 0.4 0.8 -0.8 -0.4 0 0.4 0.8 -0.8 -0.4 0 0.8 0.4
	8 -0.4 0 0.4 0.8 Figure 11: 29 0 2 4 6 a 1 -0.8 -0.4 0 0.4 0.8 0 2 4 6 a 1 Figure 12: -0.8 -0.4 0 0.4 0.8 0 10 20 a 1 Figure 13: -0.8 -0.4 0 0.4 0.8 0 2 4 6 a 1 0.8 0.4 1 0 a -0.4 0 2 4 6 Figure 14: -0.8 Figure	t t 30 t t t	8 8 8 8	10 10 40 10 10	12 12 50 12 12	14 14 14 14	60

Let us note this Gaussian assumption is formally true only for linear dynamics.

Appendix B Wake behind a cylinder: Large Eddy Simulation

The code used for this study solved the incompressible Navier-Stokes equations:

where the forcing term f represents the Immersed Boundary Method used to impose the no-slip condition at the cylinder surface. The principle of the present forcing is to create an internal flow preserving the no-slip boundary at the cylinder surface while avoiding strong discontinuities to prevent spurious oscillations. Basically, Large Eddy Simulation method consists in separating, through a spatial filtering operation, the great scales and the small scales of the turbulence. Spatial filtering (denoted by the overbar) of the Navier-Stokes equations is:

where T ij = u i u ju i u j is the subgrid-scale stress. According to the Boussinesq hypothesis, the subgrid-scale stress can be defined as:

where S is the shear stress tensor. Here, the turbulent viscosity ν t ∼ ∆ c ν ∆ is evaluated by the Structure Function model proposed by Métais and Lesieur [START_REF] Métais | Spectral large-eddy simulation of isotropic and stably stratified turbulence[END_REF], where ∆ c is the filter width estimated as the cubic root of the mesh volume and ν ∆ is the subfilter length scale.

The incompressible Navier-Stokes equations are solved on a regular Cartesian grid in non-staggered configuration. Sixth-order compact centered difference schemes are used to evaluate all spatial derivatives, except at the outlet and outflow boundaries where single sided schemes were employed for the x-derivative calculation. Time integration is performed with the second-order Adams-Bashforth scheme. A constant flow is imposed at the entrance of the domain and a simple convection equation is applied at the exit. Periodic conditions are used in the two transverse direction y and z.

Our computational domain extends over 20D in the streamwise and normal directions. The center of the cylinder is located at x cyl = 5D downstream of the inflow. The spanwize extent of the domain was chosen to be L z = πD, which corresponds to the size used by most previous authors. For the square which contains the cylinder, we only used 48 × 48 points in the streamwise and normal directions. A uniform streamwize flow is imposed in the initial conditions, with no perturbation in the other directions. For the regime considered, the boundary layer is laminar. The simulation was carried out with a constant time step size of ∆t = 0.003D/U e , which ensured that the Courant number be 0.15. The Reynolds number is Re D = 3 900, the domain size is L x × L y × L z = 20D × 20D × πD and the corresponding number of point is n x × n y × n z = 961 × 960 × 48. LES was performed for a time duration of ≃ 52 vortex shedding cycles to ensure a convergence of the 2nd order statistics.