Firstly, I thank Pr. Jörn Sesterhenn and Pr. Olivier Thual for the time they gave me to examine this thesis. I also thank the rest of my thesis committee : Jocelyn Erhel, Alejandro Gronskis, and Roger Lewandowski, for their encouragement, insightful comments, and interesting questions.

I would like to acknowledge my gratitude to my Ph.D supervisor Etienne Mémin for his guidance and I am very gratful for his welcome to the FLUMINANCE team. I thank him for offering me two wonderful opportunities to travel abroad to work with the department of fluid dynamics in Buenos Aires.

I would also like to thank my supervisor Dominique Heitz for his constant support and our insightful discussions over the entire course of this Ph.D. He showed one of the best examples of leadership and humanism I have ever seen.

I warmly thank Alejandro Gronskis for sharing generously his expertise and helping me to challenge and improve myself. I also thank Guillermo Artana for his warm welcome in his team at the department of fluid dynamics in Buenos Aires.

I express my gratitude to Sylvain Laizet for his expertise and his availability despite the distance. His programming skills and rigor pushed me to improve my own skills.

I thank Charles Deltel and Tina for their precious patience and their technical support. They were always available and willing to help whenever I ask.

I thank my co-worker Yin Yang for our friendly exchanges during this Ph.D. I thank him for our countless discussions and his valuable suggestions. It was a real pleasure to work with him for three years and to I thank Ioana Barbu for our balcony breaks and for sharing her interesting perspectives in our discussions. I thank her for moral support during the hard times. And I thank all the members of the FLUMINANCE team for the years we spent together.

I would like to thank sincerely all my dear friends from Toulouse who supported me despite the distance and their own issues, including Nadia, Amy, Rémi, Raphael, Jimmy, Matt, Julien and Will. I also thank my partner's family for their support and kind thoughts over the years. And finally, I owe my deepest gratitude to my partner Stéphan who constantly supported me and never failed to cheer me up even I was sometimes about to give up.

I dedicate this thesis to my parents Olivier and Ioly Robinson. 3 The classical PIV system configuration setup is composed of a doublepulsed laser and CCD camera. We deduce the velocity field by performing the analysis of the correlation between a sequence of frames. 1.4 The tomographic PIV system configuration setup is composed of a double-pulsed laser and CCD camera. We deduce the velocity field by performing the analysis of the correlation between a sequence of frames.

Table of contents

2.1 Sequential data assimilation representation : we seek for the analysis state x a k+1 that corrects the forecast state x f k+1 generated by the dynamics model given the observation Y k+1 at the instant t k+1 . The dynamics model start from the analysis state x a k of a previous assimilation process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Ensemble Kalman Filter representation : we consider a set, called ensemble, of possible state samples, each referred to as an ensemble member, noted {x f,(1)

k+1 , • • • , x f,(N )
k+1 } (blue dots). The latter is generated from the randomization of the non linear dynamics model given the result x a k of a previous assimilation. The analysis ensemble is obtained by correcting the forecast ensemble (red dots) with the observation Y k+1 (green dot) and is associated with the analysis covariance error matrix P e (red ellipse). . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Standard 4DVar assimilation representation : we aim at determining the optimal assimilation trajectory (red curve) from the background trajectory (blue curve) and a set of observations (green). We minimize the distance between the background trajectory and the observations by a gradient-descent algorithm. The cost function gradient is deduced by performing an adjoint procedure backwards in time (orange curve) which takes into account the observations correction. . . . . . . . . . 2.4 Incremental variational assimilation representation : we seek for the optimal increment δx

(1) 0 of the linearized convex cost function around the initial state x

(1) 0 = x b 0 (blue curve) and unpdate to x

(2) 0 = x

(1) 0 + δx (1) 0 . We repeat the procedure until we reach the optimal state x 0 which minimizes the nonlinear cost function J (red curve). . . . . . . 2.5 Construction of the differentiation code df(x) associated to the analytic function y = f (x) by the discretize-then-differentiate (solid line) and the differentiate-then-discretize (dashed line) strategies. . . . . . 56

3.1 Spatial discretization of the FVM grid : the cell C ij = [x i-1/2 ; x i+1/2 ]× [y j-1/2 ; y j+1/2 ] is the cell associated to the quantity x i,j . We aim at computing the the flux crossing each inter-cell flux F n i+1/2,j , F n i-1/2,j , G n i,j+1/2 and G n i,j-1/2 in each direction east, west, south and north respectively. . 4.11 Assimilation with 3D synthetic observations : comparison of the RMSE evolution of the background (in blue), the observations (in green), the analysis in case a (in red) and the analysis in case b (in black) with respect to the reference state for each velocity components over the assimilation window of T a = 50∆t within the sub-domain Ω obs . . . . . 119 4.12 Assimilation with 3D synthetic observations : snapshots of the velocity field for each component u, v and w at the beginning of the assimilation window at t 1 = 170 in the stream-wise plane z = 3D. The row correspond to the reference, background, 4DVar-all and 4DVar-rampe flow, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 4.13 Assimilation with 3D synthetic observations : snapshots of the velocity field for each component u, v and w at the end of the assimilation window at t 2 = 170.25 in the stream-wise plane z = 3D. The row correspond to the reference, background, 4DVar-all and 4DVar-rampe flow, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 4.14 Assimilation with 3D synthetic observations : comparison of the energy spectrum evolution of the reference flow (in blue), the background flow (in green), the 4DVar-all flow (in red) and the 4DVar-rampe black (in black) at t 1 = 170 (on the left) and at t 2 = 170.25 (on the right) in a cross-section plane x = 4D. The comparison of the spectra of the error for the same data with respect to the reference are shown in inset.

3.3

Comparison of the computation resources required for the classic 4DVar method (4DVar), the parameter perturbation approach (4DEnVar-PP), Gaussian perturbation with a local covariance approach for N = 32 ensemble members (4DVar-LC) and the Gaussian perturbation with a local ensemble approach for N = 32 ensemble members. . . . .

4.4

Characteristics of each run performed to reconstruct the three dimensional flow from 3D synthetic observations. Each data assimilation cases considered a static background error covariance matrix B -1 = 1. The 4DVar and 4DVar-rampe runs fix the observation error covariance matrix to R -1 = 1 and to the formulation given by (4.23), (4.24) and (4.25), respectively. . . . . . . . . . . . . . . . . . . . . . .

4.5

Assimilation with synthetic orthogonal-plane observations : characteristics of each run performed to reconstruct the three dimensional flow from 3D synthetic observations. The empty box case corresponds to the simulation performed to obtain an acceptable initial background from an initial null velocity field and the injection of an inlet every ∆t PIV = 20∆t . The data assimilation method 4DVar considered a static background error covariance matrix B -1 = 1 and the observation error covariance matrix R -1 built by algorithm 29. . . . . . . .

4.6

Assimilation with real orthogonal-plane observations : characteristics of each run performed to reconstruct the three dimensional flow from 3D piv observations. The empty box case corresponds to the simulation performed to obtain an acceptable initial background from an initial null velocity field and the injection of an inlet every ∆t PIV = 20∆t . The data assimilation method 4DVar considered a static background error covariance matrix B -1 = 1 and the observation error covariance matrix R -1 built by algorithm 29. . . . . . . . . . . . . . 133

Introduction

Motivation

The study of turbulent flows plays a prominent role in a wide range of applications. The aeronautics, naval and automobile industries are constantly in need for an accurate representation of the three dimensional wake of an unknown obstacle to build reliable and eco-friendly vehicles that can face the adversity of climate change. At a human scale, the study of turbulent flows can also bring significant improvements to pre-existing systems. For instance, the reconstruction and control of fundamental turbulent flows carried out by the ACTA team (Irstea, Rennes) initiated the development of the FROILOC c system which blows ultra-clean and cold air over a localized work space without altering the ambient air temperature. This system is very attractive for the food and pharmaceutical industries as it ensures the quality of the manufactured products and the well-being of the workers at the same time. All in all, we are most of the time surrounded by turbulent flows which influence human activity at various scales. Therefore, the more we understand the behavior of the turbulence phenomenon in fundamental flows, the more efficiently we can apprehend real life problems and propose efficient technological solutions.

Over the last century, two main fields emerged from the fluid mechanics context to study the behavior of fundamental turbulent flows. In the one hand, the Computational Fluid Dynamics (CFD) field developed a large panel of dynamical models to reconstruct numerically the evolution of steady or unsteady flows from an initial state. On the other hand, the Experimental Fluid Dynamics (EFD) field developed several measurement techniques to observe two dimensional, and more recently three dimensional, observations of each component of the velocity field. During the last decade, the fluid mechanics community witnessed the surge of techniques combining both approaches in order to benefit from the assets and overcome the limitations of each approach [START_REF] Druault | Generation of three-dimensional turbulent inlet conditions for largeeddy simulation[END_REF]Suzuki, 2014).

In recent studies, there is a trend towards the implementation of data assimilation techniques in a fluid mechanics context [START_REF] Cuzol | A low dimensional fluid motion estimator[END_REF][START_REF] Cuzol | A stochastic filtering technique for fluid flows velocity fields tracking[END_REF][START_REF] Kato | Approach for uncertainty of turbulence modeling based on data assimilation technique[END_REF][START_REF] Papadakis | Variational assimilation of fluid motion from image sequences[END_REF][START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF]. Data assimilation emerged from the weather forecasting and the geophysics fields which are characterized by large-sized systems. This thesis follows the work of [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF] and carries out the application of a 4DVar data assimilation technique, to the 3D reconstruction of the characteristics of a given flow from image data and dynamical models. Two different configurations are studies in this thesis.

Free-surface flow configuration

Following the work of [START_REF] Combés | A particle filter to reconstruct a freesurface flow from a depth camera[END_REF] with particle filter, we investigate with a variational data assimilation method the combined use of Kinect depth observations with a shallow water dynamical model to reconstruct free surface flow geometry and dynamics. As mentioned in its name, the Shallow Water Equations (SWE) are used to model a low-depth three dimensional fluid at its surface. The SWE describes the height and the velocity fields at its surface, and are often used as a prototype in geophysics applications. In the framework of this thesis, we built a MATLAB dynamical code from scratch and implemented the 4DVar technique with the help of an automatic differentiation tool. At first, the reconstruction is performed with partial and complete synthetic observations. In a second step, we consider the real measurements of the evolution of the free surface height observed with a Kinect depth-sensor.

In this application, we seized the opportunity to compare the performance of the 4DVar technique with various strategies of an emergent hybrid data assimilation technique. The comparison was published in Computers and Fluids [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF].

3D cylinder wake at Reynolds 300 configuration

This study is an extention of the work of [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF] to the reconstruction of a three dimensional cylinder wake at Reynolds 300 by combining a highly accurate and parallelized code Incompact3d [START_REF] Laizet | High-order compact schemes for incompressible flows : A simple and efficient method with quasi-spectral accuracy[END_REF] with a sequence of two dimensional observations. An important part of this application was focused on the construction and validation of the discrete adjoint parallelized code necessary to the implementation of the 4Dvar method. We performed a first reconstruction of a purely synthetic flow generated by Incompact3d, using three dimensional observations. We then performed the reconstruction of a fully three dimensional flow from the alternated synthetic observations of orthogonal stereo PIV like observations (inflow and stream-wise plane). We investigate the possibilities of the reconstruction with the real observations obtained by orthogonal stereo PIV measurements.

This thesis is articulated around four chapters. At first, we overview in chapter 1 the numerical and experimental approaches developed so far in the fluid mechanics context to reconstruct a fundamental turbulent flow. We present the limitations of each approach and shed light on the attempts to combine both methods in the fluid mechanics context. We overview in chapter 2 the data assimilation techniques which are well suited to our problem and focus in particular on the variational data assimilation methods. We describe in chapter 3 the construction and the implementation of the the 4DVar technique to a free surface flow problem. We compare the reconstruction of the free surface characteristics by different data assimilation techniques [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF] with synthetic and real observations. Finally, we carry out in chapter 4 the implementation of the 4DVar technique to a highly accurate parallel code. The 4DVar code was validated on a purely synthetic experimental case and puts in perspective the reconstruction of a fully three dimensional flow from a sequence of two dimensional observations.

Chapitre 1 Turbulent 3D flow reconstruction

The study of turbulent flows are crucial for a wide variety of academic and industrial applications. While fluid dynamics problems are often associated to aerodynamics and geophysics, they reveal to have a decisive influence in many other fields. From the modeling of the blood flow within a vessel, to the refrigeration process for the quality and safety of food products ; fluid dynamics has played a pivotal role and opened the door to many applications. The study and modeling of fundamental flows is a cornerstone of industrial applications. By reconstructing the evolution of fundamental fluid flows characteristics (velocity and pressure fields, temperature etc...), we can thereafter enhance complex systems used in the industry. Traditionally, there are two main approaches to the 3D reconstruction of turbulent flows that we discuss in the following.

A first approach, described in § 1.1, consists in modeling the behavior of a set of variables that characterize fundamental flows based on physical laws. The Computational Fluid Dynamics (CFD) field offers a wide range of numerical methods to solve the dynamics equations describing the flow evolution. We also introduce the main techniques used in turbulence modeling. Finally, we discuss the inherent limitations of the CFD methods to model a realistic phenomenon, like the choice of the initial and boundary conditions or the tuning of the model's parameters.

A second approach, described in § 1.2, consists in observing directly the evolution of a fundamental flow in a controlled space. In this approach, many efforts have been devoted to measure a fluid motion in a non-intrusive and accurate manner. The experimental approach was formalized in the beginning of the XXth century and grew along with the technical progress of measurement sensors and computer vision techniques. The recent development of Experimental Fluid Dynamics (EFD) offers a bright perspective for the reconstruction of a full 3D turbulent flow within a square meter volume. Despite the remarkable technological progress, the EFD approach can only provide sparse and often noisy information that require an additional model to deduce the quantities of interests, namely the forces exerted on the fluid.

Considering the inherent limitations of the numerical and experimental approaches, the investigation for a combination of both approaches has aroused a lot of enthusiasm in the fluid dynamics community. We discuss the attempts from each community to combine a dynamics model with a set of observations in § 1.3. The most promising approach is given by the data assimilation techniques which stem from the weather forecasting field and are well established in every weather forecasting centers and geophysics applications as well. Data assimilation techniques have captured the attention of many other communities, and recent efforts were notably carried out to implement them in the fluid dynamics context. We overview these endeavors and introduce the approach we followed within this PhD work.

Computational Fluid Dynamics

The simulation of steady and unsteady flow problems rose long before the growth of numerical methods and the advent of computers in the 1970s. Centuries of fluid mechanics research established the governing equations that describe the evolution of an infinitesimal particle in time and space (see § 1.1.1). The latter can be solved exactly by hand for simplified or idealized flows, or approximated by numerical methods. The increasing need to design atmospheric entry vehicles in the 1940s and beyond (for ballistic, then for orbital and lunar motivations), urged for the resolution of complex fluid dynamics problems. In parallel, numerical methods to solve ODEs and PDEs were known and conceived by paper since the time of Newton in the 1700s. Centuries of progress led to a wide range of concepts and numerical methods, however it was a tedious task to apply these methods by hand. Hence, the advent of computers in the 1970s boosted their growth and development ; and while the application of numerical methods go beyond fluid dynamics problem, there was a complete link between the development of numerical methods for PDEs and the simulation of fluid flow problems. It is within this context that the Computational Fluid Dynamics (CFD) field emerged and grew significantly. Nowadays, the literature provides us a wide choice of numerical methods that have been implemented and improved throughout the history of CFD. We briefly overview these methods in § 1.1.3. The development and improvement of the CFD field allowed us to study more and more complex flows, in particular the reconstruction of turbulent flows which is the key to many industrial applications. We briefly discuss the case of turbulence flow modeling principles in § 1.1.2. Finally, we discuss in § 1.1.4 the inherent limitations of the CFD approach.

Governing equations

We generally describe a fluid motion through three conservation laws that link the mass density ρ, the pressure field p = p(x, y, z), the velocity field in each direction u = (u, v, w) and the energy E. The first conservation law guarantees the conservation of mass, the sum of the mass flow rates entering into a control volume is equal to the sum of the mass flow rates getting out of the control volume. The second conservation law comes from the application of Newton's second law of motion to a given fluid particle, the change of momentum ρu = (ρu, ρv, ρw) equals the sum of all the external forces applied to each fluid particle. Finally, the third conservation law, which comes from the first law of thermodynamics, guarantees the conservation of energy. The rate of change of energy equals the sum of the rate of heat and the work done on each fluid particle.

Within the fluid dynamics context, we generally consider the Navier-Stokes equations which describe the relation between the velocity and the pressure of a moving fluid. In an academic context, we generally assume that the variations of pressure doesn't influence the particles' density. This ideal case is specified by the incompressibility hypothesis, which states that ∇•u = 0. In the Newtonian fluids approximation, we assume that the relation between stress and the corresponding rate of strain is linear. The viscous stress in the fluid is the sum of a diffusion viscous term (proportional to the gradient of velocity) and a pressure term. Turbulence, and the generation of boundary layers, are the result of diffusion in the flow. Under the previous hypotheses, the momentum equation is rewritten into the Navier-Stokes equations for an incompressible fluid. Hence, for an incompressible flow, the fluid dynamics equation reads,

∇•u = 0,
(1.1)

∂ t u + 1 2 (∇(u ⊕ u) + (u∇)u) -ν∇ 2 u + ∇p = 0. (1.2)
The Navier-Stokes equations, coupled with the adequate initial and boundary conditions, are primarily used in fluid dynamics to describe fundamental flows and will be considered within this thesis in chapter 4.

Most of the dynamics fluid equations encountered in many other configurations are generally derived from the latter. We will also consider in this work the shallow water equations (SWE) which are deduced from the depth-integration of the Navier-Stokes equations. These equations are often used in an academic context to model the free surface of a fluid which depth is negligible with respect to its surface dimensions. While a vertical velocity term is not present in the shallow water equations, note that this velocity is not necessarily zero and can be recovered via the continuity equation. The standard 2D SWE are formulated as :

   ∂ t h + ∂ x (hu) + ∂ y (hv) = S h , ∂ t (hu) + ∂ x (hu 2 ) + ∂ y (huv) + g∂ x h 2 = S hu , ∂ t (hv) + ∂ x (huv) + ∂ y (hv 2 ) + g∂ y h 2 = S hv ,
where g is the gravity acceleration, S h , S hu and S hv are the source terms such as, for instance, the Coriolis force S hu = f hu, the friction S hu = c f u √ u 2 + v 2 or the topology at the "bottom" of the fluid.

Introduction to turbulence modeling

While the formal definition of turbulence is still in the center of numerous debates, there is a common ground on the definition of its main characteristics. In the turbulence context, the smallest scales, referred to as the Kolmogorov scales η, have a crucial impact on the characteristic length of the smallest scales of the flow beyond which kinetic energy is transformed into heat (sometime called eddies). Therefore, an ideal simulation of a turbulent flow would consist in representing the contributions of all the spatial and temporal scales of the turbulence spectrum (see figure 1.1). This complete representation can be achieved by the Direct Numerical Simulation (DNS). Once we define a spatial and temporal discretization of the physical domain that can represent the smallest scales, we solve the Navier Stokes equations by the means of one or a combination of numerical methods presented in § 1.1.3. As one could expect, the DNS requires a tremendous amount of computational resources which increase with the size of the problem. We can assess the computational cost by considering a flow characterized by a given Reynolds number Re, the number of points necessary to describe all the turbulent flow scales depends on Re such as, N p ∝ Re 9/4 .

In addition, the number N t of time steps ∆t necessary for the simulation on a time interval [0; T ] is given by

N t ≈ T ∆t ≈ Re 1/2 .
Therefore, the cost of the DNS on a [0; T ] time interval is estimated to

N p × N t ∝ Re 3 .
The higher the Reynolds number, the more expensive the DNS. In general, engineering applications require high Reynolds number, like for instance Re = 10 8 and Re = 10 9 in automobile and aerodynamic applications, respectively. While the DNS is clearly inapplicable in these applications, it has greatly contributed to the development of turbulence models that we introduce in this section. We outline three notable turbulence models and refer the reader to [START_REF] Wilcox | Turbulence modeling for CFD[END_REF]Spalart (2000); [START_REF] Lesieur | Turbulence in fluids : stochastic and numerical modelling[END_REF] for a deeper description of these models.

Reynolds Average Navier Stokes (RANS) consists in dividing the velocities into mean components averaged in time and fluctuating components, by introducing them into the Navier Stokes equations, then modeling the new product terms. The RANS model has been implemented among commercial softwares such as ANSYS (former FLUENT). We refer the reader to Spalart (2000) for a review of RANS models.

Large Eddy Simulation (LES) consists in modeling the small eddies by a subgrid scale model, pioneered by Smagorinsky (1963) and improved by [START_REF] Germano | A dynamic subgridâĂŘscale eddy viscosity model[END_REF] for near-wall region, while the large eddies are solved directly by a transient calculation. The reader can find in [START_REF] Sagaut | Large eddy simulation for incompressible flows : an introduction[END_REF] a broader review of LES methods.

Mathematically, the LES comes to separating the velocity field into a resolved part, representing the large eddies, and a sub-grid part, representing the small scales.

More recently, [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] has proposed an alternative to the classical Reynolds decomposition by separating a differential drift component to a time uncorrelated uncertainty random term. The latter isn't differentiable with respect to time, therefore it is handled through stochastic calculus. The dynamics associated with the differentiable drift component is derived from a stochastic version of the Reynolds transport theorem. It includes in its general form an uncertainty dependent subgrid bulk formula that cannot be immediately related to the usual Boussinesq eddy viscosity assumption constructed from thermal molecular agitation analogy like in LES subgrid models. This formulation, considering uncertainties into the dynamical turbulence flow model is promising for the coupling with sparse data image. The learning of spatio-temporal model parameters taking into account data uncertainties and model error is potentially of huge interest, as for the moment no satisfying solution exists in complex real world cases encountered in industrial problems.

Numerical methods for CFD

Exact solutions of the Navier Stokes equations (1.1), (1.2) can be found in [START_REF] Crane | Flow past a stretching plate[END_REF]; [START_REF] Brady | Steady flow in a channel or tube with an accelerating surface velocity. an exact solution to the navier-stokes equations with reverse flow[END_REF]; [START_REF] Wang | Exact solutions of the steady-state navier-stokes equations[END_REF]. These solutions are often restricted to a steady-state flow and boundary conditions, and they eventually require additional simplifications of the equations. On the other hand, numerical solvers of ODEs and PDEs were theoretically established in the literature. After over a half-century of research, we have at our disposal a wide range of numerical methods and discretization schemes to solve a set of PDEs at a certain level of accuracy. In the CFD framework, we commonly use the Finite Difference Method (FDM), the Finite Element Method (FEM), the Finite Volume Method (FVM) and the Spectral Element Method (SEM).

Finite Difference Method (FDM) is the first numerical method to be applied to fluid dynamics problems. In the same way as the definition of a derivative, the FDM replaces the partial derivatives of the governing equations by algebraic difference quotients. This yields to a system of algebraic equations from which we deduce the value of the variables of interest at each grid points. The two sources of error in the FDM are round-off errors, the loss of precision due to computer rounding of decimal quantities, and truncation error or discretization error, which is the difference between the exact solution of the finite difference equation and the exact quantity assuming no round-off errors. The easy implementation of the FDM is an important asset which allowed its fast improvement and spread to various domains of application. In the context of turbulent flows, it is necessary to deal with highly accurate schemes which are capable to represent all the ranges of scales while ensuring a stability. The compact finite difference schemes [START_REF] S K Lele | Compact finite difference schemes with spectral-like resolution[END_REF] keeps the versatility of the FDM and provides a spectral-like accuracy. These methods are ideal for DNS as they provide a sufficiently high accuracy and robustness.

Finite element method (FEM) stems from structural mechanics problems which focus on the deformation, deflections and internal stresses within structures. These problems are characterized by complex geometry that can't be represented by a regular structured grid. Therefore, the FEM offers the flexibility to subdivide the domain into cells, called elements, of flexible shape. The elements have typically a triangular form (in 2D) but they can also be rectangular or curved which is an important asset with respect to FDM. The reader can find more description in Szabo and Babuvska (1991). In opposition to the FDM, the FEM doesn't seek for the solution of the PDE itself, but aims at solving the integral form of the PDE instead. This asset leads to more flexibility in terms of boundary conditions implementations and allows the construction of higher order accurate methods. However, the implementation of the FEM itself is a tedious task and is actually worth the effort for a solution that is assumed to depend a priori on a complex geometry. In this perspective, FEM is well established in industrial applications and is commonly implemented entirely or partially in most commercial tools such as ANSYS (formely Fluent) or FLOW-3D.

Finite Volume Method (FVM) is well suited to simulate various types (parabolic, elliptic, hyperbolic etc...) of conservation laws and is usually favored to solve CFD problems in particular [START_REF] Versteeg | An introduction to computational fluid dynamics : the finite volume method[END_REF]. In opposition to the FDM, the FVM is based on the variational formulation of the PDE. Instead of determining the value at a local grid point, the FVM quantifies the fluxes that move from a discretized cell -referred to as a control volume -to its neighbor. By construction, the FVM ensures the numerical conservation of the flux which is an attractive asset to the fluid dynamics context. Similarly to the FEM, the FVM may be used on complex geometries, using structured or unstructured meshes but requires less implementation effort.

Spectral Method (SM) are established as the most accurate numerical method to solve fluid dynamics PDEs. Spectral methods [START_REF] Gottlieb | Numerical Analysis of Spectral Methods : Theory and Applications[END_REF] consists in using Fourier or Chebyshev representation and solving the governing equations in the spectral space which guarantees the high accuracy of the method. However, these methods are usually limited to simple geometries and for complex configurations, conventional finite difference methods are used. Nevertheless, the accuracy of the SM is a very attractive asset which pushed the development of hybrid methods. For instance, the Spectral Element Method, introduced by [START_REF] Patera | A spectral element method for fluid dynamics : Laminar flow in a channel expansion[END_REF], decomposes the domain into a series of elements similarly to the FEM and the velocity in each element is represented as a high-order Lagrangian interpolant through Chebysev collocation points. Similarly, the Spectral Volume Method [START_REF] Wang | Spectral (finite) volume method for conservation lazs on unstructured grids. basic formulation[END_REF][START_REF] Wang | Extension of the spectral volume method to high-order boundary representation[END_REF] decomposes each spectral volume into a series of control volumes. The cell-averaged data from these control volumes are used to reconstruct a highorder approximation in the spectral volume. The Riemann solvers are then used to compute the fluxes at spectral volume boundaries.

In the framework of this PhD, our goal is to perform a DNS model for two fluid flow configurations using the appropriate numerical schemes. In a first academic application, we consider a simplified 2D Shallow Water model defined on a simple geometry form (see chapter 3). We implemented the FVM, as described in detail in § 3.1, which provides an accurate solution for a reasonable implementation effort. On the other hand, we perform a DNS for a 3D turbulent cylinder wake for Re = 300 based on the incompressible Navier Stokes equations. The DNS code employed, Incompact3d, combines a spectral method and high order compact FDM as described in [START_REF] Laizet | High-order compact schemes for incompressible flows : A simple and efficient method with quasi-spectral accuracy[END_REF].

Limitations of computational fluid dynamics

As depicted throughout this section, the final goal of CFD methods is to provide a numerical approximation of the actual physical laws that governs a realistic flow. While the CFD field provides accurate methods and tools to represent ideal flows, there are still inherent limitations that may lead to a poor representation of a turbulent realistic flow in three dimension. We discuss in this section the main limitations we can encounter in a purely CFD approach.

Model errors are inherent to the CFD approach. Indeed, we don't always have enough of knowledge of the physical laws governing the turbulent flow of interest itself. Typically, in many fundamental and industrial problems, it is nearly impossible to know the exact boundary conditions. Despite we achieve a perfect representation of the flow within the volume of interest, a bad choice of the boundary conditions model might lead to the spread of errors within the whole volume and compromise the final solution of the problem. In a channel flow configuration for instance, it is nearly impossible to model a realistic inflow for a turbulent flow which is crucial to define the evolution of the flow throughout the channel. In a three dimensional cylinder wake configuration in particular, it is custom to impose a periodicity condition along the cylinder axis. Indeed, this choice is based on experimental measurements of [START_REF] Zdravkovich | Flow around Circular Cylinders[END_REF], which shows a periodic distribution (2 < λ z /D < 9) depending on the Reynolds number.

Additional errors also rise when we consider turbulence models such as (LES, RANS). These errors are induced by the model itself and the adjustment of the parameters of the model. These errors can be eventually included in the dynamical model and the parameters can be estimated in situ (see § 1.1.2).

The initial condition definition is also a major limitation to the CFD approach. Even if we possess a perfect dynamics model and we input a misinformed initial condition, we might obtain a completely different flow. This comes from the nonlinear terms of the Navier Stokes equations and the deterministic nature of the CFD model. In addition, as we mentioned previously, a turbulent flow is sensible to the finest scales of the flow. Therefore, in a turbulent flow with a high Reynolds number, a slight change in the finest scale structures can impact the large scale structures.

Numerical errors induced by the numerical schemes employed to discretize the dynamics model ( § 1.1.3) can also compromise the representation of realistic turbulent flows. Despite the remarkable progress in the CFD field and the constant growth of computational resources, we are inevitably exposed to round-up and truncation errors. Such errors can be lowered by the representation of the flow on a finer grid, however the latter requires a considerable increase of computation time and storage.

In a nutshell, CFD models provide an essential pillar of the reconstruction of turbulent flows as they provide time continuity and dynamical consistency among the velocity and pressure fields. Nevertheless, these models are only able to reconstruct accurately ideal flows and we still lack knowledge on the physical model of realistic three dimensional turbulent flows. This lack mainly concerns the initial and boundary conditions, and the turbulence model errors and the tuning of their parameters. Additionally, we have to find the best compromise between the computational resources and the accuracy of the numerical schemes used to discretize the dynamics model.

Experimental fluid dynamics

One of the early and most prominent example of fluid mechanics experiment has been carried out by L. Prandtl in 1904 who observed and studied the aspects of unsteady separated flows behind several objects in a water tunnel(figure 1.2). The flow was visualized by distributing a suspension of mica particles on the surface of water. The latter set up a revolutionary experimental methodology that has systematically been used until today, and paved the way for measurement techniques expansion. Likewise the CFD field, the Experimental Fluid Dynamics (EFD) field expansion has been closely linked with technological progress. Within the fluid dynamics context, the intrusive Hot Wire Anemometry (HWA) and the non intrusive Laser Doppler Anemometry (LDA) were popularly used and developed to measure an accurate and punctual point in space. The estimation of the whole velocity field could be obtained by an additional post-processing technique. While these techniques are well established in the community, they were progressively supplanted by Particle Image Velocimetry (PIV) and later on 3D Particle Image Velocimetry which are briefly discussed throughout this section.

Particle Image Velocimetry for 2D observations

Nowadays, by combining Prandtl's attempt to trace particles, a contemporary laser and image processing techniques, we are able to measure the velocity of micronsized particles following the flow by the so-called Particle Image Velocimetry (PIV) technique, illustrated in figure 1.3. PIV is a well established non-intrusive 2D measurement technique for measuring velocity fields of fundamental flows from a sequence of instantaneous snapshots [START_REF] Bibliographie | Particle-imaging techniques for experimental fluid mechanics[END_REF]. The displacement of the particles is estimated statistically by correlating two sequences of particle images as described in [START_REF] Meinhart | A parallel digital processor system for particle image velocimetry[END_REF]. If we consider for instance an airflow at Re = 300, i.e. ) ≈ 72, the PIV experiments is able to provide the same spatio-temporal resolution than a DNS. However, if we consider a higher Reynolds number, the corresponding dynamics range is L η >> 100 whereas the PIV is limited to a velocity dynamics range of range of 100 :1 [START_REF]Limiting resolution of particle image velocimetry for turbulent flow[END_REF]. Therefore, at the PIV resolution we can't retrieve the smallest scales yielding a coarse representation of a turbulent flow. The reader can find a complete description of the PIV technical background and applications in [START_REF]Twenty years of particle image velocimetry[END_REF] and in [START_REF] Raffel | Particle image velocimetry : a practical guide[END_REF].

Figure 1.3 -The classical PIV system configuration setup is composed of a doublepulsed laser and CCD camera. We deduce the velocity field by performing the analysis of the correlation between a sequence of frames.

The measurement of a two-dimensional three-component velocity field is a major stake for fluid mechanics and is even indispensable for a good understanding of turbulent flows. By combining PIV and stereoscopy, [START_REF] Arroyo | Stereoscopic particle image velocimetry[END_REF] developed a set-up which takes two stereoscopic images of the flow simultaneously with only one camera, yielding a measurement of the three components velocity flow in one plane. [START_REF] Kähler | Fundamentals of multiple plane stereo particle image velocimetry[END_REF] developed a multi-plane stereo PIV system which consists of a four-pulse laser system and two pairs of progressivescan cameras in an angular imaging configuration. This system provides the three components velocity field based on a classical PIV setting.

3D Particle Image Velocimetry

The progress of PIV techniques led to satisfying three components reconstruction of the velocity fields. Nevertheless, in order to have a better understanding of turbulent flows, the three dimensional three component velocity field is of crucial need. Several techniques to measure three dimensional quantities were proposed in the literature, such as scanning-PIV [START_REF] Brücker | 3d scanning piv applied to an air flow in a motored engine using digital high-speed video[END_REF], 3D Particle Tracking Velocimentry (3D-PTV) [START_REF] Virant | 3d ptv and its application on lagrangian motion[END_REF] or holographic-PIV [START_REF] Hinsch | Holographic particle image velocimetry[END_REF][START_REF] V S S Chan | Application of holography to fluid flow measurements using bacteriorhodopsin (br)[END_REF]. Tomographic PIV (tomoPIV), pioneered by [START_REF] G E Elsinga | Tomographic particle image velocimetry[END_REF], offers one of the most interesting perspective for three dimensional unsteady flow measurements. In analogy to the PIV technique, tomoPIV consists in seeding the measured 3D volume with tracer particles and illuminate the 3D region by a pulsed light source. The particles are then recorded simultaneously by CCD cameras similar to stereo-PIV as illustrated by figure 1.4. The tomoPIV measurement technique has recently emerged thus we are able to retrieve very small volumes using thin tracer particles. The first tomo-PIV measurements carried out by [START_REF] G E Elsinga | Tomographic particle image velocimetry[END_REF] retrieved the three components of the a cylinder wake velocity field within a 40 mm × 40 mm × 10 mm volume in air using 1µm droplets as tracer particles. Most of the tomo-PIV experiments achieved so far remain in that scale order (Schröder et al., 2008;[START_REF] Humble | Three-dimensional instationary structure of a shock wave turbulent boundary layer interaction[END_REF] in order to provide an accurate representation of the particle motion.

Figure 1.4 -The tomographic PIV system configuration setup is composed of a double-pulsed laser and CCD camera. We deduce the velocity field by performing the analysis of the correlation between a sequence of frames.

More recently, there is an emergence of large-scale tomo-PIV systems to retrieve a cubic-meter sized volume by using Helium-Filled Soap Bubbles (HFSB) [START_REF] Kuhn | Large-scale tomographic particle image velocimetry using helium-filled soap bubbles[END_REF][START_REF] Scarano | On the use of helium-filled soap bubbles for large-scale tomographic piv in wind tunnel experiments[END_REF]. In this configuration, we seed the turbulent flow with HFSB of which are much larger than the 1µm particle used so far in PIV and tomoPIV techniques. Due to their size, the particle distribution is smaller and more sparse thus they are esentially able to measure large scale structures. Other measurement techniques also such as tomoPTV and 4D-TomoPIV "Shake the box" also propose to retrieve the three dimension three component velocity fields in for wide volumes [START_REF] Schanz | shake the box' : A highly efficient and accurate tomographic particle tracking velocimetry (tomo-ptv) method using prediction of particle positions[END_REF].

Limitation of Experimental Fluid Dynamics

The EFD techniques presented throughout this section are able to capture phenomena that we might not be able to model yet and that improves our knowledge on the flow of interest. Despite the technological and methodological progress in the measurement techniques, there are still inherent limitations that may restrict severely the representation of the realistic turbulent flows. We discuss in this section the limitations that we can encounter in the 2D and 3D measurement techniques.

2D measurements are first and foremost limited in size, PIV measures are typically limited to 1000 × 1000px 2 which comes to 8 × 8cm 2 . In terms of spatial resolution, as mentioned previously the PIV is limited to a velocity dynamics range of 100 : 1 [START_REF]Limiting resolution of particle image velocimetry for turbulent flow[END_REF]. Therefore, at a high Reynolds number, we are only able to capture large scales and obtain the same resolution as a LES (see § 1.1.2). In addition to the resolution limitation, 2D measurements are inherently sparse in space, as they only provide the velocity components within planes, and in time. We can eventually set an experimental configuration to retrieve the velocity fields at the inlet and the stream-wise planes alternatively in time (see § 4.4.3), however we still lack the whole 3D contribution of the measured flow.

3D measurements are even more limited in size due to the tremendous cost of energy of the laser used to illuminate the flow. As described in § 1.2.2, tomoPIV measurements are limited to a small volume of a few cm 3 and the velocity dynamics range is smaller than the 2D measurement case. Therefore, similarly to the 2D measurements, we can't retrieve the smaller scales and obtain observations which are LES like. The recent large-scale tomo-PIV techniques can eventually retrieve a larger domain of a few m 3 [START_REF] Kuhn | Large-scale tomographic particle image velocimetry using helium-filled soap bubbles[END_REF][START_REF] Scarano | On the use of helium-filled soap bubbles for large-scale tomographic piv in wind tunnel experiments[END_REF], however the obtained observations are sparse in order to avoid occlusions of the particles.

In a nutshell, EFD techniques provide important characteristics of a realistic flow. However, the measurements size and resolution are limited by the experimental settings and are often sparse in time and in space. Therefore, we obtain a LES-like observations which captures the large scales of the flow of interest. Let us note this LES-like representation is not model-free since strong smoothing functions are used in the estimation process.

Coupling CFD with EFD

Traditionally, experimental observations are used as a reference for a numerical method result as in [START_REF] Dong | A combined direct numerical simulation-particle image velocimetry study of the turbulent near wake[END_REF] and in [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at reynolds number 3900[END_REF], while the model can be used to deduce non-observable variables characterizing the fluid dynamics. Recently, the CFD and the EFD communities have both shown great interest in combining a dynamical model with observations to obtain a better representation of a given flow motion. We present an overview of the early attempts to introduce measurements into a model and vice versa. A particular focus is brought to the promising data assimilation techniques that couple a model with the observations of a system of interest. We describe the recent attempts to apply a data assimilation technique in the fluid mechanics context and finally present the approach proposed by this thesis.

Early attempts to combine measurements with a model

Early attempts to introduce measurements into a model in the fluid mechanics community stem from the application of the Proper Orthogonal Decomposition (POD). Initially, POD was used in the CFD community to reduce large CFD models. The popular snapshot method proposed by Sirovich (1987) provides a set of modes from a set of instantaneous flow solution obtained from a simulation of the CFD method. The governing equations of the large CFD models are then projected onto the reduced space spanned by the POD modes. The snapshot method was thereafter improved by [START_REF] Everson | Karhunen-loeve procedure for gappy data[END_REF] in the image processing context, to reconstruct of images -human face -from a sequence of incomplete or gappy images. The first application of the POD in the aerodynamics context was carried out by [START_REF] Bui-Thanh | Aerodynamic data reconstruction and inverse design using proper orthogonal decomposition[END_REF] to reconstruct flow fields corresponding to a given airfoil shape, using a set of incomplete experimental PIV measurements. In the fluid mechanics context, [START_REF] Druault | Generation of three-dimensional turbulent inlet conditions for largeeddy simulation[END_REF] use a POD model which combines a LES model with HWA and PIV measurements to generate realistic inflow conditions for a numerical simulation. This procedure can also be used as a dynamical data post-processing which aims at improving the experimental data, for instance by completing partial information or restoring the data deterioration occurred during the measurement acquisition. More recently, Suzuki (2014) applied a POD model by combining with time-resolved PIV/PTV measurements for a turbulent flow at Reynolds Re = 2000.

Optimal design problems aim at determining the optimal input parameters corresponding to measured forces exerted on the fluid. In this context, we seek for an optimal design parameter, for instance an initial condition [START_REF] Lundvall | Iterative methods for data assimilation for burger's equation[END_REF] or a material characteristic [START_REF] Avril | Further observations on the mean velocity distribution in fully developed pipe flow[END_REF], such that the error between the model-based solution and the true measurements reaches a minimum. To that purpose, optimal design problems are formulated as an optimal control problem [START_REF] Lions | Optimal control of systems governed by PDEs[END_REF] which is solved by an adjoint-based control optimization, as pioneered by [START_REF] Jameson | Optimum aerodynamic design using cfd and control theory[END_REF]. In this aspect, optimal design problems are similar to the variational approach of the data assimilation techniques that we present in the following sub-section.

Towards data assimilation techniques

Data assimilation (DA) is a whole branch of mathematical techniques that emerged from the weather forecasting and geophysics fields. One of the most seducing assets of the data assimilation techniques is their ability to retrieve the state of interest in a full spatial and temporal resolution by combining the dynamics model describing the system of interest with low resolution and partial observations. The application of DA techniques to the fluid mechanics context was first investigated in D 'Adamo et al. (2007) and in [START_REF] Cuzol | A low dimensional fluid motion estimator[END_REF] to provide low representations of the flow characteristics. We usually distinguish two main approaches to data assimilation techniques, a sequential data assimilation (SDA) approach which estimates an optimal state from a first guess state and its observation ; and a variational data assimilation (VDA) approach which seeks for an optimal trajectory which is a compromise between the evolution of a system and a set of observations in a given assimilation window.

Application of the SDA approach In the SDA approach, Colburn et al. ( 2011) employed the Ensemble Kalman Filter (EnKF) technique, introduced by [START_REF] Evensen | The ensemble kalman filter : Theoretical formulation and practical implementation[END_REF], to estimate the 3D near-wall flow at a low Reynolds number (Re=100). This application combines a standard spectral second order finite difference approach with measurements of the skin friction and pressure only on the wall. The results are a first step to determine the dynamics law of turbulent near-wall flows with respect to the distance from the wall. [START_REF] Kato | Approach for uncertainty of turbulence modeling based on data assimilation technique[END_REF] employed the EnKF technique to retrieve the optimal parameters of the Spalart-Allmaras turbulence model. These parameters were initialized by random values within a defined range and were compared to the original values proposed by the Spalart-Allmaras turbulence model. [START_REF] Combés | A particle filter to reconstruct a freesurface flow from a depth camera[END_REF] proposed to reconstruct the 3D geometry dynamics of a free surface flow by coupling Kinect depth measurements with SWE dynamics model via stochastic filtering.

Application of the VDA approach In the VDA approach, a first study was carried out by [START_REF] Papadakis | Variational assimilation of fluid motion from image sequences[END_REF] to couple turbulent flow image sequences with DNS. This study was followed by [START_REF] Cuzol | A stochastic filtering technique for fluid flows velocity fields tracking[END_REF] who coupled turbulent flow image sequences with vortex particle simulations. Suzuki et al. (2009a,b) propose a hybrid unsteady-flow simulation technique which combines DNS and PTV which is very similar to the variational approach formulation. The latter is used to and is compared to a 2D URANS solution. These first attempts are reviewed in [START_REF] Heitz | Variational fluid flow measurements from image sequences : synopsis and perspectives[END_REF]. The VDA approach has later on been used by [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF] to reconstruct a 2D accurate full resolution inflow and initial condition by coupling PIV measurements and a highly accurate DNS model at Reynolds Re = 400. This reconstruction was performed by controlling the initial and inlet conditions. In the same vein, [START_REF] Dovetta | Assimilation de données pour la reconstruction d'écoulements et l'identification de turbulence[END_REF] recovered a turbulent pipe flow statistics at high Reynolds (Re = 37155) by performing the VDA approach with a RANS model and mean velocity measurement taken from [START_REF] Mckeon | Further observations on the mean velocity distribution in fully developed pipe flow[END_REF]. In the same vein of [START_REF] Fujisawa | Evaluation of pressure field and fluid forces on a circular cylinder with and without rotational oscillation using velocity data from piv measurement[END_REF] and Van Oudheusden (2005), a variational assimilation technique was carried by [START_REF] Lemke | Adjoint based pressure determination from piv-data -validation with synthetic piv measurements[END_REF] to reconstruct pressure fields for acoustic fluctuations. The latter combines an adjoint-based formulation of the fully 2D compressible Navier Stokes equation, with synthetic PIV observations to determine an improved full resolution of the velocity fields. The pressure field is thereafter deduced by solving the Poisson equation with the obtained velocity fields. The resulting pressure fields revealed to be more accurate than those obtained directly from solving the Poisson equation on the observed PIV field. [START_REF] Foures | A data-assimilation method for reynolds-averaged navier-stokes-driven mean flow reconstruction[END_REF] combined a 2D RANS model with several measures of the mean flow for a low Reynolds number (Re = 100). These studies laid out the VDA asset to interpolate, reconstruct and extrapolate the flow statistics from partial observations, in regions where measurements are difficult or impossible to obtain. The VDA approach was also used in [START_REF] Mons | Is isotropic turbulence decay governed by asymptotic behavior of large scales ? an eddy-damped quasi-normal markovian-based data assimilation study[END_REF] to reconstruct the initial kinetic energy spectrum by coupling an eddy-damped quasi-normal Markovian model with the observation of ah homogeneous isotropic turbulence decay.

Hybrid approach More recently, there is an emergence of hybrid methods which combine the SDA and the VDA approaches which is introduced in § 2.5. A complete review of such method can be found in [START_REF] Yang | Study of variational ensemble methods for image assimilation[END_REF]. In chapter 3, we considered the application of the hybrid technique En4DVar and compare its results with a classic VDA approach.

In the framework of this thesis, we carried out two applications of the VDA approach. In the one hand, we reconstruct a free surface height and velocity field by coupling a 2D SWE model and depth map observations provided by the Kinect sensor. On the other hand, we combine a DNS model with stereo-PIV observations to reconstruct the three components three dimension velocity fields of a cylinder wake at Re = 300. In extension to [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF], we perform the reconstruction by controlling the initial and inflow conditions. One of the biggest challenge of this second application is to retrieve the whole volume in the wake of an unknown obstacle from 2D observations. Chapitre 2

Data assimilation techniques

The concept of data assimilation (DA) comes from the numerical weather prediction (NWP) field. The goal is to combine an a priori knowledge of a given system with observations of the system as formulated in section § 2.1. When we apply data assimilation techniques to a given dynamic system (i.e. we consider the evolution in time), we ought to fulfill three main purposes : improve our knowledge of the present state from the past observations, predict a better state in the future from the past and present observations and provide a better estimate of the evolution of the state, referred to as the trajectory, of the system of interest. Early attempts to combine a model with a set of observations in the NWP context were carried out through the interpolation of the observations to the dynamic numerical grid points. Pioneering works by [START_REF] Bergthorsson | Numerical weather map analysis[END_REF] and [START_REF] J P Cressman | An operational objective analysis system[END_REF], later combined with a statistical interpretation by [START_REF] Gandin | Objective analysis of meteorological fields : Gidrometeorologicheskoe izdatel'stvo (gimiz), leningrad[END_REF], together gave birth to a methodology called optimal interpolation (OI). The estimation theory offers a wide range of parametric estimators for a state of interest given a set of observations, which is discussed in section § 2.2. Thereafter, two main approaches to data assimilation techniques stood out and were developed mostly independently from one another. In the one hand, the sequential approach, discussed in § 2.3, resorts to the Kalman Filter equations which require the calculation of a single and expensive matrix (called the Kalman gain), and on the other hand, the variational approach, discussed in section § 2.4, is based on the application of optimal control theory elements to the data assimilation problem. Finally, we discuss in § 2.5 the emergence of the hybrid data assimilation techniques which combines the concepts of the sequential and varitional approaches.

General context

First and foremost we define the mathematical context common to the data assimilation techniques described throughout this chapter. The common goal of data assimilation techniques is to provide an optimal estimate of a true initial state x t 0 and control parameter u t (t, x) that are in practice unknown. The obtained optimal state is referred to as the analysis state, denoted x a 0 . Data assimilation techniques require two main ingredients. In the one hand, we need a dynamic model that describes the evolution of the system of interest, named the background trajectory, from an a priori initial state x b 0 . On the other hand, we need a set of measurements Y of the system. In the fluid mechanics context, the state variables of interest are typically the velocity and pressure fields and the observations are obtained by PIV or tomoPIV measurements. The data assimilation problem is formulated by the following system of equations, ∂ t x(t, x) + M(x(t, x)) = q(t, x),

(2.1)

x(t 0 , x) = x b 0 (x) + η(x), (2.2) Y(t, x) = H(x(t, x)) + ε(t, x).
(2.

3)

The first equation (2.1) describes through a differential operator M the evolution of the state x in time and space t, x. The differential operator M is often nonlinear ; in the fluid mechanics context for instance, (2.1) corresponds to the general form of the Navier Stokes equations (1.1),(1.2). Function q(x, t) takes into account potential model errors, such as dynamical parameter errors or errors induced by the discretization of the continuous dynamics equations into the numerical solver. In the framework of this thesis, we assume that the dynamics model is perfect, thus q(t, x) = 0.

The second equation (2.2) describes the initial condition for the dynamics equation which is defined by the background state x b 0 up to the background error η(x). The latter is supposed to have a zero mean and is associated to a so-called background covariance matrix,

B = E((x b 0 -x)(x b 0 -x) T ) = E(ηη T
). The last equation (2.3) links the measurements Y of the system with the state variable x for each occurrence of the observation t. This relation is characterized by a differential observation operator H which can be non-linear. The measurements Y are often noisy and can be completely or partially defined in time and space at a different scale than the state of interest x. The discrepancy error between the measurement and the state in the observation space, denoted ε(t, x) = Y(t, x) -H(x(t, x)) and is assumed to be a zero mean Gaussian random field associated to the observation covariance tensor R. The observation error ε(t, x) covers at the same time the measurement errors induced by the measurement sensors and by the observation operator H. In the former case, we usually have an information on the mean and standard deviation corresponding to the adopted measurement sensor. Whereas, in the latter case, it is not trivial to estimate the errors induced directly by the observation operator which depends on a given model. In both applications studied within this thesis, we will consider pseudo observations which are equivalent to the state of interest and which are obtained from a specified pre-process, therefore the observation operator will be set to an identity matrix.

Linear estimation theory

Estimation theory is a branch of statistics which provides a basis for data assimilation techniques. We consider a noisy observation Y ∈ R m of a given physical system, linked to a state variable x ∈ R n by an observation operator. We will consider throughout this section a linear observation operator noted H such as

Y = Hx + ε.
The discrepancy error ε between the observations and the state variable is assumed to be unbiased and is associated to the following covariance matrix

R = E((Y -Hx)(Y -Hx) T ) = E(εε T ).
Estimation theory aims at determining and studying the properties of an optimal estimator x of the true state x t , based on the knowledge of the observation Y. In the Numerical Weather Prediction (NWP) context, we have an a priori knowledge of the data behavior, thus, we are interested into parametric estimators which are very dependent on the adjustment goodness of the data to the density functions. A wide range of estimators can be found in the literature, we present two main approaches to solve the given linear estimation problem.

Linear least squares estimation

Given the estimation problem above, the linear Least Squares (LS) estimation formulation is straightforward. In this approach, we seek for the optimal estimator x that minimizes the following cost function

J(x) = 1 2 (Y -Hx) T (Y -Hx).
The optimal state is obtained by canceling the cost function gradient :

∇J = H T Hx -H T Y = 0,
supposing H T H is reversible, we deduce the optimal linear LS estimator :

x = (H T H) -1 H T Y.
This formulation can be generalized by taking into consideration a weight, described by the observation covariance matrix R, into the cost function :

J(x) = 1 2 (Y -Hx) T R -1 (Y -Hx).
We deduce, the generalized optimal linear LS estimator is given by

x = (H T R -1 H) -1 H T R -1 Y.
We note that the LS approach is purely deterministic technique, which relies on a Gaussian a priori of the discrepancy errors between the data and the state variables. Such a statistic analysis can be pushed forward within a Bayesian context.

Bayesian approach

In the Bayesian estimation, we consider the state vector x as a random variable associated to the probability distribution P (x). We assume that the observation Y is also a random vector. In this approach, we seek for an estimation of the true (unobservable) state x t knowing the observed data Y associated to the conditional probability P (x|Y). By applying the Bayes' rule, the latter is written as

P (x|Y) = P (Y|x)P (x) P (Y) , (2.4)
where P (Y|x) is the distribution of the data given the un-observable state x t . As the observation data come from imperfect measurements, P (Y|x) quantifies the distribution of measurement errors, it can reflect possible biases as well as instrument errors. P (x) is the prior distribution which provides an a priori understanding of the quantities of interest. Within this PhD framework, the prior distribution will be based on the dynamics model defined by physics laws. The marginal distribution P (Y) is in general not needed as it defines a normalizing constant that can be replaced by the following definition,

P (Y) = P (Y|x)P (x)dx.
Given all the previous elements, a given Bayesian estimator x is associated to the conditional probability distribution :

P (x|Y) = P (Y|x)P (x)

P (Y|x)P (x)

.

A Bayesian estimator x is built by minimizing the Bayes risk noted E(C(x, x)), where C(x, x) is associated to a certain loss function. For instance, the Minimum variance estimator (MV), noted x a M V , considers the Bayes risk as the mean square error,

C(x, x) = ||x -x|| 2 = T r (x -x)(x -x) T ,
where T r is the trace of a given matrix. Therefore, we aim at minimizing the cost function J M V given by the expectation of ||x -x|| 2 knowing Y, which reads after further simplification :

J M V = E ||x -x|| 2 |Y = (x -x) T (x -x)P (x|Y)dx.
(2.5)

The estimator is obtained by nullifying the cost function gradient with respect to x,

∇J M V = 2(x -x)P (x|Y)dx = 0.
The unique estimator that minimizes J M V is given by :

x a M V = xP (x|Y)dx = E(x|Y).

Best Linear Unibiased Estimator (BLUE)

We discuss in this section another approach which is based on the knowledge of the error covariance matrices without specifying first the probability densities P (x|Y) and P (x). We consider the background error η = x tx b associated with the background error covariance matrix

B = E(ηη T ),
(2.6)

and we consider the observation error ε = Y -Hx t associated with the observation error covariance matrix

R = E(εε T ), (2.7)
We assume that the background and observation error are unbiased, i.e. E(η) = 0 and E(ε) = 0. In this context, we seek for the optimal estimate x a defined by the following linear combination :

x a = Lx b + KY.
(2.8)

First and foremost, we assume that the analysis error ε a = x ax t is unbiased. By introducing the definition of the background and observation errors into (2.8) we obtain :

E(ε a ) = E(x a -x t ) = LE(η) + KE(ε) + (L + KH + Id)E(x t ) = 0.
As the background and the observation errors are unbiased, we easily establish the value of L and deduce the expression of the analaysis state x a :

x a = x b + K(Y -Hx b ),
(2.9)

where the vector Y -Hx b is called the innovation vector and the matrix K is referred to as the gain matrix. The latter is unknown and has to be defined along with the analysis covariance matrix noted A. If we assume in addition that x and Y are both Gaussian random vectors and the background and observation errors are non-correlated, i.e. E(εη T ) = 0, then the analysis state (2.8), referred to as the Best Linear Unbiased Estimator (BLUE), is exactly equivalent to the MVE described in § 2.2.2. Under the previous hypotheses, we seek for the gain matrix K which minimizes the variance E(||ε a || 2 ) by using the information provided by covariance matrices (2.6) and (2.7). We consider the following cost function by injecting (2.8) into (2.5) :

J BLU E = E(||x a -x t || 2 ) = (x b -x t +K(Y -Hx b )) T (x b -x t +K(Y -Hx b ))P (x|Y)dx,
then we cancel its gradient with respect to the gain matrix K

∂J BLU E ∂K = (Y -Hx b ) T (x b -x t + K(Y -Hx b ))P (x|Y)dx = 0.
We use the background and observation errors and their proprieties definition to deduce the gain matrix (2.10) and the analysis covariance matrix (2.11). The reader can find the details in the appendice A. The complete procedure to determine the BLUE is summarized in algorithm 1.

K = BH T (R + HBH T ) -1 , (2.10) A = E((x a -x t )(x a -x t )) = B -BH T (HBH T + R) -1 HB (2.11)
Algorithm 1 Best Linear Unbiased Estimator (BLUE)

1: Compute the gain matrix

K = BH T (R + HBH T ) -1
2: Deduce the analysis and its covariance error matrix

x a = x b + K(Y -Hx b ) A = E((x a -x)(x a -x)) = (Id -KH)B
We point out that so far we didn't take into account the state variable dynamics. At this point, two main data assimilation approaches emerged. On the one hand, the sequential approach, presented in § 2.3, couples the results of the Bayesian estimation, in particular the BLUE, with a dynamic model. On the other hand, the variational approach, described in § 2.4, uses the results of optimal control theory by imposing the dynamic model as a constraint and aims at minimizing the least squares cost function written in terms of covariance matrices :

J(x) = (x -x b )B -1 (x -x b ) + (Y -H(x))R -1 (Y -H(x)).

Sequential approach

The sequential approach can be seen as an extension of the Bayesian minimum variance estimation which propagates the background state forward in time by a dynamical model. As the name suggests, the goal of sequential data assimilation consists in correcting the background state each time a new observation is available. We discuss throughout this section the Kalman filter and its extensions.

Kalman filter

The Kalman Filter introduced by [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] can be seen as the direct extention of the BLUE for dynamics model (see § 2.2.3). In this approach, we seek for the analysis state x a k+1 that corrects the background state, that we refer to as the forecast state x f k+1 , given the observation Y k+1 at the instant t k+1 as illustrated in figure 2.1. We consider the following dynamic linear model and observation defined by equations (2.12) and (2.13),

x f k+1 = M k,k+1 x f k , (2.12) Y k+1 = H k x t k+1 + ε k , (2.13)
where the linear discrete operator M k,k+1 describes the flow operator that links the state evolution between the instant t k and t k+1 . This model imperfections are represented by the model error q k associated to the error covariance matrix Q k such as,

x t k+1 = M k,k+1 x t k + q k .
In the same way, the observation operator H k is assumed linear and the error observation ε k+1 at the instant t k+1 is associated to the observation covariance matrix R k+1 .

We observe immediately that the forecast error is propagated through the following forecast covariance matrix P f k+1 ,

P f k+1 = E((x f -x t ) T (x f -x t )) = M k,k+1 P a k M T k,k+1 + Q k .
The Kalman filter approach considers the same hypotheses as those of the BLUE given in § 2.2.3, and assumes in addition that the model and observation errors aren't correlated, in other words,

E(ε k q T k ) = 0.
The analysis state x a is also associated to an error covariance matrix P a ,

P a = E((x a -x t )(x a -x t ) T ).
Under all the previous hypotheses, the Kalman Filter algorithm, given by algorithm 2, consists into an analysis step, similar to the BLUE algorithm, and an additional forecast step, that invokes the dynamical model M k,k+1 and computes the forecast covariance matrix.

Algorithm 2 Kalman filter

1: Analysis step similar to the BLUE algorithm

K = (H k P f k ) T H k (H k P f k ) T + R k -1
(2.14)

x a k = x f k + K k (Y -Hx f k ) (2.15) P a k = (Id -K k H k )P f k (2.

16)

2: Forecast step performed by using the model

x f k+1 = M k,k+1 x a k
(2.17)

P f k+1 = M k,k+1 P a k M T k,k+1 + Q k (2.18)

Extended Kalman Filter

The Kalman filter can be extended to non-linear dynamics and observation operators. In this case, we consider the linear tangent operator around the analysis state x a k , noted ∂ x M k , of the non-linear dynamics operator M k,k+1 and the linear tangent operator around x f k , noted ∂ x H k , of the non-linear observation operator H k+1 . The tangent linear operators are obtained by the application of the Taylor series :

M k,k+1 (x a (t k )) = M k,k+1 (x f (t k )) + ∂ x M k (x a (t k ) -x f (t k )) + O( ε k | 2 ), H k (x a (t k )) = H k (x f (t k )) + ∂ x H k (x a (t k ) -x f (t k )) + O( ε k | 2 ).
This extension, called the Extended Kalman Filter (EKF), is illustrated by algorithm 3.

The major drawback of the classic and extended Kalman filter lays in the computation of the Kalman gain filter matrix K and the error covariance matrices P a k and P f k , which requires a significant computation and storage resources of the squared size of the state space vector. Therefore, these techniques are only usable for systems of low size which isn't suitable for NWP, geophysics and in particular, fluid mechanics applications. Besides, we point out that the EKF provides a good estimation for slightly nonlinear models, which is restrictive for many applications. If we consider for instance the Navier Stokes equations as the dynamics model of interest, its associated tangent linear model operator is likely to neglect the non-linear effects and compromise the optimality of the analysis.

Algorithm 3 Extended Kalman filter

1: Forecast step 2: Compute the forecast using the nonlinear dynamic operator

x f k+1 = M k,k+1 x a k (2.

19)

3: Compute the forecast covariance matrix P f k+1 depending on the tangent linear and adjoint dynamic model

P f k+1 = ∂ x M k P a k ∂ x M * k + Q k (2.

20)

4: Analysis step 5: Compute the Kalman gain matrix and deduce the analysis state

K k+1 = P f k+1 ∂ x H * k ∂ x H k P f k+1 ∂ x H * k + R k -1
(2.21)

x a k+1 = x f k+1 + K k+1 (Y k+1 -∂ x Hx f k+1 ) (2.

22)

6: Update the analysis error covariance matrix P a k+1

P a k+1 = (Id -K k+1 ∂ x H k+1 )P f k+1 (2.23)

Ensemble Kalman Filter

The Ensemble Kalman Filter (EnKF), first introduced in Evensen (1994), proposes an efficient alternative to the EKF using a Monte-Carlo approach. In this context, we consider a set, called ensemble, of possible state samples, each referred to as an ensemble member, noted {x f,(1) (t + 1), • • • , x f,(N ) (t + 1)}. This ensemble is generated from the randomization of the non linear dynamics model given the result x a (t) of a previous assimilation. Instead of directly computing and storing the error covariance matrices, the EnKF uses an empirical expression of the ensemble mean and covariance matrix to deduce a low-rank approximation of the gain matrix and the covariance matrices P f and P a . The EnKF process is illustrated in figure 2.2.

The more ensemble members are employed, the more accurate the approximations are and the more computational resources is required. A balance policy must be clearly chosen. However, ours is always reduced to choose a number of ensemble members N much smaller than the size of the state of interest, noted n. Similarly to the sequential data assimilation techniques presented throughout this section, the EnKF consists in the following forecast and analysis steps.

Ensemble forecast

In the EnKF context, the ensemble forecast consists in generating and propagating the ensemble of N ensemble members by the nonlinear dynamics. Each ensemble member is obtained by perturbing the given forecast trajectory associated to the initial forecast state x f (t 0 ). As we deal with nonlinear dynamics, each ensemble

f,(1) k+1 , • • • , x f,(N )
k+1 } (blue dots). The latter is generated from the randomization of the non linear dynamics model given the result x a k of a previous assimilation. The analysis ensemble is obtained by correcting the forecast ensemble (red dots) with the observation Y k+1 (green dot) and is associated with the analysis covariance error matrix P e (red ellipse). member can be generated by adding a perturbation δ (i) to the initial forecast, such that

x f,(i) (t 0 ) = x f 0 + δ (i) , ∀i ∈ {1, • • • , N },
and propagating each x f,(i) (t 0 ) by the time integration in [t 0 , t] of the nonlinear model operator M t . We point out however that more and more ensemble systems are generated by perturbing the model parameters, which seems to provide a more accurate ensemble spread as studied in [START_REF] Wei | The impact of initial spread calibration on the relo ensemble and its application to lagrangian dynamics[END_REF] and in [START_REF] Wu | A comparison study of data assimilation algorithms for ozone forecasts[END_REF].

The ensemble members are then gathered into the following matrix,

X f t := x f,(1) (t), • • • , x f,(1) (t) ∈ R n×N ,
and the empirical ensemble mean of the N ensemble members is noted

< X f t >:= 1 N N i=1 x f,(i) t
.

The forecast step ends with the estimation of the following ensemble covariance matrix P e , which implicitly captures the dominant directions of uncertainty of the forecast error,

P f ≈ P e = 1 N (X f t -< Xt f >)(X f t -< X f t >) T .
(2.24)

Ensemble analysis

We saw previously that the analysis procedure consists in computing the gain matrix K k , deducing the analysis state x a k and finally updating the analysis covariance matrix P a k from the forecast P f k . In the ensemble context, there are two main approach to deduce the analysis estimate and evaluate the analysis errors.

In the same vein as the pioneering work of [START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics[END_REF] and taking into account the work of [START_REF] Burgers | Analysis scheme in the ensemble kalman filter[END_REF], a first approach consists in generating an ensemble of measurement samples in the same way as the ensemble forecast (x

f,(i) t ) i=1,••• ,N .
Therefore, each ensemble observation is defined as

Y (i) t = Y t + ε (i) t , i = 1, • • • , N.
(2.25)

At this point, the most intuitive process consist to replace P f k by the matrix P e (2.24) to compute the Kalman gain matrix K,

K t = P e H(HP e H T + R) -1 , (2.26)
In this context, the observation operator is nonlinear, therefore the terms P e H and HP e H T are approximated by

P f H = 1 N -1 (X f t -< X f t >)(H(X f t )-< H(X f t ) >) T , (2.27)
HP f H T = 1 N -1 (H(X f t )-< H(X f t ) >)(H(X f t )-< H(X f t ) >) T .
(2.28)

We then deduce the analysis state

x a,(i) t = x f,(i) t + K t (Y (i) t -H(x f,(i) t
).

(2.29)

At last, the update of the analysis matrix (2.23) is replaced by the mean of the variance of the ensemble forecast errors given by

P a,(i) t = 1 N -1 (X a t -< X a t >)(H(X a t )-< H(X a t ) >) T .
(2.30)

Variational approach

The variational approach was first introduced by [START_REF] Sasaki | An objective analysis based on the variational method[END_REF] and expanded rapidly in the 80s with the works of Le Le [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations : theoretical aspects[END_REF] and [START_REF] Courtier | Variational assimilation of meteorological observations with the adjoint vorticity equation. ii : Numerical results[END_REF]. In opposition to the sequential approach, the variational approach considers a given assimilation window in which we have a so-called background trajectory, generated by the dynamics from an initial background state x b 0 , and a set of observations Y. The goal of this approach consists in changing gradually the initial background state x b 0 to draw the trajectory closer to the observations by solving an optimal control problem as described by [START_REF] Lions | Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles[END_REF][START_REF] Lions | Optimal control of systems governed by PDEs[END_REF]. In practice, this approach comes to minimize an objective functional which measures the mismatch between the model predictions and complete or partial measurements of the system, weighed by inverse covariance matrices. This optimization is conducted under the constraint that the optimal state obeys to the model dynamics. Therefore, we expect the assimilation results to satisfy two basic requirements. On the one hand, the assimilation trajectory has to be close to the observations, ideally to the accuracy of the observations themselves. On the other hand, the trajectory has to be consistent with the physical relations between each variables and the evolution of the flow imposed by the model.

Variational formulation

The data assimilation problem can be formulated as the minimization of a least square cost function expressed in terms of the background and the observation covariance matrices, respectively B and R. In this context, we consider the minimization of such function over an assimilation window [t 0 ; t f ] such as

J(x 0 ) = 1 2 (x 0 -x b 0 ) T B -1 (x 0 -x b 0 ) + 1 2 t f t 0 (H(x t ) -Y(t)) T R -1 (H(x t ) -Y(t)).
A complete integration of the dynamics model on a given assimilation widow [t 0 ; t f ], provides a trajectory of the state, associated to an initial guess x 0 and a flow map ϕ t (x 0 ) such as :

x t := ϕ t (x 0 ) = x 0 + t t 0 M(x(s))ds.
We place the data assimilation problem within the control optimal theory [START_REF] Lions | Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles[END_REF][START_REF] Lions | Optimal control of systems governed by PDEs[END_REF] framework. Formally, let V be a Hilbert space identified to his dual space U. We assume that the state variable x t ∈ V and the observations Y ∈ O are square integrable functions in Hilbert spaces identified to their dual. The norms corresponds to the Mahalanobis distance defined from the inner products < R -1 ., . > O , < B -1 ., . > V of the measurements and the state variable spaces respectively. They involve covariance matrices R, and B related to the measurement error and the error on the initial condition.

The minimization is performed by a gradient-based algorithm which, as its name suggests, requires the evaluation of the cost function gradient ∇J(x 0 ) with respect to the initial state x 0 . The optimization algorithm usually requires the computation of the gradient for several values of x 0 to converge to the optimal solution. Therefore one of the major issue of the variational approach lays in the computation of an accurate numerical gradient. There are basically three approaches to obtain numerically such gradient.

The finite difference schemes consist in approximating the gradient by varying each variables of the cost function. Let's consider a first order finite difference scheme is given by

∇J(x 0 ) ≈ J(x 0 + αd) -J(x 0 ) α ,
where d is an arbitrary direction. In order to meet the accuracy requirements, we need to fix a sufficiently small differentiation step α. Another alternative consists in increasing the numerical scheme order p which costs 2p -1 computation and storage of the cost function. This comes to 2p -1 temporal integration of the nonlinear model for each iteration of the gradient based algorithm. Therefore, this approach is clearly unsuitable to our problem.

The sensitivity equation gives the information on how the flow will react to a given perturbation dx = ∂x ∂η δη. In opposition to the finite difference approach, we don't need several evaluations of the cost function for different inputs. The sensitivity equation is obtained by the linearization of the cost function J around x 0 in the direction δη. The differentiation of the perfect nonlinear equations (2.1) and (2.2) in the direction δη reads :

∂ t dx(t, x) + ∂ x M(x(t, x))dx(t, x) = 0, (2.31) dx 0 (x) = δη(x),
(2.32)

where ∂ x M(x) denotes the tangent linear operator associated to the non-linear operator M defined as :

∂ x M(x)dx = lim β→0 M(x + βdx) -M(x) β .
(2.33)

We deduce the differentiation of the cost function in the direction δη :

∂J ∂η , δη V = B -1 (x 0 -x b 0 ), δη V - t f t 0 R -1 (Y(t) -H(x t )), ∂ x H(x t ) ∂x ∂η δη O dt.
(2.34) Introducing the adjoint of the linear tangent observation operator (∂ x H(x)) * , defined as : (2.35) the relation (2.34) can be reformulated as :

∀(u, v) ∈ (V, O), < (∂ x H(x)) u, v > O =< u, (∂ x H(x)) * v > V ,
∂J ∂η , δη V = B -1 (x 0 -x b 0 ), δη V - t f t 0 (∂ X H) * R -1 (Y(t) -x t ), ∂x ∂η δη V dt.
(2.36) Equation (2.36) describes the functional gradients in the directions δη. For the whole gradient, if we consider p the number of design variables, we have to solve p sensitivity equations. In the framework of our PhD, for instance, we ought to control all three components of the velocity fields in a p ≈ 10 6 volume, therefore this approach is too expensive.

The adjoint formulation is an elegant alternative based on the work of [START_REF] Lions | Optimal control of systems governed by PDEs[END_REF] and Le Le [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations : theoretical aspects[END_REF]. We consider a so-called adjoint variable λ ∈ V and carry out the integration over the range [t 0 , t f ] of the inner product between λ and the relation (2.31) such as :

t f t 0 ∂dx ∂t (t), λ(t) V dt + t f t 0 ∂ x Mdx(t), λ(t) V dt = 0.
(2.37)

An integration by parts of the first term yields :

-

t f t 0 - ∂λ ∂t (t) + (∂ x M) * λ(t), dx(t) V dt = λ(t f ), dx(t f ) V -λ(t 0 ), dx(t 0 ) V ,
(2.38) where the adjoint of the tangent linear operator (∂ x M) * : V → V has been introduced. At this point no particular assumptions nor constraints were imposed on the adjoint variable. However, we are free to specify the set of adjoint variables of interest by setting a particular evolution equation or a given boundary conditions simplifying the computation of the functional gradient. As we will see, imposing the adjoint variable λ as the solution of the system

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

As a matter of fact, injecting this relation into equation (2.38) with dx(t 0 ) = δη and dx = (∂x/∂η)δη allows us to identify the right hand second terms of the functional gradients (2.36) and we get :

∂J ∂η , δη V = -λ(t 0 ), δη V + B -1 (δx(t 0 ) -δx 0 ), δη V .
(2.40)

From these relations, one can now readily identify the expression of the cost function derivative with respect to the control variable

∂J ∂η = -λ(t 0 ) + B -1 (δx(t 0 ) -δx 0 ).
(2.41)

In the adjoint-based approach, we only require to solve the adjoint equation once. We note that if the nonlinear dynamics equation is time dependent, then the adjoint equation is time dependent and is solved backwards in time.

Standard 4DVar technique

The aim of the standard 4DVar technique consists in determining the initial optimal state x a (t 0 ) from the so called background trajectory, associated to an initial guess x b (t 0 ) and generated by the dynamics model, and a set of observations defined within a given assimilation window. The optimal state leads to an analysis trajectory which is a compromise between the background trajectory and the observations. The standard 4DVar technique is illustrated by figure 2.3.

A complete integration of the dynamics model on a given assimilation widow [t 0 ; t f ], provides a trajectory of the state, referred to as x t , associated to an initial guess x 0 and a flow map ϕ t (x 0 ) such as :

t := ϕ t (x 0 ) = x 0 + t t 0 M(X(s))ds.
We ought to minimize the distance between the background trajectory defined by (2.1) and (2.2), and the observations (2.3), by controlling the initial state x 0 , as formulated by the following cost function :

J(x 0 ) = 1 2 x b 0 -x(t 0 , x) 2 B -1 + 1 2 t f t 0 Y(t) -H(x t ) 2 R -1 dt, (2.42)
where the L 2 norm is defined as

f 2 A = Ω f (x)A(x, y)f (y)dxdy.
The associated minimization problem is referred to in the literature as the strong constraint variational assimilation formulation.

The background and observation covariance matrices, respectively B and R, define the variation of the background and observations through the variance and the dependence of each variable to one another through the covariance. The inverted matrices used in the cost function aim at normalizing each term and define the weight of the observation with respect to the background. These matrices guide the minimization process by determining which data should be given more confidence at each mesh point of the grid.

In order to solve the optimal control problem defined by (2.42), (2.1), (2.2) and (2.3), we resort to an iterative optimization method. The literature provides a large range of iterative methods to solve an optimization problem under constraint. In this thesis, we opt for the limited storage gradient-based optimization method LBFGS proposed by [START_REF]Updating quasi-newton matrices with limited storage[END_REF].

We now have all elements in hand to establish the standard 4DVar assimilation algorithm 4.

Algorithm 4 Standard Variational assimilation technique

Initialization: k = 0 and x k 0 = x b 0 1: repeat 2:
Compute x t with the forward integration (2.1) and the cost function (2.42)

3:

Compute the initial adjoint variable λ(t 0 ) by the backward integration (2.39)

4:

Update the initial state

x k+1 0 = x k 0 + Bλ(t 0 ) 5: until x k+1 0 -x k 0 < tol 6: From x a 0 = x k+1
0 deduce the analysis trajectory x a t

Incremental Variational Assimilation

As the previous standard 4DVar assimilation cost function (2.42) implies the use of the nonlinear model, there is a chance that the optimization procedure leads to a local minima. Whereas, if we deal with a convex cost function, we are guaranteed to find the global minima. It is in this perspective that [START_REF] Courtier | A strategy for operational implementation of 4d-var, using an incremental approach[END_REF] introduced the incremental strategy to improve the 4DVar technique. This optimization strategy consists in performing a linearization of the dynamics around the current trajectory and minimizes a convex cost function with respect to an increment δx 0 instead of the initial state x 0 . The optimal increment is then used to update the initial state and we carry out the linearization and optimization again for a fixed number of iterations. The incremental 4DVar strategy is illustrated in figure 2.4. This procedure is similar to a Gauss-Newtoon minimization.

The initial increment δx 0 is propagated by the following tangent dynamics equation linearized around a given state x,

∂ t δx(t, x) + (∂ x M)δx(t, x) = 0, (2.43) δx 0 (x) = (x b 0 (x) -x 0 (x)) + η(x), (2.44)
where (∂ x M) is the tangent operator associated to the nonlinear dynamics operator M (2.33). Consequently, the cost function is expressed in terms of the initial increment δx 0 ,

J(δx 0 ) = 1 2 δx 0 2 B -1 + 1 2 t f t 0 H(x t ) -Y(t) 2 R -1 dt, (2.45)
where we consider the perturbed trajectory

x t = ϕ(x b t ) + ∂ x ϕ t (x b t )δx 0 .
The cost function (2.45) can be alternatively formulated as : (2)

J(δx 0 ) = 1 2 δx b 0 -δx 0 2 B -1 + 1 2 t f t 0 ∂ x Hδx(t, x) -D(t, x) 2 R -1 dt, (2.46)
0 = x (1) 0 + δx (1)
0 . We repeat the procedure until we reach the optimal state x 0 which minimizes the nonlinear cost function J (red curve). where

x t = ∂ x ϕ t (x b t )δx 0 , and the innovation vector D(t, x) is defined by

D(t, x) = Y(t, x) -H(ϕ t (x b t )).
As in the non-incremental case, we perform a gradient descent optimization with the LBFGS algorithm and we compute the incremental cost function gradient with the adjoint equation which reads in this case :

-∂ t λ(t) + (∂ x M) * λ(t) = (∂ x H) * R -1 (∂ x H∂ x ϕ(x 0 )D(t, x)) λ(t f ) = 0.
(2.47)

The gradient is then deduced by :

∂ δx 0 J(δx 0 ) = -λ(t 0 ) + B -1 δx 0 .
(2.48)

Once we determine the optimal increment, we update the former initial state x 0 and deduce its corresponding trajectory by solving the nonlinear dynamics equations. The updated trajectory can then be used as a linearization state and we can perform a second outer loop optimization around the new state. The incremental variational assimilation technique proposed by [START_REF] Courtier | A strategy for operational implementation of 4d-var, using an incremental approach[END_REF] is illustrated by algorithm 5.

In the standard case we seek to update directly the background trajectory and in the incremental approach we only use the tangent procedure to update an incremental trajectory. The nonlinear dynamics is used for the outerloop only.

Algorithm 5 Incremental 4DVar assimilation technique [START_REF] Courtier | A strategy for operational implementation of 4d-var, using an incremental approach[END_REF] Initialization:

x

(1) 0 = x b 0
Forward integration of the nonlinear dynamics (2.1)

for k = 1, N k do
Set the initial increment as : δx

(k) 0 = δx b 0 repeat
Compute δx t with the tangent linear dynamics (2.43) and (2.44) and the cost function

Compute the adjoint variable λ(t 0 ) with the backward integration of relation (2.47)

Update the initial increment value δx

(k) 0 = Bλ(t 0 ) until x k+1 0 -x k 0 < tol Update the background : x (k+1) 0 = x (k) 0 + δx (k) 0 end for From x a 0 = x (k+1) 0
deduce the analysis trajectory x a t

Inflow and initial conditions control

In the framework of this thesis, we aim at reconstructing the downstream of a cylinder wake by combining a DNS model with three velocity components measurements of the flow in the inflow and in the stream-wise planes (see § 4.1). In this configuration, noisy or incomplete inlet conditions (with no or erroneous small and even large scales information) may lead to an unrealistic background trajectory which might compromise the convergence of the variational assimilation technique employed. This issue can be tackled by gradually modifying the inflow condition at each time step by including the latter as a control variable in the variational assimilation formulation. This approach has been successfully carried out by [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF] to generate inflow and initial conditions of a 2D turbulent flow. The idea of controlling the inflow was also carried out in [START_REF] Habert | Estimation of lateral inflows using data assimilation in the context of real-time flood forecasting for the marne catchment in france[END_REF] with an Extended Kalman Filter in the context of real-time flood forecasting.

Within this context, we introduce the control parameter u(t) which is defined in the space of square integrable function on spatio-temporal domain U T . Hence, we consider the data assimilation system,

∂ t x(t, x) + M(x(t, x), u(t)) = 0,
(2.49)

x 0 (x) = x b 0 (x) + η(x).
(2.50)

In addition to the background error η(x) and observation error ε(t, x), associated to the covariance error matrices B and R respectively, we consider the deviation θ(t) between the control u(t) an its a priori value u c (t) associated to the covariance error matrix B c . Therefore, our goal consists in minimizing the cost function (2.51) with respect to the control variable γ(t,

x) = {η(x), θ(t)} = {x(x, t 0 ) -x 0 (x), u(t) - u c (t)}. J(γ) = 1 2 t f t 0 Y(t) -H(x t ) 2 R -1 dt + 1 2 x(t 0 , x) -x b 0 (x) 2 B -1 1 2 t f t 0 u(t) -u c (t) 2 B -1 c dt, (2.51) with f 2 B -1 c = Γ f (x)B -1 c (x, y)f (y)
dxdy where Γ refers to the inflow section. Given a perturbation dx = ∂x ∂u δu(t) + ∂x ∂η δη, the differentiation of the equations (2.49) and (2.50) in the direction δγ = (δu, δη) reads :

∂ t dx(t, x) + ∂ x M(x(t, x), u(t))dx(t, x) + ∂ u M(x(t, x), u(t))du(t) = 0, (2.52) dx 0 (x) = δη(x).
(2.53)

The differentiation of the cost function in the direction δγ :

∂J ∂u , δu U T = t f t 0 B -1 c (u(t) -u c (t)), δu(t) U - t f t 0 R -1 (Y(t) -H(x t )), ∂ x H(x t ) ∂x ∂u δu(t) O dt, (2.54) ∂J ∂η , δη V = B -1 (x 0 -x b 0 ), δη V - t f t 0 R -1 (Y(t) -H(x t )), ∂ x H(x
U T = t f t 0 B -1 c (u(t) -u c (t)), δu(t) U - t f t 0 (∂ x H(x t )) * R -1 (Y(t) -H(x t )), ∂x ∂u δu(t) V dt (2.56) ∂J ∂η , δη V = B -1 (x 0 -x b 0 ), δη V - t f t 0 (∂ x H(x t )) * R -1 (Y(t) -H(x t )), ∂x ∂η δη V dt. (2.57)
As it is unfeasible to compute directly the cost function gradients, we resort to a similar adjoint approach than in § 2.4.1. To that end, we proceed to the integration over the range [t 0 , t f ] of the inner product between an adjoint variable λ ∈ V and the relation (2.52) is given by :

t f t 0 ∂dx ∂t (t), λ(t) V dt + t f t 0 ∂ x Mdx(t), λ(t) V dt + t f t 0 ∂ u Mδu(t), λ(t) V dt = 0.
(2.58)

An integration by parts of the first term yields :

-

t f t 0 - ∂λ ∂t (t) + (∂ x M) * λ(t), dx(t) V dt = λ(t f ), dx(t f ) V -λ(t 0 ), dx(t 0 ) V + t f t 0 δu(t), (∂ u M) * λ(t) U dt.
(2.59)

As in § 2.4.1, we impose the adjoint variable λ as the solution of (2.39) and inject the latter into (2.59) with dx(t 0 ) = δη and dx = (∂x/∂η)δη to identify the right hand second terms of the functional gradients (2.56) and (2.57). Hence, we get ∂J ∂u , δu

U T = t f t 0 δu(t), B -1 c (u(t) -u c (t)) + (∂ u M) * λ(t) U dt = B -1 c (u -u c ) + (∂ u M) * λ, δu U T ∂J ∂η , δη V = -λ(t 0 ), δη V + B -1 (x 0 -x b 0 ), δη V .
Thus, we deduce the expression of the cost function derivative with respect to δγ = (δη, δu) :

∂J ∂u = -λ(t 0 ) + B -1 c (u -u c ) + (∂ u M) * λ, (2.60)
∂J ∂η = -λ(t 0 ) + B -1 (x(t 0 ) -x 0 ). (2.61)
All in all, the optimal control problem defined by (2.51), (2.49), (2.50) and (2.3) is solved by the LBFGS method (see § 2.4.1), the assimilation algorithm is summarized in 7.

Algorithm 6 4DVar assimilation with inflow and initial condition control [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF] Initialization: k = 0 and x k 0 = x b 0 repeat Compute x t with the forward integration (2.49) and the cost function (2.51) Compute the initial adjoint variable λ(t 0 ) by the backward integration (2.39) Update the initial state

x k+1 0 = x k 0 + Bλ(t 0 ) Update the parameter model u k+1 = u k + B c (∂ u M) * λ(t 0 ) until x k+1 0 -x k 0 < tol From x a 0 = x k+1
0 deduce the analysis trajectory x a t

Incremental variational assimilation with inflow and initial condition control

In the framework of this thesis, we implemented the incremental version of the previous variational formulation which includes the control of the inflow and the initial conditions (see § 2.4.4). In the same vein as in § 2.4.3, we aim at minimizing the following convex cost function with respect to an increment δγ = {δx, δu}

J(δγ 0 ) = 1 2 δx 0 2 B -1 + 1 2 t f t 0 δu(t) 2 Bc -1 dt + 1 2 t f t 0 H(x t ) -Y(t) 2 R -1 dt, (2.62)
where the increment vector δγ verify (2.52) and (2.53).

As in the previous case, we call upon the adjoint equation (2.39) to compute the incremental cost function gradient with respect to the increment vector δγ such that

∂ δx 0 J = -λ(t 0 ) + B -1 δx 0 , (2.63) ∂ δu J = -λ(t 0 ) + B -1 c δu + (∂ u M) * λ. (2.64)
The present 4DVar technique is described in algorithm 7 and is referred to as the 4DVar technique in the following chapters.

Algorithm 7 Incremental 4DVar assimilation technique with inflow and initial condition control

Initialization: x (1) 0 = x b 0 Forward integration of the nonlinear dynamics (2.49) for k = 1, N k do
Set the initial increment as : δx Update the initial increment value δx

(k) 0 = Bλ(t 0 ) until x k+1 0 -x k 0 < tol Update the background : x (k+1) 0 = x (k) 0 + δx (k) 0 Update the parameter model u k+1 = u k + B c (∂ u M) * λ(t 0 ) end for From x a 0 = x (k+1) 0
deduce the analysis trajectory x a t

Numerical adjoint construction strategies

As described previously, the variational assimilation techniques employs a gradientdescent method to solve an optimal control problem. The computation of the cost function gradient is deduced from the backward integration of the adjoint equation (2.39). We discuss in this section the construction and implementation of the numerical adjoint equation corresponding to the dynamics.

There are two main strategies for the construction of a numerical adjoint associated to a given dynamics model. On the one hand, the differentiate-then-discretize strategy consists in determining the analytic adjoint from the nonlinear dynamics model, then discretizing the adjoint model to obtain a numerical adjoint procedure, referred to as the continuous adjoint. On the other hand, discretize-then-differentiate consists in constructing the discrete adjoint directly from the discretized non-linear dynamics model. The latter can be performed by the mean of an Automatic Differentiation tool (AD), which is described in appendix C.1. Both strategies are illustrated in figure 2.5 and were compared in different contexts [START_REF] Nadarajah | A comparison of the continuous and discrete adjoint approach to automatic aerodynamic optimization[END_REF][START_REF] Griesse | Evaluating gradients in optimal control : continuous adjoints versus automatic differentiation[END_REF]. In general, as iterative gradient-descent algorithms are sensible to numerical errors, we seek for the most accurate numerical adjoint associated to a given problem. The discretize-then differentiate strategy provides an accurate numerical adjoint operator up to the machine precision, whereas the continous adjoint obtained by the differentiate-the-discretize strategy is accurate up to the discretization errors of the numerical scheme on a given mesh. Therefore, if we possess beforehand a sufficiently accurate original code, the discrete adjoint will remain at the same accuracy up to the machine precision. In addition, the continuous case generates boundary terms, commonly by integration by parts, that have to be dealt separately. On the other hand, the construction of discrete adjoint is straightforwardly derived from the original code and doesn't require additional thoughts to discretizations or boundary conditions for either operator.

The variational approach has been applied so far on 2D problems which require a reasonable computational cost, and most of the authors opted for the construction of an continuous adjoint [START_REF] Lemke | Adjoint based pressure determination from piv-data -validation with synthetic piv measurements[END_REF][START_REF] Mons | Is isotropic turbulence decay governed by asymptotic behavior of large scales ? an eddy-damped quasi-normal markovian-based data assimilation study[END_REF] Within the framework of this thesis, we build the discrete adjoint by applying the AD on the highly accurate 3D DNS code Incompact3d (see § 4.2). Unlike the continuous adjoint, the discrete adjoint takes directly advantage of the Incompact3d's optimal parallelization structure and conserves the high accurate of the original code. We also we use the same strategy to build the discrete adjoint associated to the 2D SWE model for the reconstruction of the free surface height and velocity fields (see § 3.3). It is indeed difficult to argue this choice as the analytic adjoint is often different from the real one.

The AD implementation can also be eased by the use of a automatic differentiation (AD) tools. In terms of AD implementation, the literature offers several AD tools that employs either the source code transformation (SCT) or the operator overloading (OO) strategy (see appendix C.1.1). In this thesis, we opt for TAPE-NADE which is an AD tool developed by the Tropics team in Inria Sophia Antipolis (Hascoët, 2004). This tool employs the SCT strategy and was formerly successfully used in the work of [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF].

Validation of the adjoint procedure

The discrete tangent and adjoint procedure built by the means of the AD has to be validated by a series of tests. Before presenting the different validation tests carried out, we point out the importance of the choice of the inputs into each validation process. Indeed, in the strictly numerical sense, the validation test should be verified for a random value of the state. However, in practice the dynamics model obey to a physical model, and it has been tested that inputting a random state x to the original direct model f leads often to a numerical explosion. Therefore, in order to validate a dynamic model on several time steps, we have to input state x has to obey to the physics. We detail in the following validation tests the degrees of freedom that are adopted to this respect.

Tangent validation

The tangent validation is based on the Taylor formula applied to the direct method operator f around a given state x 0 given a direction d,

f (x 0 + αd) ≈ f (x 0 ) + α < ∇f (x 0 ), d > .
(2.65)

This validation test consists in studying the ratio

R 1 := f (x 0 + αd) -f (x 0 ) α < ∇f (x 0 ), d > .
(2.66)

It is particularly interesting if we obtain the tangent operator ∇f (x 0 ) with an AD tool.

We point out that the ratio accuracy depends on the choice of α. A too large value of α might result to accuracy errors as the dynamics operator f is non-linear, whereas a too small value of may lead to the accumulation of rounding errors. In the applications within this thesis, we study the ratios for α in a range of α ∈ [10 -3 ; 10 -9 ]. We can reasonably expect the ratio R 1 to be more or less close to 1 up to O(α). In the latter case, it would mean that the tangent procedure would have the same accuracy that a finite difference which isn't suitable.

Dot product

The results of the reverse adjoint mode are validated by using the previously validated tangent. We recall the definition of the adjoint operator (∂ x M) * :

< (∂ x M) Ẋ, Ẏ >=< Ẋ, (∂ x M) * Ẏ >, ∀ Ẋ, Ẏ (2.67)
We set in particular Ẏ = (∂ x M) Ẋ and define X :

= (∂ x M) * (∂ x M) Ẋ = (∂ x M) * Ẏ .
We then develop the dot product of X and Ẋ :

< Ẋ, X >=< Ẋ, (∂ x M) * Ẏ >=< (∂ x M) Ẋ, Ẏ >=< Ẏ , Ẏ > . (2.68)
This comes to reduce the following ratio to 1 :

R 2 := < (∂ x M)X, (∂ x M)X > < X, (∂ x M) * Y > = < Ẏ , Ẏ > < Ẋ, Ȳ > ≈ 1.
(2.69)

The adjoint operator (∂ X M * ) has to be the most accurate as possible. In other words, we have to obtain a ration of R 2 = 1+O( m ) where M represents the machine accuracy.

Adjoint validation

In the 4DVar context, we will particularly evaluate the theoric equality between the adjoint variable λ(t 0 ) at time t 0 and the cost function gradient ∇J(x 0 ) at time t 0 for a given direction d, such that

< ∇J(x 0 ), d >=< λ(t 0 ), d > .
(2.70)

We recall that the adjoint variable λ(t 0 ) is obtained after the backwards adjoint model integration. Similarly to the tangent validation, we apply the Taylor formula to the cost function J and evaluate the ratio (2.71) given by :

R 3 = J(x 0 + αd) -J(x 0 ) α < ∇J(x 0 ), d > ≈ J(x 0 + αd) -J(x 0 ) α < λ(t 0 ), d > .
(2.71)

In analogy to the accuracy discussions for the ratio R 1 , the value of the ratio R 3 depends on the value of α. Although, we have to keep in mind that the first order finite difference scheme provides a rough approximation of the cost function gradient. Therefore, we expect don't expect the ratio to be accurate up to 10 -α .

Hybrid approaches

During the last decade, there is a trend in the development of hybrid data assimilation techniques which combine the advantages of ensemble methods and varitional assimilation approaches. Each technique differs according to the context in which they are developed, therefore there are two main groups of hybrid techniques :

A first group of hybrid techniques, based on the variational formalism, aims at incorporating an ensemble-based background covariance into a variational system. [START_REF] Hamill | A hybrid ensemble kalman filter-3d variational analysis scheme[END_REF] proposed the EnKF-3DVar method which basically combines several 3DVar data assimilation process with the EnKF update process. The background error covariace is set as a linear combination of the ensemble covariance and the static covariance. This idea was extended to the 4DVar formalism by [START_REF] Lorenc | The potential of the ensemble kalman fi filter for nwpâĂŤa comparison with 4d-var[END_REF] who also introduced a localization procedure. Several authors proposed hybrid methods that express either the background error covariance as a linear combination of the ensemble covariance and the static covariance or the solution as a linear combination of the square root of the ensemble-based covariance [START_REF] Liu | An ensemble-based four-dimensional variational data assimilation scheme. part ii : Observing system simulation experiments with advanced research wrf (arw)[END_REF]Buehner et al., 2010a,b;[START_REF] A M Clayton | Operational implementation of a hybrid ensemble/4d-var global data assimilation system at the met office[END_REF][START_REF] Desroziers | 4denvar : link with 4d state formulation of variational assimilation and different possible implementations[END_REF].

A second group of hybrid techniques, based on the EnKF formalism, incorporates the 4DVar formulation to overcome the limitations of the sequential approach. In the one hand, the standard EnKF assimilates sequentially the observation whenever they are available [START_REF] Evensen | The ensemble kalman filter : Theoretical formulation and practical implementation[END_REF]. The hybrid 4D-EnKF technique proposed by [START_REF] Hunt | Four-dimensional ensemble kalman filtering[END_REF] seeks for an observation operator which could link the state to future or past observations. This idea is then extended in the 4D-LETKF technique introduced by [START_REF] Hunt | Efficient data assimilation for spatiotemporal chaos : A local ensemble transform kalman filter[END_REF] and [START_REF] Fertig | A comparative study of 4d-var and a 4d ensemble kalman filter : Perfect model simulations with lorenz-96[END_REF]. [START_REF] Sakov | Asynchronous data assimilation with the enkf[END_REF] introduced the AEnKF which aims at finding the explicit form of temporal evolution of the increment by propagating the correction along the forecast system trajectory.

On the other hand, the explicit form of the update equations in EnKF (or EKF) are difficult to calculate when dealing with high dimensional problems. A group of hybrid techniques propose to replace the EnKF update by the minimization of an arbitrary cost function by an iterative gradient descent procedure similarly to the 4DVar technique which handles efficiently large dimension problems [START_REF] Zupanski | Maximum likelihood ensemble filter : Theoretical aspects[END_REF][START_REF] Sakov | An iterative enkf for strongly nonlinear systems[END_REF]Solonen et al., 2012).

In the framework of this thesis, we compare an ensemble-based 4DVar, referred to as 4DEnVar, a classical 4DVar and the 4D-LETKF techniques to reconstruct a 2D free surface in chapter 3. The 4DEnVar technique is closely related to the strategy proposed by [START_REF] Buehner | Ensemble-derived stationary and flow-dependent background-error covariances : Evaluation in a quasi-operational nwp setting[END_REF] and [START_REF] Liu | An ensemble-based four-dimensional variational data assimilation scheme. part i : Technical formulation and preliminary test[END_REF]. This technique introduces in its objective function an empirical ensemble-based background error covariance which avoids the use of tangent linear and adjoint model. We describe the 4DEnVar in § 3.2.

Chapitre 3

A first application : reconstruction of a free surface flow

The shallow water equations (SWE) -also called Saint Venant equations -describes the height and the velocity of a thin layer of an inviscid fluid at its free surface. These equations are often used as a prototype in geophysics applications, such as ocean motion, river sedimentation and under certain hypotheses, they can also provide a satisfying atmospheric flow approximation model. The simplicity of both numerical model implementation and experimental setting (for an academical context) is an ideal framework for a first application of the 4DVar technique. In this work, we also seized the opportunity to compare the reconstruction of the dimensional free surface flow characteristics (height and 2D velocity fields) using a classic formulation of the 4DVar and an enhanced version of an ensemble based technique [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF].

We present in § 3.1 the dynamical model used to simulate the free surface of the flow of interest. The data assimilation techniques used throughout this chapter are overviewed in § 3.2. We describe in particular the construction of the adjoint procedure necessary to implement the classical 4DVar technique in § 3.3. In a first step, we consider in § 3.4.1 the reconstruction of the free surface characteristics from partial and complete synthetic observations. In a second step, we consider the real measurements of the evolution of the free surface height observed with a Kinect depth sensor as described in § 3.4.2. We finally draw the conclusions of the comparison between the assimilation techniques considered for this application and layout some perspectives in § 3.5.

Dynamics

Hypotheses and governing equations

As its name suggests, the thickness (or height when we consider a flat bottom topology) of the fluid of interest is neglected with respect to its length and width. Under this hypothesis, the mass conservation implies that the vertical velocity of the fluid is small, hence, we only take into consideration the 2D velocity field u := (u, v) T . In the framework of this chapter, we consider a simplified version of the SWE under the following hypothesis :

-The bottom of the domain is flat and non tilted.

-We neglect the Coriolis effect and the bottom friction.

As introduced in § 1.1.1, the SWE are deduced by the depth-integration of the Navier-Stokes equations and are used to model the height and the velocity field at the free surface of the fluid. The simplified 2D SWE used throughout this chapter reads

∂ t h + ∂ x (hu) + ∂ y (hv) = 0, (3.1) ∂ t (hu) + ∂ x (hu 2 ) + ∂ y (huv) + g∂ x h 2 = 0, (3.2) ∂ t (hv) + ∂ x (huv) + ∂ y (hv 2 ) + g∂ y h 2 = 0.
(3.

F (x) =   hu hu 2 + 1 2 gh 2 huv   and G(x) =   hv huv hv 2 + 1 2 gh 2   . The conservative form is
given by :

∂ t x(t, x) + ∂ x F (x(t, x, y)) + ∂ x G(x(t, x, y)) = 0. (3.4)
The literature provides many numerical schemes that can be used to solve the SWE model (see § 1.1.3). In this work we created from scratch a numerical code based on the work of [START_REF] Vreugdenhil | Numerical methods for shallow-water flow[END_REF] who studied numerical solutions for 2D SWE. This code is described thereafter.

Spatial discretization and Roe solver

As described in § 1.1.3, the Finite Volume Method (FVM) is well suited to simulate conservation laws and is usually favored to solve CFD problems. We describe in the present section the discretization of the advection terms F (x) and G(x) of the conservative equation (3.4).

Finite volume formulation

The solution domain Ω is a partition of many control volumes, that are referred to as cells C ij = [x i-1/2 ; x i+1/2 ] × [y j-1/2 ; y j+1/2 ] with ∆x i = x i+1/2 -x i-1/2 and ∆y j = y j+1/2 -y j-1/2 where i ∈ {1, . . . , nx} and j ∈ {1, . . . , ny}, as illustrated in figure 3.1.

Within our context, we assume that each cells have the same size, i.e. ∆x := ∆x i = ∆y j , ∀i, j. On a given cell C ij , the integral form of (3.4) reads :

d dt C ij x(t, x) + 1 ∆x F (x(t, x i+1/2 , y j )) -1 ∆x F (x(t, x i-1/2 , y j )) + 1 ∆y G(x(t, x i , y j+1/2 )) -1 ∆y G(x(t, x i , y j-1/2 )) = 0.
(3.5)

Figure 3.1 -Spatial discretization of the FVM grid : the cell C ij = [x i-1/2 ; x i+1/2 ]× [y j-1/2 ; y j+1/2
] is the cell associated to the quantity x i,j . We aim at computing the the flux crossing each inter-cell flux F n i+1/2,j , F n i-1/2,j , G n i,j+1/2 and G n i,j-1/2 in each direction east, west, south and north respectively.

Given the following averages for each cell C ij at a given time t n ,

x n i,j = 1 ∆x∆y C ij x(t n , x, y)dxdy, F n i+1/2,j = 1 ∆t t n+1 t n F (x(t n , x i+1/2 , y j ))dt, F n i-1/2,j = 1 ∆t t n+1 t n F (x(t n , x i-1/2 , y j ))dt, G n i,j+1/2 = 1 ∆t t n+1 t n G(x(t n , x i , y j+1/2 ))dt, G n i,j-1/2 = 1 ∆t t n+1 t n G(x(t n , x i , y j-1/2 ))dt,
the finite volume formulation reads

x n+1 i,j = x n i,j - ∆t ∆x F n i+1/2,j -F n i-1/2,j - ∆t ∆y G n i,j+1/2 -G n i,j-1/2 . (3.6)
Godunov (1959) proposed a conservative numerical scheme which assumes that the finite volume solution is piece-wise constant at each cell, as illustrated in figure 3.2. Under this hypothesis, we seek for an estimation of the inter-cell fluxes F n i+1/2,j , F n i-1/2,j , G n i,j+1/2 and G n i,j-1/2 by solving a local Riemann problem at each inter-cell.

1D Riemann problem

Given a time interval [t n ; t n+1 ], the 1D Riemann problem at a the inter-cell (x i+1/2 , y j ) is formulated as the following initial value problem :

∂ t x + A.∂ x x = 0, x(t n , x, y j ) = x L := x n i,j if x < x i+1/2 x R := x n i+1,j if x > x i+1/2 , (3.7)
where A = ∇F n i+1/2,j is the Jacobian matrix of F n i+1/2,j .

C i = [x i-1/2 ; x i+1/2 ].
If the flux F n i+1/2,j is linear, the solution is straightforward. We perform a diagonalization A = R -1 ΛR, where the eigenvalues e 1 < • • • < e m correspond to the linearly independent eigenvectors r 1 , • • • , r m of the Jacobian A. Each wave p ∈ {1, • • • , m} carries a jump discontinuity with a velocity propagation e p , also called characteristic speed. By introducing an auxiliary variable W := Rx, the problem (3.7) comes to a system of m scalar hyperbolic problems :

∂ t W (p) (t, x) + Λ∂ x W (p) (t, x) = 0, for p = 1, • • • , m.
(3.8)

The solution of (3.8) is obtained by introducing the following matrices

A + = RΛ + R -1 , A -= RΛ -R -1 , |A| = R|Λ|R -1 ,
where Λ + and Λ -gather the positive and negative characteristic speed respectively, and |Λ| contains the absolute value of the characteristic speed. Given these matrices, we apply the first-order upwind method to (3.8) for positive an negative characteristic speeds such that

(W + ) n+1 i = (W + ) n i -Λ + ∆t ∆x ((W + ) n i -(W + ) n i-1 ), (W -) n+1 i = (W -) n i -Λ -∆t ∆x ((W -) n i+1 -(W -) n i ).
As W is by the definition the sum of W + and W -, the numerical solution of the hyperbolic equation (3.8) reads :

W n+1 i = W n i -Λ ∆t ∆x (W n i -W n i-1 ) + |Λ| ∆t ∆x W n i+1 -2W n i + W n i-1
2 .

Thus, we deduce the numerical solution of the 1D Riemann problem (3.7) for a linear flux F with respect to its Jacobian matrix A :

x n+1 = x n i -A ∆t ∆x (x n i -x n i-1 ) + |A| ∆t ∆x x n i+1 -2x n i + x n i-1 2 .
(3.9)

In the non-linear case, the literature provides several solvers that approximate the Riemann solution [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF][START_REF] Roe | Approximate riemann solvers, parameter vectors, and difference schemes[END_REF]. The Roe solver, proposed by [START_REF] Roe | Approximate riemann solvers, parameter vectors, and difference schemes[END_REF], is commonly chosen when it comes to solve numerically the SWE. This solver consists in linearizing the Jacobian matrix A ≈ Ã(x L , x R ), called Roe matrix, which has to verify the following proprieties :

-Ã(x L , x R ) is assumed to be constant between the left and the right cells, -Ã(x L , x R ) has to be diagonalizable with real eigenvalues {e p } p associated to linearly independent eign vectors {r p } p . This propriety is consistent with the exact Jacobian and the conservation of

Ã(x R -x L ) = F (x R ) -F (x L ).
Once the linearized Jacobian matrix is obtained, the solution of the 1D Riemann problem follows the linear case (3.9). The Roe matrix associated to the 2D SWE was proposed by [START_REF] Priestley | Roe Type Schemes for the 2-D Shallow Water Equations[END_REF] and also [START_REF] Glaister | Approximate riemann solutions of the shallow water equations[END_REF]. We describe the construction of the Roe solver in the following.

Application of the Roe solver to the 2D SWE

In order to construct the Roe matrix, [START_REF] Priestley | Roe Type Schemes for the 2-D Shallow Water Equations[END_REF] introduced the following parameter vectors :

w L :=   w 1 L w 2 L w 3 L   =   √ h L u L √ h L v L √ h L   , w R :=   w 1 R w 2 R w 3 R   =   √ h R u R √ h R v R √ h R   .
At each cell interface we compute the difference between the left flux and the right flux, denoted by the subscripts L an R respectively,

F (x L ) -F (x R ) = Ã(x L -x R ) = 3 i=1 α k e k r k .
For each variable, we denote ∆w := w L -w r and w := 1 2 (w L +w R ). The eigenvalues of the linearized Jacobian Ã(x R , x L ) are defined by

e 1 = w2 w1 -w2 1 , e 2 = w2 w1 , e 3 = w2 w1 + w2 1 ,
and their corresponding eigenvectors are

r 1 =   w1 w2 -w1 w2 1 w3   , r 2 =   0 0 w 1   , r 3 =   w1 w2 + w1 w2 1 w3   ,
with the coefficients

α 1 = ∆w 1 - w1 ∆w 2 -w2 ∆w 1 2 w1 w2 1 , α 2 = w1 ∆w 2 -w2 ∆w 1 w1 , α 3 = ∆w 1 + w1 ∆w 2 -w2 ∆w 1 2 w1 w2 1 .
Hence, the fluxes along the x-direction are described by :

F n i+1/2,j = 1 2 (F (x n L ) + F (x n R )) - 1 2 3 k=1 α k |e k |r k , for x L = x n i,j and x R = x n i+1,j ,
(3.10)

F n i-1/2,j = 1 2 (F (x n L ) + F (x n R )) - 1 2 3 k=1 α k |e k |r k , for x L = x n i-1,j and x R = x n i,j .
(3.11)

Similarly, in the y-direction, we compute the flux difference between the north flux and the south flux, denoted by the subscripts N an S respectively,

G(x S ) -G(x N ) = Ã(x N -x S ) = 3 k=1 β k e k r k .
For each variable we note ∆w = w N -w S and w = 1 2 (w N + w S ). The eigenvalues of the approximated Jacobian Ã(x N , x S ) are defined by

e 1 = w3 w1 -w2 1 , e 2 = w3 w1 , e 3 = w3 w1 + w2 1 ,
and their corresponding eigenvectors are

r 1 =   w1 w2 w3 -w1 w2 1   , r 2 =   0 0 w 1   , r 3 =   w1 w2 w3 + w1 w2 1   ,
with the coefficients

α 1 = ∆w 1 - w1 ∆w 3 -w3 ∆w 1 2 w1 w2 1 , α 2 = w1 ∆w 3 -w3 ∆w 1 w1 , α 3 = ∆w 1 + w1 ∆w 3 -w3 ∆w 1 2 w1 w2 1 .
Finally, the fluxes along the y-direction are given by :

G i,j+1/2 = 1 2 (G(x L ) + G(x R )) - 1 2 3 k=1 α k |e k |r k , for x S = x ij and x N = x i,j+1 , (3.12
)

G i,j-1/2 = 1 2 (G(x L ) + G(x R )) - 1 2 3 k=1 α k |e k |r k , for x S = x i,j-1 and x N = x ij .
(3.13)

Boundary conditions

As we study a fluid contained in a non-porous rectangular tank (see § 3.4.2), we expect the fluid to return in the opposite direction whenever it hits the walls. This comes to impose a reflecting boundary condition in every direction. In this configuration, we impose the values of the height and velocity field at ghost cells {C 1,j , j = 1, . . . ny}, {C i,1 , i = 1, . . . nx}, {C nx,j , j = 1, . . . ny} and {C i,ny , i = 1, . . . nx} around the boundaries. The height and velocity field in these ghost cells are set to the same height as the boundaries and the same velocity in the opposite direction than the incoming wave respectively. The reflecting boundary conditions are illustrated in figure 3.3.

Time integration

We describe in the present section the formulation of the time integration scheme. The general formulation of Runge-Kutta techniques is given by :

x n+1 = x n + ∆t r i=1 b i k i , k i = f (t n + ∆tc i , x n i + ∆t r j=1 a ij k i ),i = 1, ..., r
where the function f computes the fluxes given a state variable

x n such that f (x n ) = F n i,j+1/2 -F n i,j-1/2 + G n i,j+1/2 -G n i,j-1/2
. The 4th order Runge-Kutta method (RK4) is usually privileged in many CFD applications for its robustness and convergence proprieties. The 3rd order method (RK3) also provides a good compromise between accuracy and computation time. Both techniques were implemented in the SWE code and performed very similar results. Thus, we decided to consider the RK3 method for the implementation of the 4DVar code. The RK3 method is defined by fixing r = 3 and the coefficients

(c 1 , c 2 , c 3 ) = (0, 0.5, 1), (b 1 , b 2 , b 3 ) = (1/6, 2/3, 1/6) and (a ij ) 1≤i,j≤3 =   0 0 0 1/2 0 0 -1 2 0   .

Implementation into MATLAB

We created the SWE code based on the dynamical model described previously. The code computes the trajectory of the non-conservative variable state on a given assimilation window [0; N ∆t] as described in algorithm 8. The auxiliary routines u2q and q2u perform a simple transformation from a given non conservative state variables x = (h, u, v) T to the conservative variables x c := (h, q x , q y ) T = (h, hu, hv) T and vice versa. The RK3 routine performs the time integration procedure which deduces the conservative state x n+1 c at time t n+1 = (n+1)∆t from the previous state

x n c at time t n = n∆t (see § 3.1.4). The latter performs several calls of the boundary conditions routine BC defined by § 3.1.3, and the FVM routine which computes the fluxes at each inter-cells (see § 3.1.2).

Algorithm 8 SWE(h 0 ,u 0 ,v 0 )

Initialization: n = 0, h 0 , u 0 , v 0 while n < N do call u2q(h n ,u n ,v n ,q n x ,q n y ) call RK3(h n ,q n x ,q n y ,h n+1 ,q n+1 x ,q n+1 y ) call q2u(h n+1 ,u n+1 ,v n+1 ,q n+1 x ,q n+1 y ) save (h n+1 ,u n+1 ,v n+1 ) n = n + 1 end while return h N , u N , v N

Ensemble-based 4DVar technique

We compare in this chapter the performance of a purely variational data assimilation technique (see § 2.4.7) and an enhanced version of an ensemble based technique [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF]. The latter is based on the work of [START_REF] Liu | An ensemble-based four-dimensional variational data assimilation scheme. part i : Technical formulation and preliminary test[END_REF] and [START_REF] Liu | An ensemble-based four-dimensional variational data assimilation scheme. part ii : Observing system simulation experiments with advanced research wrf (arw)[END_REF], and is referred to as the 4DEnVar in the following. This technique is specifically defined within the framework of preconditioned incremental variational formulation (see § 3.2.1) while handling an empirical approximation of the background covariance matrix. The low rank approximation of the background covariance matrix is directly inspired from the Ensemble Kalman filter techniques where the covariance terms are empirically estimated from an ensemble of realizations.

Preconditioning of the incremental 4DVar assimilation

We recall the incremental variational assimilation formalism defined by the assimilation system

∂ t δx(t, x) + ∂ x Mδx(t, x) = 0, (3.14) δx(t 0 , x) = x b 0 -x 0 + η, (3.15) Y(t, x) = H(ϕ t (x 0 (x)) + ε(t, x), (3.16)
where ∂ x M(x) denotes the tangent linear operator of the dynamics operator M. The variable x b 0 denotes an arbitrary background condition, x 0 is the initial guess and η is some error of background covariance B. The quantity

ϕ t (x 0 (x)) = x t = x 0 + t t 0 M(x s )ds,
denotes the flow map. The state variables x t and the measurements Y are linked by an observation operator H, up to an observation error ε that is assumed to be a zero mean i.i.d. (independent and identically distributed)Gaussian random field with covariance tensor R. In this study, for simplification purpose, we will consider a linear observation operator set to the identity or to an incomplete identity operator when only a part of the state is observable. The optimal increment δx 0 at time t 0 is obtained by minimizing the following cost function :

J(δx 0 ) = 1 2 δx b 0 -δx 0 2 B -1 + 1 2 t f t 0 ∂ x Hδx(t, x) -D(t, x) 2 R -1 dt, (3.17)
where

x t = ∂ x ϕ t (x b t
)δx 0 , and the innovation vector D(t, x) is defined by

D(t, x) = Y(t, x) -H(ϕ t (x b t )
). In order to minimize the objective function J(δx 0 ) with respect to δx 0 we have to cancel the gradient ∇J(δx 0 ). This comes to solve the following linear system :

∇J(δx 0 ) = B -1 δx 0 + t f t 0 ∂ x ϕ * t (∂ x H) * R -1 (∂ x H∂ x ϕ t δx 0 -D(t))dt = 0, (3.18) ⇔ B -1 + t f t 0 ∂ x ϕ * t (∂ x H) * R -1 ∂ x H∂ x ϕ t dt δx 0 = t f t 0 ∂ x ϕ * t (∂ x H) * R -1 D(t)dt, (3.19)
which can be compactly written as Hδx 0 = b, denoting the left hand operator as the Hessian matrix

H = B -1 + t f t 0 (∂ x ϕ t ) * (∂ x H) * R -1 ∂ x H∂ x ϕ t dt, (3.20)
and the right hand as the vector

b = t f t 0 (∂ x ϕ t ) * (∂ x H) * R -1 D(t, x)dt. (3.21)
The possibly ill-conditioned nature of (3.19) depends on the condition number of the Hessian matrix H (3.20). The condition number is defined by the ratio between the maximal eigenvalue over the minimal eigenvalue of the matrix of interest. The larger the condition number, the more sensitive the system is to errors both in the b vector and in the estimate. A poorly-conditioned variational incremental system can be improved by applying a change of variable with the matrix square-root of the background error covariance matrix B. This change of variable is commonly achieved by the so-called control variable transform (CVT),

δx t = B 1 2 δz t , (3.22)
and leads to a new Hessian matrix H = B 1 2 HB 1 2 possessing a lower condition number [START_REF] S A Haben | Conditioning and preconditioning of the variational data assimilation problem[END_REF][START_REF] Yang | Study of variational ensemble methods for image assimilation[END_REF].

Thus, we obtain a modified cost function by injecting the CVT (3.22) into the objective function (3.17) such as

J(δz 0 ) = 1 2 δz 0 2 B -1 + 1 2 t f t 0 ∂ x H∂ x ϕ t (x 0 )B 1 2 δz 0 -D(t, x) 2 R -1 dt, (3.23)

Low rank approximation of the background error covariance matrix

The En4DVar assimilation is designed on such preconditioned incremental variational system and introduces an empirical approximation of the background covariance matrix B. This low rank approximation of the background covariance matrix is directly inspired from the Ensemble Kalman filter where the covariance terms are estimated from the spread of an ensemble of samples.

Denoting < f (t) >= 1 N N i=1 f (i) (t)
as an empirical ensemble mean of a quantity f (t) through N samples, the unbiased empirical background covariance matrix reads

B ≈ 1 N -1 N i=1 x (i),b -< x b > x (i),b -< x b > T . (3.24) Noting A b := 1 √ N -1 (x (1),b 0 -< x b 0 >, • • • , x (N ),b 0 -< x b 0 >
), the perturbation matrix gathering the N zero mean centered background ensemble members as a lowdimensional approximation of the background matrix. Introducing the background covariance approximation in the preconditioned cost function (3.23), we get

J(δz 0 ) = 1 2 δz 0 2 B -1 + 1 2 t f t 0 ∂ x H∂ x ϕ t (x 0 )A b δz 0 -D(t, x) 2 R -1 dt, (3.25)
we define the square root of the background error covariance matrix as the perturbation matrix :

B 1 2 ≈ A b := 1 √ N -1 (x (1),b 0 -< x b 0 >, • • • , (x (N ),b 0 -< x b 0 >)).
Let us note B1 2 := ∂ x ϕ t (x 0 )A b the propagation of ensemble perturbation matrix A b which corresponds to a forecast by the dynamical model of the centered squareroot background covariance matrix. As we rely here on an empirical description of this matrix from a set of samples, we can observe that integrating these samples in time provides us immediately an empirical expression of a low-rank approximation of the background covariance trajectory and of its square root. This avoids thus the employment of the adjoint operator.

The gradient of the cost function is now given by

∇J(δz 0 ) = δz 0 + t f t 0 ( B1 2 T (∂ x H) * R -1 (∂ x H B1 2 δz 0 -D(t, x))dt. (3.26)
and its Hessian reads

H = I + t f t 0 B1 2 T (∂ x H) * R -1 ∂ x H B1 2 dt. (3.27)
Once the minimizer δz 0 estimated, the analysis reads

x a 0 = x b 0 + A b δz 0 . (3.28)
Let us emphasize that, as the empirical covariance matrix B is at most of rank N -1, the control variable has at most N -1 non null components in the eigenspace. Compared to the full 4DVar approach, the control variable's degrees of freedom are thus considerably lowered and the minimization computational complexity is significantly decreased. Indeed, this ensemble version has a lower computation cost if the ensemble forecasting step is distributed on a grid computing.

Preconditioning matrix update

As mentioned earlier, in this study we focus on situations where the background state is only poorly known. It is hence essential to allow in the estimation process a substantial deviation from the background state. So unlike typical incremental ensemble-based variational methods which keep a fixed background covariance and apply a single outer loop of the Gauss-Newton minimization, we propose to update the approximation of this associated error covariance between two consecutive outer loops. The update of the error covariance can be either derived from the ensemble of analysis based on perturbed observations or by a direct transformation of the background ensemble perturbations. The first method relies on a perturbed ensemble of observations, generated with an additional noise with the same standard deviation than in (3.16) :

Y j = Y + j , j = 1, ..., N.
(3.29)

At the kth outer loop iteration, the innovation vector of the member of the initial ensemble x (j),k 0 is defined as,

D (j),k (t, x) = Y j (t, x) -H(ϕ t (x (j),k 0
)), j = 1, ..., N, (3.30) with x

(j),0 0 = x (j),b 0
. Thus a parallel realization of minimization with regard to each member of initial ensemble is conducted

δz (j),k 0 = (A k b ) -1 δx (j),k 0 , j = 1, ..., N. (3.31) J(δz (j),k 0 ) = 1 2 δz (j),k 0 2 + 1 2 t f t 0 ∂ x H∂ x ϕ t (x 0 )A k b δz (j),k 0 -D (j),k (t, x) 2 R -1 dt. (3.32)
Finally the updated initial ensemble field and its perturbation matrix read :

x (j),k+1 0 = x (j),k 0 + A k b δz 0 (j),k
, j = 1, ..., N.

(3.33)

A (k+1) b = 1 √ N -1 x (1),k+1 0 -< x k+1 0 >, ..., x (N ),k+1 0 -< x k+1 0 > . (3.34)
The direct transformation approach corresponds to a linear transformation of the initial error's x b -x t covariance. This approach can take many forms as the transformation matrix is not unique, here we opt for a mean preserving transformation as used in Ensemble Transform Kalman filter, the updated background ensemble perturbations reads,

A (k+1) b = A k b I + t f t 0 B1 2 T t (∂ x H) * R -1 ∂ x H B1 2 dt 1 2 ,k V. (3.35)
It corresponds to the Hessian square root computed from previous perturbation matrix at outer loop iteration k. As the minimization algorithm LBFGS relies on an approximation of the inverse Hessian matrix H -1 , we can use this byproduct to evaluate (3.35). At the initial time, the background matrix is fixed from the initial random conditions chosen. The arbitrary orthogonal matrix V is used to center the posterior ensemble on the updated initial condition/analysis. In this approach a single minimization process is conducted with respect to the background state in opposition to previous cases where the minimization has to be done with respect to each member of the ensemble plus the background state. Finally the updated initial ensemble fields are,

x (j),k+1 0 = xk 0 + √ N -1A (k+1) b , j = 1, ..., N, (3.36)
where xk 0 corresponds to the updated initial state at the kth outer loop,

xk 0 = x k 0 + A k b -1 δz 0 k . (3.37)
Both variants of the update will be assessed in the experimental § 3.4.2.

Localization issues

The previous ensemble method relies on a low rank approximation of the background matrix. This empirical approximation built from only very few samples, compared to the state space dimension, leads in practice to spurious correlations between distant points. For ensemble Kalman techniques, it is customary to remove these long distance correlations through the so-called localization procedure. There are generally two methods to filter the pseudo-correlations.

The first approach introduces a Schur element-wise product between the background correlation matrix and a local isotropic correlation function

P b = C B.
The spatial correlation function can be simply defined as a matrix C( x -y /L) in which we set C xy = 0 when the distance between x and y exceeds the cutoff distance L. Polynomial approximations of a Gaussian function with compact support and a hard cutoff are often employed to that end [START_REF] Gaspari | Construction of correlation functions in two and three dimensions[END_REF]. They lead to sparse correlation matrices, which is computationally advantageous. In order to incorporate the localized background error matrix into our system, we approximate the square root of P b by a spectral decomposition of the isotropic correlation function and keep only the r first leading Fourier modes

C = E n×r λ 1/2 r×r .
(3.38)

The modified perturbation matrix is then provided by

P b = (diag(A (1) b )C, ..., diag(A (N ) b )C). (3.39)
Here the diag operator sets the vector x b,k as the diagonal of a matrix. This localized perturbation matrix is used to precondition the assimilation system associated with (3.25). Remark that this approach is incompatible with the direct transform update of the background error covariance matrix. This is due to the inconsistency of matrix dimensions when updating the background ensemble based on (3.36). As the dimension of P b is n × N × r instead of n × N , the ensemble perturbation matrix cannot be recovered from its localized counterpart.

A variant of covariance localization approach based on Lorenc ( 2003) is implemented in [START_REF] Desroziers | 4denvar : link with 4d state formulation of variational assimilation and different possible implementations[END_REF] along with other localization schemes. [START_REF] Desroziers | 4denvar : link with 4d state formulation of variational assimilation and different possible implementations[END_REF] also highlights the importance of the dimension of the control vector which is directly related to the cost of the minimization algorithm.

Another localization technique proposed by [START_REF] Hunt | Four-dimensional ensemble kalman filtering[END_REF] employs local ensemble. This approach involves a transformation M l from state space R n to local space R l , the local vector is defined as :

x l = M l x n .
(3.40)

Then the analysis process is done in local space R n around each grid point only incorporating the model points and observations within a certain range. This certain range, denoted as l, which corresponds to the concept of cut-off distance aforementioned, determines the size of the local space. This localization strategy is ideally compatible with the method of direct transformation approach associated with the update of background error covariance matrix. The great advantage of this combination is its low computational cost when implemented with properly parallelized minimization procedures. However, good performances can only be reached with a small local space.

All these elements (i.e. CVT, localization and incremental/background error covariance matrix update) associated with a LBFGS minimization strategy constitute the proposed ensemble method. The algorithm descriptions of the overall methods are presented in Figs. 2 and3. We point out that this assimilation system -composed of perturbed observations, one outer loop and localization via a modified covariance matrix -is equivalent to the 4DEnVar method proposed by [START_REF] Liu | An ensemble-based four-dimensional variational data assimilation scheme. part i : Technical formulation and preliminary test[END_REF] and [START_REF] Liu | An ensemble-based four-dimensional variational data assimilation scheme. part ii : Observing system simulation experiments with advanced research wrf (arw)[END_REF]. If direct ensemble transformation update is used, it is hence close to the 4D-LETKF method [START_REF] Fertig | A comparative study of 4d-var and a 4d ensemble kalman filter : Perfect model simulations with lorenz-96[END_REF] and the IEnKS method [START_REF] Bocquet | An iterative ensemble kalman smoother[END_REF] but with a minimization performed on a variational basis.

Construction of the SWE discrete tangent and adjoint

The construction of the discrete adjoint corresponding to the SWE model described in algorithm 8 is the cornerstone of the 4DVar implementation. While the choice of constructing the discrete adjoint over the analytic adjoint in this application can be debated, our goal is to explore the implementation of a variational technique in preparation to the 3D reconstruction (see § 4.2). In this application, we need to construct the most accurate as possible adjoint represenetation of the discrete SWE dynamics (i.e. the finite volume scheme with the Roe approximation).

z 0,l = (A b,l ) -1 δx 0,l
The cost function minimization is done only once in terms of the background. The updated background is obtained by (3.37)

Update the initial ensemble and ensemble perturbation matrix by (3.36) and (3.35)

Reconstruct the analysis and the ensemble form local space to state space end loop Evolve the analysis state x a (t 0 ) to the beginning of the next cycle through the nonlinear dynamics (3.14). The forecast state and the forecast enesemble are used to initialize the next assimilation cycle.

To that end, we build the discrete adjoint procedure associated to the SWE code with the help of the AD tool TAPENADE (Hascoët, 2004). We note that we had to perform a conversion of the SWE code from MATLAB to FORTRAN 90 in order to implement the AD tool TAPENADE. The reader can find the description of the AD principles and the different implementation strategies in appendix C.1. We also provide in appendix C.1.2 an example of the detailed construction of the tangent and adjoint procedure corresponding to a simple program.

In terms of the notations, we inherit from the formalism imposed by TAPENADE such that :

-the tangent procedures are ended by the subscript _d, the differential tangent variables associated to a given variable x are denoted ẋ, -the adjoint procedures are ended by the subscript _b, the differential adjoint variables associated to a given variable x are denoted x. We describe in § 3.3.1, § 3.3.2 and § 3.3.3 the construction of the adjoint procedure of the boundary conditions, the Roe solver and the time discretization respectively. We deduce the discrete adjoint of the SWE code in § 3.3.4. The tangent and adjoint procedures of each subroutine and the entire code are systematically verified by the validation procedures presented in § 3.3.5.

Adjoint of the boundary conditions

The imposed reflecting boundary conditions consist in adding ghosts cells in every directions, and defining the value of non-conservative variable state according to figure 3.3. In the perspective of the computation of the fluxes F E , F W , G N and G S in the east, west, north and south directions with respect to a given cell C ij respectively, the BC routine generates the auxiliary variables x E , x W , x N and x S as described in the left side of figure 3.4. The values at the ghost cells for each variables x E , x W , x N and x S are determined according to figure 3.3.

The adjoint routine BC_b corresponding to the BC routine outputs the adjoint variable x from each adjoint xE , xW , xN and xS which are the output of the adjoint flux\_b routine (see § 3.3.2). The adjoint variable x is the sum of the values at the ghost cells with the values at the boundaries of the adjoint variable state x for each adjoint xE , xW , xN and xS . The reader can find in appendix D.1 the complete adjoint routine of the boundary conditions.

Adjoint of the flux terms

Adjoint of the Roe solver

As described in § 3.1.2, we solve a 1D Riemann problem with a Roe solver to determine the flux crossing every inter-cell in a given direction (see § 3.1.2). In practice, the Roe solver computes the flux crossing from the left to the right cell from the input x L and x R . The values of x L and x R depend on the direction of the flux we aim at computing. We assume beforehand that we computed the quantities x W , x E , x N and x S by the BC routine. The horizontal fluxes F W := F i-1/2,j and F E := F i+1/2,j require the couple (x W , x) and (x, x E ) respectively. In the same way, the vertical flux F S := G i,j-1/2 and F N := G i,j+1/2 require the couple (x S , x) and (x, x N ) respectively.

At first, the Roe solver calls the Priestley routine which computes the parameter variables w L and w R from x L and x R . The parameters w L and w R are then used to compute the eigenvalues, eigenvectors and the coefficient of a linearized Jacobian Ã(x R , x L ). These quantities are computed by the eigenvalue, eigenvector and coefficient routines respectively. Finally, the latter are combined if the nonlinear flux F (x L ) and F (x R ) or G(x L ) and G(x R ) to deduce flux of interest in the given direction. We recall that the nonlinear flux of a given state x in the horizontal and vertical direction is given by F

(x) =   hu hu 2 + 1 2 gh 2 huv   and G(x) =   hv huv hv 2 + 1 2 gh 2  
respectively. The latter is discretized in the algorithm fluxH and fluxV respectively. The adjoint procedure corresponding to the auxiliary routines Priestley, eigenvalue, eigenvector, coefficient, fluxH and fluxV are systematically built with the TA-PENDADE tool and are detailed in appendix § D.2. We assume that we determined these discrete adjoint and deduce the discrete adjoint Roe_b associated to Roe in algorithm 11.

Adjoint of the FVM

The FVM routine combines the fluxes terms F E , F W , G S and G N in each direction according to the FVM formulation (3.6). Each flux term was computed beforehand by the Roe solver described previously. The original FVM routine is described in

Algorithm 11 Roe_b( F ,h L ,u L ,v L ,h R ,u R ,v R ,f lag) ! Forward step =================> ! Preprocessing w L = Priestley (h L , u L , v L ,f lag) w R = Priestley (h R , u R , v R ,f lag) ∆w = w L -w R w = 1 2 (w L + w R ) ! Characteristics of the Roe matrix Ã(X R , X L ) [e 1 , e 2 , e 3 ] = eigenvalues( w,f lag) [r 1 , r 2 , r 3 ] = eigenvectors( w,f lag) [α 1 , α 2 , α 3 ] = coefficients( w,∆w,f lag) ! Backward step <================= ! Compute the left and right fluxes FL = FL + 1 2 F FR = FR + 1 2 F if f lag then [ hL , ūL , vL ] = fluxH_b ( FL ,h L , u L , v L ) [ hR , ūR , vR ] = fluxH_b ( FR ,h R , u R , v R ) else [ hL , ūL , vL ] = fluxV_b ( FL ,h L , u L , v L ) [ hR , ūR , vR ] = fluxV_b ( FR ,h R , u R , v R ) end if [ w, ∆w] = coefficients_b(ᾱ1, ᾱ2 , ᾱ3 , w,∆w,f lag) w = eigenvectors_b(r 1 , r2 , r3 , w) w = eigenvalues_b(ē 1 , ē2 , ē3 , w) wL = wL + 1 2 w wR = wR + 1 2 w wL = wL + ∆w wR = wR -∆w [ hR , ūR , vR ] = Priestley_b ( wR ,h R , u R , v R ) [ hL , ūL , vL ] = Priestley_b ( wL ,h L , u L , v L ) return hL , ūL , vL , hR , ūR , vR algorithm 12. Algorithm 12 FVM(h,q x ,q y ) call q2u(h,u,v,q x ,q y ) call BC(h, u, v, h W , u W , v W , h E , u E , v E , h N , u N , v N , h S , u S , v S ) F E = Roe(h, u, v, h E , u E , v E , f lag) F W = Roe(h, u, v, h W , u W , v W , f lag) G S = Roe(h, u, v, h S , u S , v S , f lag) G N = Roe(h, u, v, h N , u N , v N , f lag) F = ∆t ∆x (F W -F E ) -∆t ∆y (G S -G N ) return F
The q2u routine transforms the conservative variable x c = (h, q x , q y ) = (h, hu, hv) to the non-conservative variable x = (h, hu, hv). The adjoint procedure of the latter is described in appendix D.3. Given the adjoint procedure of the boundary condition and the Roe solver discussed previously, we deduce the discrete adjoint FVM_b described in algorithm 13.

Algorithm 13 FVM_b( F ) ! Forward step =================> call q2u(h,u,v,q x ,q y ) call BC(h, u, v, h W , u W , v W , h E , u E , v E , h N , u N , v N , h S , u S , v S ) ! Backward step <================= FW = ∆t ∆x F , FE = -∆t ∆x F ḠS = ∆t ∆y F , ḠN = -∆t ∆y F [ h, ū, v, hN , ūN , vN ] = Roe_b( FN ,1) [ hS , ūS , vS , h, ū, v] = Roe_b( FS ,1) [ hW , ūW , vW , h, ū, v] = Roe_b( FW ,0) [ h, ū, v, hE , ūE , vE ] = Roe_b( FE ,0)
call BC_b( h, ū, v, hW , ūW , vW , hE , ūE , vE , hN , ūN , vN , hS , ūS , vS ) call q2u_b(h, h,u,ū,v,v,q x ,q x ,q y ,q y ) return h, ū, v

Adjoint of the time integration scheme

The time integration is simply carried out by the RK3 method recalled below,

k 1 = f (x n c ), (3.41) k 2 = f (x n c + 0.5k 1 ), (3.42) k 3 = f (x n c + 2k 1 -k 2 ), (3.43) x n+1 c = x n c + ∆t 6 (k 1 + k 2 + k 3 ), (3.44)
where x n c is the conservative vector x n c = (h n , q n x , q n y ) = (h n , h n u n , h n v n ). The function f calculates the fluxes in each direction as described in § 3.1.2, in other words, f (x n c ) = FVM(h n , q n x , q n y ). [START_REF] Hager | Runge-kutta methods in optimal control and the transformed adjoint system[END_REF] provided the general formulation of the Runge-Kutta methods discrete adjoint procedure and [START_REF] Sandu | On the properties of runge-kutta discrete adjoints[END_REF] showed that the Runge-Kutta discrete adjoint method has the same order of accuracy as the original forward method. Hence, we apply below the general formulation of the discrete adjoint to the RK3 method :

k 1 = f (x n c ), (3.45) k 2 = f (x n c + 0.5k 1 ), (3.46) l 1 = ∂ x f * (x n ), (3.47) l 2 = ∂ x f * (x n + 0.5k 1 ), (3.48) l 3 = ∂ x f * (x n + 2k 1 -k 2 ), (3.49) xn-1 = xn + ∆t 6 (l 1 + 4l 2 + l 3 ), (3.50)
where ∂ x f * is the analytic adjoint operator corresponding to f . In our context, this operator corresponds to the adjoint routine FVM_b such that f (x n c ) = FVM( hn , qn x , qn y ). We deduce the adjoint routine RK3_b summarized in algorithm 14.

Algorithm 14 RK3_b(h n-1 , hn-1 ,q n-1 x ,q n-1 x ,q n-1 y ,q n-1 y ,h n , hn ,q n x ,q n x ,q n y ,q n y ) ! Forward procedure =================> x n = (h n ,q n x ,q n y ) k 1 = FVM(h n ,q n x ,q n y ) tmp = x n + 0.5 * k1 k 2 = FVM(tmp) ! Backward procedure =================> xn = ( hn ,q n x ,q n y ) l 1 =FVM_b(x n ) l 2 =FVM_b(x n + 0.5 * k1) l 3 =FVM_b(x n -k1 + 2 * k1) hn-1 = hn + 1 6 (l 1 (1) + 4 * l 2 (1) + l 3 (1)) qn-1 x = qn x + 1 6 (l 1 (2) + 4 * l 2 (2) + l 3 (2)) qn-1 y = qn y + 1 6 (l 1 (3) + 4 * l 2 (3) + l 3 (3))

Adjoint of the main code

Given the previous adjoint routines, we deduce the adjoint procedure corresponding to the main SWE routine in algorithm 15. As we work with a rather small sized problem, we choose to implement a store-all strategy which consists in saving every quantities h n , u n , v n for each time step t n and invoke them when needed.

Code validation

We validated previous tangent and adjoint procedures by applying the validation tests presented in § 2.4.7. We recall that the original SWE code studied in this chapter

Algorithm 15 SWE_b

Initialization: n = N , hN = 0, ūN = 0 ,v N = 0 while n > 0 do call q2u_b(h n , hn ,u n ,ū n ,v n ,v n ,q n x ,q n x ,q n y ,q n y ) call RK3_b(h n-1 , hn-1 ,q n-1 x ,q n-1 x ,q n-1 y ,q n-1 y ,h n , hn ,q n x ,q n x ,q n y ,q n y ) call u2q_b(h n-1 , hn-1 ,u n-1 ,ū n-1 ,v n-1 ,v n-1 ,q n-1 x ,q n-1 x ,q n-1 y ,q n-1 y ) n = n -1 end while return h0 , ū0 , v0 was implemented in MATLAB, therefore we performed beforehand the conversion to an equivalent FORTRAN code to use the AD tool TAPENADE. The obtained tangent and adjoint procedures are then converted back to MATLAB procedure. Consequently, the following test were performed on the procedures converted to MATLAB.

Tangent process validation

We recall the expression of the tangent validation procedure for the operator f of interest, which represents in this case the original code SWE, and its corresponding tangent denoted ∇f ,

R 1 := f (x 0 + αd) -f (x 0 ) α < ∇f (x 0 ), d > . (3.51)
We evaluate the value of 1-R 1 for different differentiation step α over a temporal window of 10∆t then 100∆t . The results are summarized in table 3.1. The tangent code accuracy meets our expectation as the ratio R 1 is close to 1 up to 10 -α . We obtain the same accuracy for both temporal windows. 10 -4 2.561910 -6 4.116810 -6 10 -5 4.397310 -7 3.567710 -7

10 -6 4.445610 -8 3.537610 -8 10 -7 4.260110 -9 4.368910 -9

10 -8 4.305410 -10 8.092410 -10 10 -9 1.140810 -9 1.369510 -9

Table 3.1 -Evaluation of 1 -R 1 results for different values of the differentiation step α over 10∆t and 200∆t . These results are calculated in MATLAB's default accuracy.

Dot product validation

By definition, the adjoint procedure must correspond exactly to the tangent procedure associated to original code. Therefore, this validation test is carried out in double precision O(10 -16 ) and we expect the results to reach the same accuracy. We computed the ratio R 2 over 200∆t obtained a satisfying result (3.52). Hence, the discrete adjoint corresponds accurately to the previous tangent procedure.

R 2 = < (∂ X M)X, (∂ X M)X > < X, (∂ X M # )Y > = < Ẏ , Ẏ > < Ẋ, Ȳ > = 1.000000000000220. (3.52)

Adjoint process validation

Finally, we evaluated the validity of the entire 4DVar assimilation implementation by computing the ratio R 3 for several differential step α,

R 3 = J(x 0 + αd) -J(x 0 ) α < λ(t 0 ), d > . (3.53)
We recall that we don't require the latter to be close to 1 as the first order finite difference provides a lowly accurate estimate of the cost function. The results are summarized in table 3.2.

Application to the reconstruction of a free surface characteristics

In this chapter, we seek for the complete characteristics (height and velocity field) of the free surface of a fluid contained in a rectangular flat bottom tank of size L x × L y = 250mm × 100mm. Since we have a numerical code which describes the evolution at the free surface in time (see § 3.1.1) and a set of observations, we apply and compare the following data assimilation strategies :

-An incremental 4DVar technique, based on Courtier et al. (1994), which requires the construction of the discrete adjoint code described in § 3.3. We fixed the background error co-variance matrix as a static diagonal matrix B = σ 2 b I for both experiments. σ b was optimally tuned as the standard deviation between the true solution and the a priori experimental initial state.

-Several hybrid data assimilation strategies based on the 4DEnVar developed in [START_REF] Liu | An ensemble-based four-dimensional variational data assimilation scheme. part i : Technical formulation and preliminary test[END_REF], labeled by the suffix "Liu-et-al" in the following figures. The strategies with several outer loops, perturbed observation or direct transformation are indicated by suffix "OL", "OP" and "DT" respectively. Covariance localization and local ensemble are indicated by suffix "LC", "LE" accordingly.

In terms of background error covariance matrix, it is crucial that the initial ensemble represents correctly the background errors, therefore the latter are specified in each case. The reader can find complete study of these strategies in Yang (2014).

Synthetic experiment

General context

In this purely synthetic experiment, we define a reference initial condition as a smooth slope tilted along the x-axis and the y-axis by 21% and 10% respectively. The initial reference velocity field was fixed as a Gaussian field with a standard deviation of 1mm/s. The data assimilation techniques are run on a meshgrid of n x × n y = 101 × 41 nodes. The background is set as a smooth 20% slope along the x-axis only. The slight slope in the y direction induced in the reference height contributes to the divergence between the reference and the background trajectory. The initial reference and background configurations are illustrate in figure 3.5.

Figure 3.5 -Representation of the initial configuration of the reference (on the left) and the background (on the right). The reference is obtained by a smooth 21% slope along the x-axis and a 10% slope along the y-axis. The background height is set as a smooth slope at 20% slope along the x-axis.

In terms of the initialization of the ensemble members, we set up two kind of perturbations of the initial condition. In the one hand, we opt for a Gaussian perturbation of the background. On the other hand, we opt for a parameter perturbation strategy which consists in generating heights with random smooth slope.

The synthetic observations are generated by adding i.i.d Gaussian noise perturbations to the reference free surface height and velocity fields at each grid point. In this experiment, we observe either the free surface height only, or the free surface velocity field only or both at the same time. The observations are acquired every ∆t obs = 50∆t•U/L x of the dynamic step ∆t = 1e -3 s. The common temporal setup of both experiments is illustrated in figure 3.6. The height and velocity observation errors are specified as 0.25mm and 1mm/s respectively. In this experiment, we define the characteristic height ∆h = 50mm as the difference between the maximum and minimum surface height of the initial background state. We also define the characteristic velocity U = √ g∆h = 0.34 m/s as an approximation of the wave phase velocity. For every data assimilation techniques, we performed a single assimilation window of duration t f U/L x = 0.8951. All the following results are adimensionalized with respect to these characteristic values.

We analyze the following results by computing the Root Mean Square Error (RMSE) defined for the given state of interest x with respect to the reference state x ref . We assess display the evolution in time of the RMSE for the background, observation and the resulting analysis states with respect to the values of the reference state. The RMSE function is defined as

RM SE(x) = 1 N N i=1 (x -x ref ) 2 , (3.54)
where N represents the total size of the state x. We also consider the RMSE comparison on a semi-logarithmic graph. In general, a lower RMSE at the end of the assimilation window indicates better performance in terms of the forecast, and a lower temporal RMSE indicates better performance in terms of reconstruction.

Height observations

In this configuration, we only take into account noisy free surface height observations. The En4DVar techniques provided comparable results with the 4DVar for an ensemble of 8 members. Figure 3.7 shows that in terms of the RMSE comparison, the 4DEnVar yields to the best reconstruction of each variables than the standard 4DVar assimilation techniques. Both methods performed strong corrections of the height whereas the standard 4DVar seems to have difficulties to reconstruct the non-observed velocity fields.

Figure 3.7 -RMSE comparison between an incremental 4DVar and 4DEnVar assimilation approaches with a a partial height observation through a noisy free surface height.

Concerning the 4DEnVar methods, the parameter perturbation strategy (4DEnVar-Para) provided the best reconstruction of the complete state. In order to obtain comparable results with the Gaussian perturbation strategy, a localization is necessary and the cutoff distance was set as the optimal cutoff distance minimizing the average RMSE (the same value was used to define the size of local space). The distance increases as the ensemble member increases. We also observe that several outer loops clearly improve the results (OL-LC against Liu-et-al and OL-LE against LE), which highlights the pertinence of the background covariance update. The method with local ensemble (LE) provided slightly better results than the method with localized covariance (LC).

Velocity observations

This time, we only take into account noisy free surface velocity observations. Figure 3.8 reveals a similar behavior than the previous case as the 4DEnVar techniques globally provided a better reconstruction than the standard 4DVar. As in the previous case, the parameter perturbation strategy (4DEnVar-Para) provides the best reconstruction of the complete state for a low ensemble member N = 8. On the other hand, the Gaussian perturbation strategy required a moderately higher ensemble number N = 16 to obtained better comparable results to the 4DVar. The increase of the ensemble size to N = 32 slightly improves the quality of the results at the expense of a higher computational cost.

Height and velocity observations

Finally, we performed the reconstruction with complete observations of the system of interest. The comparison of the each data assimilation techniques provided a very similar behavior than in the previous case. Indeed, the ensemble techniques provided a better reconstruction than the 4DVar in terms of the RMSE evolution. However the advantage of the ensemble methods towards the 4DVar is less enhanced than in the previous case. As we possess the complete observations of the system, we assume that each variable component is mainly corrected by the corresponding observations than indirectly inferred from other observed components.

In terms of the ensemble methods comparison, we obtain globally the same conclusions as previously. The parameter perturbation strategy provided a slightly better result than the Gaussian perturbation. The localization ensemble (LE) approach generally led to better results than the localization covariance (LC) approach, however the more we increase the ensemble number N the more the results of both approaches are alike.

As the data assimilation techniques are comparable in terms of RMSE, we compare their computational needs in table 3.3. The parameter perturbation approach (4DEnVar-PP), requires the less CPU time among all the data assimilation techniques compared. The 4DVar approach requires the most CPU time as it performs several forward and backward time integrations, however in the upside, this technique requires a smaller memory resources. Concerning the Gaussian perturbation, for an equal number of ensemble member, the local covariance approach costs more in time and in memory than the local ensemble approach.

Data assimilation method

CPU time Memory demands Table 3.3 -Comparison of the computation resources required for the classic 4DVar method (4DVar), the parameter perturbation approach (4DEnVar-PP), Gaussian perturbation with a local covariance approach for N = 32 ensemble members (4DVar-LC) and the Gaussian perturbation with a local ensemble approach for N = 32 ensemble members.

Conclusion

The ensemble based data assimilation techniques showed great potential in handling incomplete and noisy observations and each 4DEnVAr strategies clearly outperformed the 4DVar reconstruction in terms of RMSE. Regarding the ensemble initialization, a physical noise, in our case the generation of smooth random slopes, is preferred over non-physical noise which requires a localization procedure. However, in practice we often don't have a profound knowledge of both the dynamical model and the observations to generate a good parameter sampling.

When we consider a set of complete observations, the ensemble-based methods provide similar results with respect to the 4DVar technique. Table 3.3 showed that the 4DVar technique requires the highest CPU time for a relatively small memory cost. Although, these results assume that the En4DVar methods are parallelized beforehand which can represent a tedious task more complex dynamical systems.

Real image data experiment

Experimental observations

The experimental setup simply consists in measuring the height of a colored liquid contained in a tank with a Kinect depth sensor placed above the tank as illustrated in figure 3.9. The measurements were acquired on a 222×88 grid. [START_REF] Combes | Free-surface flows from kinect : Feasibility and limits[END_REF] described in details the measurement technique employed to extract the surface height from the raw data captured by the Kinect sensor. The free surface height observations are characterized by a high level of noise and exhibit large regions of missing data mostly along the boarders which are due to light reflections on the tank's wall. We performed a sequence of procedures that transforms the sequence of spatially sparse height observations into complete height observations defined on 248×98 nodes. However, the computational limitations leads us to run the data assimilation on 124×49 nodes instead. Therefore, the observations are interpolated on the coarse grid. The entire pre-procedure that transforms the height observations acquired on the 222 × 88 grid into a format suitable to run the assimilation techniques is summarized in algorithm 16.

Algorithm 16 Real observations pre-processing 1: Eliminate singular points whose value exceeds a given threshold 2: Fill the missing boundary areas located besides the long borders with the average of the all the pixels' values located in the same section 3: Extrapolate the height profile to the missing boundary areas located besides the short boarders 4: Fill the missing inner hole areas with pseudo-values computed as the mean of all the adjacent pixel's observation values

Similarly to the synthetic configuration, we consider an adimensionalized system with respect to the characteristic height, velocity and time defined as ∆h = 12mm, U = √ g∆h = 0.34 m/s and t f •U/L x = 0.4608 respectively. We estimate the observation error within an unobserved region with respect to the distance from the closest observed point. Thus, the longer the distance, the larger the error. The observation error is however bounded by a threshold of 60%∆h.

The observation time step was tuned cautiously to ∆t obs = 30∆t•U/L x such that the observation sequence is synchronized with the dynamics model.

Choice of the background

In this framework, the initial background was built from the partial observations provided by the Kinect measurements. The initial height is deduced by filtering the observation at the initial time which was obtained by the pre-processing described in algorithm 16. As we notice that the observed free surface behaved roughly as an unidrectional wave along the x-axis, we set the initial x-axis velocity field as a smooth linear slope. The velocity at the top of the wave was set as 23% of the wave velocity and the velocity at the bottom of the wave was set to 0. The transverse velocity field was set to v = 0. The initial observation and background height and velocity are illustrated in figure 3.10. In terms of ensemble sampling, we initialize the ensemble members by adding isotropic Gaussian perturbation fields with standard deviation σ b = 3.6%∆h to the background state. The cutoff distance and the size of local space are fixed as 15%L x .

The assimilation scheme was adapted to sliding assimilation windows to avoid long range temporal correlations. Each window contains five observations and we adopted five windows over nine observations, as illustrated in figure 3.11.

Figure 3.11 -Temporal setup for the sliding data assimilation techniques : we perform the assimilation over a period of t f •U/L x . The observations are injected each ∆t obs = 30∆t•U/L x . We then slide the assimilation window to the next observation and perform once again the data assimilation. This procedure is repeated until we reach the ninth observation.

Results

In opposition to the synthetic case, we don't possess a reference state in this case. We can only assume the true solution to lie within an interval around the observation. Therefore, instead of studying the RMSE with respect to the observations, we compare directly the height fields.

Figure 3.12 shows the evolution of the average surface height of the wave crest of the background, the observations and the result of the different data assimilation techniques. As the observations indicate that the free surface behaves globally as a single wave, we particularly focus on the wave crest's region rather than hte other flat regions. While the background diverges from the observations, we observe that the 4DVar and the 4DEnVar can both follow the observation trajectory tendency. The 4DVar tends to understimate the surface height at the beginning of the simulation whereas the ensemble-based methods follow the trend imposed by the 1st to the 4th observation. After the fifth observation, the En4DVar [START_REF] Liu | An ensemble-based four-dimensional variational data assimilation scheme. part i : Technical formulation and preliminary test[END_REF] diverges from the observation trajectory.

We also compared the free surface in figure 3.13. According to these free surfaces, we can see that the 4DVar solution showed some difficulties to handle the discontinuities at the boundaries of the regions in which the data have been extrapolated. On the other hand, the En4DVar provided a much smoother solution on the borders which corresponded clearly to a better compromise between the observation and the model.

Conclusion and perspectives

In this chapter, we evaluated the capacity of a classical incremental method to reconstruct the complete characteristics (height and velocity field) at the free surface of a unidirectional wave from partial and noisy observations. We had the opportunity to compare the results with the ensemble-based method 4DEnVar, proposed by [START_REF] Liu | An ensemble-based four-dimensional variational data assimilation scheme. part ii : Observing system simulation experiments with advanced research wrf (arw)[END_REF], and various strategies based on the latter [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF].

In general, both methods retrieved successfully the complete characteristics of the free surface in the synthetic and the experimental case. In opposition to the varitional approach, the major advantage of the ensemble-base methods lies in the fact that we avoid the tedious task of constructing the adjoint procedure. In general, the ensemble-based methods provided better results in terms of RMSE for partial and complete observations. In return, these techniques require a significant number of ensemble members to have comparable results with the 4DVar technique. A parametric perturbation of the initial ensemble clearly improves the assimilation process and requires the less computational ressources. However, we often don't have enough knowledge on the system of interest, especially in industrial applications, to generate an adequate physical perturbation. The parallelization of the ensemble members is also required to reach a computational comparable to the 4DVar. While the parallelization can be simply implemented for simple dynamical models, this task gets tougher when we consider a DNS or even a LES model, as the latter often require a parallelization by domain decomposition.

Finally, we note that the incremental 4DVar technique used throughout this chapter could be enhanced by applying a preconditioning technique similar to § 3. throughout this chapter, by injecting the dependence of each variables, height and velocity, from the spatial correlation.

On the other, as the experimental setup is simple to reproduce, we could also add a sensor that measures the free surface velocity field and lead a comparative study as in the synthetic case § 3.4.1. Despite these implementations are out of the scope of this present thesis, the latter can pave the way to significant improvements of the results.

In the following chapter, we chose to remain in variational framework as the gain provided by the En4DVar techniques are relatively small with respect to the implementation difficulties.

Chapitre 4 Reconstruction of a turbulent 3D cylinder wake flow

The reconstruction of three dimensional fundamental turbulent flows is a cornerstone for fluid mechanics applications. As discussed in § 1.3.2, several papers hinted towards the use of data assimilation techniques which combine a numerical model with a set of observations. In this perspective, we aim at building and validating an incremental 4DVar code which is capable to reconstruct a turbulent flow at Re = 300 by combining a DNS code with 2D or 3D observations.

The CFD field provides a wide range of models and numerical schemes to simulate a three dimensional turbulent flow by solving the Navier-Stokes equations (see § 1.1). In the present application, we seek for an accurate representation of the flow, therefore we employed the parallel DNS code Incompact3d developed by [START_REF] Laizet | High-order compact schemes for incompressible flows : A simple and efficient method with quasi-spectral accuracy[END_REF]. In § 4.1, we describe the structure and present briefly the numerical schemes employed in Incompact3d. We describe in § 4.2 the construction and the numerical validation of the adjoint procedure necessary to implement the assimilation technique. We apply and assess the 4DVar code on a ideal experiment in § 4.3 with a set of synthetic 3D observations. Finally, we investigate the possibility to reconstruct a 3D cylinder wake flow with orthogonal-plane stereo PIV observations in § 4.4.

DNS code Incompact3D

We consider the Navier Stokes equations, recalled below, as the dynamics model to describe the evolution of the incompressible fluids of interest, ∇.u = 0, (4.1)

∂ t u + 1 2 (∇(u ⊗ u) + (u∇)u) -ν∇ 2 u + ∇p = 0. (4.2)
In this work, we use the parallel DNS code Incompact3D to solve these equations [START_REF] Laizet | High-order compact schemes for incompressible flows : A simple and efficient method with quasi-spectral accuracy[END_REF][START_REF] Laizet | A numerical strategy to combine highorder schemes, complex geometry and parallel computing for high resolution dns of fractal generated turbulence[END_REF]. Incompact3d is a highly accurate code which combines high-order compact schemes with a spectral solver. It was initially designed for serial processors, then converted to vector processors and more recently converted to parallel platforms [START_REF] Laizet | A numerical strategy to combine highorder schemes, complex geometry and parallel computing for high resolution dns of fractal generated turbulence[END_REF]. The parallel computation, which is achieved by using a 2D decomposition library decomp2d [START_REF] Li | 2decompfft a highly scalable 2d decomposition library and fft interface[END_REF], is a very important asset especially for fluid mechanics applications which usually involve large dimension variables. This section provides a succinct description of the main procedures implemented in Incompact3d. First of all, we present briefly the numerical aspects such as the 2D decomposition library (see § 4.1.1) and the required computational resources (see § 4.1.2). We then present the discretization schemes of the time advancement, the convection and diffusion terms and the pressure treatment in § 4.1.3, § 4.1.4 and § 4.1.5 respectively.

2d decomposition

The 2D decomposition library decomp2d [START_REF] Li | 2decompfft a highly scalable 2d decomposition library and fft interface[END_REF] decomposes the computation domain (n x × n y × n z ) into "pencils" along a given direction x, y or z, respectively referred to as X-pencil, Y-pencil and Z-pencil, as illustrated by figure 4.1. The communication between each processor is taken in charge by "transpose" subroutines provided by 2Ddecomp library. Each processor performs the same instructions bloc to the input data. For instance, let's assume we start from the X-pencil configuration and we aim at computing a gradient ∇f = ∂f ∂x , ∂f ∂y , ∂f ∂z ≈ f i , f j , f k . We compute at first the 1D derivative in the x-axis, f i , then we swap to the Y-Pencil to compute 1D derivative in the y-axis, f j , and finally swap to the Z-Pencil to compute the last 1D derivative in the z-axis, f k . This procedure is extended to every instruction of Incompact3d with the necessary swap operations. For the sake of clarity, we voluntarily omit to mention the swap operations in the following sections.

Computational resources

In terms of computational resources, every calculation of this thesis was run on the computer grid IGRIDA which offers a 125 computer nodes (1200 cores). In terms of storage, IGRIDA offers to each member of the research teams at IRISA/INRIA a shared scratch space of 1.7 TB for temporarily storing simulations inputs and outputs. There are also 5 GPU nodes with NVIDIA GPUs which hasn't been used during this thesis. The computation resources provided by each cluster used in the

Time advancement

The general formulation of the time advancement of the equation (4.2), which deduces the state u k+1 at time t k+1 from the previous state u k at time t k , can be expressed as

u * -u k ∆t = a k F k + b k F k-1 + c k F k-2 -g k ∇p k , (4.3) u * * -u * ∆t = g k ∇p k , (4.4) u k+1 -u * * ∆t = -g k ∇p k+1 , (4.5)
with

F k = - 1 2 ∇(u k ⊗ u k ) + (u k .∇)u k + ν∇ 2 u k (4.6)
and the time averaged values on a given sub-step g k ∆t :

pk+1 = 1 g k ∆t t k+1 t k p(t)dt (4.7)
We specify the discretization scheme by fixing the coefficient triplet

{a k , b k , c k } k=1,••• ,n k (where g k = a k + b k + c k
for simplicity) on n k sub-time steps with t 1 = t n and t n k = t n+1 (∆t = t n+1 -t n being the full time step). The time discretization schemes available in Incompact3d are the Runge Kutta methods at order 3 (RK3) and 4 (RK4), and the Adams-Bashford methods at order 2 (AB2) and 3 (AB3). The values attributed to each coefficient for each technique are summarized in appendix B. In this thesis, we kept the default numerical scheme AB3 for which n k = 1 and (a 1 ,b 1 ,c 1 ) = (23/12, -16/12, 5/12).

Spatial discretization of convective and diffusive terms

The convection and diffusion terms of the Navier-Stokes equations (4.6) are directly solved in the computational domain L x × L y × L z discretized on a regular three-dimensional Cartesian mesh of n x ×n y ×n z nodes. Despite this meshing choice, the code is capable to treat complex geometries via an Immersed Boundary Method [START_REF] Parnaudeau | Combination of the immersed boundary method with compact schemes for dns of flows in complex geometry[END_REF]. The discretization is carried out by a sixth-order, and a second-order on the boundaries, compact finite difference scheme. The following description is based on [START_REF] Laizet | A numerical strategy to combine highorder schemes, complex geometry and parallel computing for high resolution dns of fractal generated turbulence[END_REF].

Let us consider a uniform 1D meshgrid consisting of the distribution of n x nodes {x i } i=1,••• ,nx defined on the domain [0; L x ]. A given derivative f (x i ) at the node x i is approximated by its numerical derivative f i ≈ f (x i ) as follow :

αf i-1 + f i + αf i+1 = a f i+1 -f i-1 2∆x + b f i+2 -f i-2 4∆x , (4.8)
where the coefficients are given by α = 1/3, a = 14/9 and b = 1/9. The second order derivative is therefore deduced by :

αf i-1 + f i + αf i+1 = a f i+1 -2f -i + f i-1 2∆x 2 + b f i+2 -2f i + f i-2 4∆x 2 + c f i+3 -2f i + f i-3 9∆x 2 ,
(4.9) where the coefficients are given by α = 2/11, a = 12/11, b = 3/11 and c = 0. The latter is a sixth-order like approximation with the same favorable properties than the first derivative in terms of spectral-like resolution.

In practice, Incompact3d offers the possibility to impose a combination of three types of boundary conditions : periodic, free-slip and open boundary conditions.

Periodic condition conserves the order of the numerical scheme and is easy to implement. The values of f, f , f at the boundaries are simply deduced by :

f 0 = f nx , f -1 = f nx-1 , f 0 = f nx , f 0 = f nx .
Free slip condition has the same assets than the periodic condition and corresponds whether to the symmetric condition

f 0 = f 2 , f -1 = f 3 , f 0 = f 2 , f 0 = f 2 ,
or the asymmetric condition

f 0 = -f 2 , f -1 = -f 3 , f 0 = f 2 , f 0 = -f 2 .
Open boundary condition doesn't make any assumptions on the flow outside the computational domain. Single sided formulations are used for the approximation of first derivatives and second derivatives for these types of boundaries using relations of the form :

f 1 + 2f 2 = 1 2∆x (-5f 1 + 4f 2 + f 3 ), f 1 + 11f 2 = 1 2∆x 2 (13f 1 -27f 2 + 15f 3 -f 4 ),
that are third order accurate [START_REF] S K Lele | Compact finite difference schemes with spectral-like resolution[END_REF]. At the adjacent nodes, because a three point formulation must be used. Padé schemes are employed with forth-order accurate schemes :

1 4 f 1 + f 2 + 1 4 f 3 = 3 2 f 3 -f 1 2∆x , 1 10 f 1 + f 2 + 1 10 f 3 = 6 5 f 3 -2f 2 + f 1 ∆x .
The open boundary conditions are thereafter specified at the 2nd step of the time advancement (4.4). In our application, as we aim at modeling the downstream of a cylinder wake without knowing the obstacle geometry. Therefore, the inlet is assumed to be known beforehand and is imposed at each time step ∆t of the DNS. The outlet is modeled by a first-order advection model in the streamwise direction that reads :

u BC nx = u * nx -g k ∆x ∆t u * nx (u * nx -u * nx-1 ) (4.10)

Pressure treatment

The tedious incompressibility condition (4.1) can be verified at the end of each sub-time step ∇•u k+1 = 0, (4.11) through solving a Poisson equation

∇ 2 pk+1 = ∇•u * * g k ∆t , (4.12)
that provides the estimation of the pressure pk+1 from the velocity u * * (obtained by (4.4)) required to perform the correction (4.5).

The inversion of a 3D Poisson equation is in itself a tedious task which requires sophisticated and expensive methods such as in [START_REF] Mercier | A multidimensional compact higher-order scheme for 3-d poisson's equation[END_REF] who use multidimensional compact high order schemes. On the other hand, the Poisson equation can be easily solved in Fourier space by the means of a Fast Fourier Transform (FFT) library. [START_REF] Laizet | A numerical strategy to combine highorder schemes, complex geometry and parallel computing for high resolution dns of fractal generated turbulence[END_REF] provides a review of the parallel FFT libraries to be implemented in Incompact3d. Initial attempts were carried out to adapt existing FFT packages, however each library has their limitations and aren't always practical in terms of implementation with the decomposition library. [START_REF] Li | 2decompfft a highly scalable 2d decomposition library and fft interface[END_REF] created from scratch a generic FFT based on the derivation of Glassmans's general N Fast Fourier Transform [START_REF] Ferguson | A simple derivation of glassman general-n fast fourier transform[END_REF]) and avalaible at http://www.jjj.de/fft/glassman-fft.f. The latter is written in FORTRAN 77 and provides an interface which supports complex to complex, real to complex/complex to real transforms.

In order to reduce any spurious oscillations on the pressure field, the pressure discretization is performed on a regular staggered mesh grid shifted by a half-mesh in each spatial direction as illustrated in figure figure 4.2. The partial staggering, initially proposed by [START_REF] Wilhelmson | Direct solutions for poisson's equation in three dimensions[END_REF] and further developed by [START_REF] Canuto | Spectral methods in fluid dynamics[END_REF] and [START_REF] Chen | Acceleration in turbulent channel flow[END_REF] using second-order schemes, is easy to implement because it does not require extensive use of mid-point interpolations which are computationally expensive when combined with high-order schemes.

The evaluation of the first order derivative f i+1/2 at the staggered nodes is given by a sixth-order finite difference scheme expressed by : (4.13) with α = 9/62, a = 63/62 and b = 17/62. The spectral behaviour of this scheme is better than its collocated counterpart (4.8), however we have to consider in addition the combination with a mid-point intepolation procedure which is given by :

αf i-1/2 + f i+1/2 + αf i+3/2 = a f i+1 -f -i 2∆x 2 + b f i+2 -f i-1 3∆x 2 ,
αf i-1/2 + f i+1/2 + αf i+3/2 = a f i+1 + f i 2 + b f i+2 + f i-1 2 , (4.14)
which provides a sixth-order estimation of f i+1/2 for α = 3/10, a = 3/4 and b = 1/20 [START_REF] Nagarajan | A robust high-order compact method for large eddy simulation[END_REF].

Initialization from a previous state

In preparation for the 4DVar implementation, we are led to restart a DNS from a pre-existing turbulent flow and the inflow computed or observed at each DNS time step. In order to ensure that the restart procedure after a given time t k doesn't induce any errors and the flow doesn't diverge from its initial trajectory, Incompact3d generates the restart file sauve.dat att k which is saves the following elements :

-each component of the velocity field at t k , -the the pressure field p k at t k , -the gradient of the pressure field ∇p k at t k , -and the previous convection and diffusion terms F k-1 , F k-2 necessary to the time advancement scheme (4.3).

Construction of Incompact3d discrete adjoint

In this chapter we implement the 4DVar technique described in § 2.4.5. The construction of the discrete adjoint corresponding to the Incompact3d code is the cornerstone of the 4DVar implementation. In this work, the discrete adjoint code was built with the help of the AD tool TAPENADE (Hascoët, 2004). The reader can find the description of the AD principles and the different implementation strategies in § C.1. We also provide in § C.1.2 an example of the detailed construction of the tangent and adjoint procedure corresponding to a simple program.

In the present chapter, we consider the complete DNS code Incompact3d which is structured into a sequence of subroutines as illustrated in figure 4.3. In terms of the notations, we inherit from the formalism imposed by TAPENADE such that :

-the tangent procedures are ended by the subscript _d, the differential tangent variables associated to a given variable x are denoted ẋ, -the adjoint procedures are ended by the subscript _b, the differential adjoint variables associated to a given variable x are denoted x. First and foremost, we present in § 4.2.1 the treatment of the parallelization routines involved in the code. We then describe in § 4.2.2, § 4.2.3, § 4.2.4, § 4.1.5 and § 4.2.6 the discrete adjoint procedures corresponding to each subroutine of the original program depicted in figure 4.3. Most of these procedures involve the first and/or second order derivatives which are defined by high order compact schemes. The latter are a algebraic combination of coefficients and values of the flow at given points, therefore we deduce their corresponding discrete adjoint by a systematic application of TAPENADE. We refer the reader to appendix E. Finally, the complete discrete adjoint of Incompact3d is summarized in § 4.2.7. The tangent and adjoint procedures of each subroutine and the entire code was systematically verified by the validation procedures presented in § 4.2.8.

Discrete adjoint of parallel-MPI codes

Incompact3d and in particular the decomp2d library are based on MPI subroutines which need to be differentiated in the same way as the other subroutines of the code. The numerical gradient of MPI-codes have been investigated in the literature by [START_REF] Faure | Extension of odyssée to the mpi library-reverse mode[END_REF], Utke et al. (2009) and [START_REF] Schanen | Interpretative adjoints for numerical simulation codes using mpi[END_REF]. We base the following description on Utke et al. ( 2012).

In the one hand, we consider the global communicators which carry out operations that involve all the data within all the processors. For instance, the global communicator MPI_allreduce(x,y,...,SUM,P,...) sums the values of x from each processor into the global variable y in the processor P. The procedure adjoint corresponding to the global communicators are summarized in figure 4.4, where P denotes a specific processor to which the result is sent to for MPI_bcast, MPI_reduce, MPI_gather and MPI_scatter.

Original Adjoint

MPI_bcast(x, ... ,P, ...) MPI_reduce (x, t, ... , SUM, P, ...) x = 0.0 ; on P : x = x + t

MPI_reduce(x,y, ... ,SUM, P, ...) on P : t = ȳ ; on P : ȳ = 0.0 MPI_bcast ( t, ... , P, ...) ; x = x + t

MPI_allreduce(x,y, ... ,SUM, ...) MPI_allreduce (ȳ, t, ... , SUM, ...) ȳ = 0.0 ; x = x + t

MPI_gather(x,y, ... ,P, ...) MPI_scatter (ȳ, ... , t, P, ...) on P : ȳ = 0.0 ; x = x + t

MPI_scatter(y, ... ,t, P, ...) MPI_gather (ȳ, ... , t( :), P, ...) on ȳ = 0.0 ; on P : x = x + t(:) On the other hand, the point-to-point communicators consists in sending or receiving the data between two given processes. The reader can find in Utke et al. (2009) the adjoining of point-to-point communication. In our context, these communicators are exclusively used within the decomp2d library to swap from a pencil to another [START_REF] Li | 2decompfft a highly scalable 2d decomposition library and fft interface[END_REF]. It is shown that there exists a duality between forward and adjoint MPI communication routines for several applications [START_REF] Faure | Extension of odyssée to the mpi library-reverse mode[END_REF][START_REF] Cheng | A duality between forward and adjoint mpi communication routines[END_REF]. This duality can be extended in our context with the swap of a configuration to another. In other words, the adjoint of the swap procedure transpose_x_to_y which swaps from the X-pencil to the Y-pencil is given by transpose_y_to_x which swaps from the Y-pencil to the X-pencil. By analogy, the adjoint of the swap procedure transpose_y_to_z is given by transpose_z_to_y.

For the sake of clarity, we avoid to mention the swap procedures in the following sections and we assume that each operation performs the adequate forward and backward swap before each procedure.

Discrete adjoint of the convection and diffusion terms

Discrete adjoint of the convective term

The convective term, namely 1 2 (∇(u ⊗ u) + (u∇)u), requires the knowledge of the first order derivatives of the velocity field u, denoted derx(u), dery(u) and derz(u). Given these derivatives, the computation of the convective term, summarized into the vector conv = (conv x , conv y , conv z ), with respect to the velocity field u = (u, v, w), is given by algorithm 17.

Algorithm 17 conv(u,v,w,conv x ,conv y ,conv z ) t x = u.u t y = u.v t z = u.w conv x = 0.5(derx(t x )+dery(t y )+derz(t z )) + 0.5(u.derx(u) + v.dery(u) + w.derz(u)) t x = v.u t y = v.v t z = v.
w conv y = 0.5(derx(t x )+dery(t y )+derz(t z )) + 0.5(u.derx(v) + v.dery(v) + w.derz(v)) t x = w.u t y = w.v t z = w.w conv z = 0.5(derx(t x )+dery(t y )+derz(t z )) + 0.5(u.derx(w) + v.dery(w) + w.derz(w))

We detail the construction of the adjoint for the first component conv x the other components conv y and conv z are then deduced by analogy.

The tangent of the intermediate terms t x , t y , t z is straightforward and reads

  ṫx ṫy ṫz   =   2u 0 0 v u 0 w 0 u     u v ẇ   . (4.15)
The adjoint is simply deduced by transposing the previous matrix such as

  ū v w   =   2u v w 0 u 0 0 0 u     tx ty tz   . (4.16)
We assume that we know beforehand the tangent procedures derx_d, dery_d and derz_d, respectively corresponding to each derivative subroutines derx, dery and derz, and we denote their associating Jacobian matrix as D x , D y and D z . Thus, given (4.15), we deduce the tangent of ċonv :

ċonv x = 1 2 D x D y D z   2u 0 0 v u 0 w 0 u     u v ẇ   + 1 2 derx(u) dery(v) derz(w)   u v ẇ   + 1 2 u.D x v.D y w.D z   u v ẇ   .
Given the transposed of the Jacobian matrices D T x , D T y and D T z and (4.16), we deduce the expression of the adjoint variable conv

x   ū v w   = 1 2   2u v w 0 u 0 0 0 u     D T x D T y D T z   conv x + 1 2   derx(u) dery(v) derz(w)   conv x + 1 2   u.D T x v.D T y w.D T z   conv x .
The discrete adjoint of each convective term conv y and conv y can easily be deduced by analogy. The discrete adjoint of the convection procedure is finally presented in algorithm 18.

Discrete adjoint of the diffusion term

The diffusion term, namely -ν∇ 2 u, requires the knowledge of the second order derivatives of the velocity field u, determined by the subroutines derxx(u), deryy(u) and derzz(u) respectively. The computation of the diffusion vector diff = (dif f x , dif f y , dif f z ) is given by algorithm 19.

Similarly to the previous case, we assume that we know beforehand the tangent and discrete adjoint associated to the second order derivative subroutine. Thus the adjoint of the the diffusive terms are given by algorithm 20. Let us note although this operator is formally auto-adjoint, the discretization of its continous adjoint doesn't respect this property. Thus the discrete adjoint has to be specified by AD.

Discrete adjoint of the convection and diffusion terms

Given the previous routines, the convection and diffusion terms are summed into the vector F := (F x , F y , F z ) = conv + diff. We easily deduce the discrete adjoint convdiff_b summarized in algorithm 22.

Algorithm 18 conv_b (u,ū,v,v,w, w,conv x ,conv x ,conv y ,conv y ,conv z ,conv z ) Initialization: ū = 0, v = 0, w = 0 ! Forward procedure call conv(u,v,w,conv x ,conv y ,conv z ) ! Background procedure tmp1 = derx_b(conv x , conv x ) tmp2 = dery_b(conv y , conv x ) tmp3 = derz_b(conv z , conv x ) ū = ū + 0.5(2u.tmp1 + v.tmp2 + w.tmp3) + 0.5 derx(u).conv x + 0.5.u.tmp1 v = v + 0.5u.tmp2 + 0.5 dery(v).conv x + 0.5.v.tmp2 w = w + 0.5u.tmp3 + 0.5 derz(w).conv x + 0.5.w.tmp3 tmp4 = derx_b(conv x , conv y ) tmp5 = dery_b(conv y , conv y ) tmp6 = derz_b(conv z , conv y ) ū = ū + 0.5v.tmp4 + 0.5 derx(u).conv y + 0.5.u.tmp4 v = v + 0.5(u.tmp4 + 2v.tmp5 + w.tmp6) + 0.5 dery(v).conv y + 0.5.v.tmp5 w = w + 0.5v.tmp6 + 0.5 derz(w).conv y + 0.5.w.tmp6 tmp7 = derx_b(conv x , conv z ) tmp8 = dery_b(conv y , conv z ) tmp9 = derz_b(conv z , conv z ) v = ū + 0.5w.tmp7 + 0.5 derx(u).conv z + 0.5.u.tmp7 w = v + 0.5w.tmp8 + 0.5 dery(v).conv z + 0.5.v.tmp8 ū = ū + 0.5(u.tmp7 + v.tmp8 + 2w.tmp9) + 0.5 derz(w).conv z + 0.5.w.tmp9

Algorithm 19 diff(u,v,w,dif f x ,dif f y ,dif f z ) Initialization: u, v, w dif f x = -ν(derxx(u)+deryy(u)+derzz(u)) dif f y = -ν(derxx(v)+deryy(v)+derzz(v)) dif f z = -ν(derxx(w)+deryy(w)+derzz(w))

Discrete adjoint of the time integration scheme

We recall the general formulation of the time advancement procedure described in § 4.1.3,

u * = u k + a k ∆tF k + b k ∆tF k-1 + c k ∆tF k-2 ,
where F k is the result of the convdiff routine. The time integration scheme can be written into the following matrix form

    u * F k F k-1 F k-2     =     1 a k ∆t b k ∆t c k ∆t 0 1 0 0 0 0 1 0 0 0 0 1         u k F k F k-1 F k-2     .
(4.17)

Algorithm 20 diff_b(u,ū,v,v,w, w,dif f x ,dif f x , dif f y , dif f y ,dif f z ,dif f z ) Initialization: ū = 0, v = 0, w = 0 ! Forward procedure Initialization: u, v, w call diff(u,v,w,dif f x ,dif f y ,dif f z ) ! Background procedure ū = ū + (derxx_b(-ν.dif f x , -ν.dif f x ) ū = ū+deryy_b(-ν.dif f x , -ν.dif f x ) ū = ū+derzz_b(-ν.dif f x , -ν.dif f x ) v = v + (derxx_b(-ν.dif f y , -ν.dif f y ) v = v+deryy_b(-ν.dif f y , -ν.dif f y ) v = v+derzz_b(-ν.dif f y , -ν.dif f y ) w = w + (derxx_b(-ν.dif f z , -ν.dif f z ) w = w+deryy_b(-ν.dif f z , -ν.dif f z ) w = w+derzz_b(-ν.dif f z , -ν.dif f z ) Algorithm 21 convdiff(u,v,w,F x ,F y ,F z ) call conv(u,v,w,conv x ,conv y ,conv z ) call diff(u,v,w,dif f x ,dif f y ,dif f z ) F x = conv x + dif f x F y = conv y + dif f y F z = conv z + dif f z Algorithm 22 convdiff_b(u,ū,v,v,w, w,F x , Fx , F y , Fy ,F z , Fz ) Initialization: ū = 0, v = 0, w = 0 ! Forward procedure call conv(u,v,w,conv x ,conv y ,conv z ) call diff(u,v,w,dif f x ,dif f y ,dif f z ) ! Backward procedure conv x = Fx ; conv y = Fy ; conv z = Fz dif f x = Fx ; dif f y = Fy ; dif f z = Fz call diff_b(u,ū,v,v,w, w,dif f x ,dif f x , dif f y , dif f y ,dif f z ,dif f z ) tmp1 = ū ; tmp2 = v ; tmp3 = w call conv_b(u,ū,v,v,w, w,conv x ,conv x ,conv y ,conv y ,conv z ,conv z ) ū = ū + tmp1 ; v = v + tmp2 ; w = w + tmp3
The corresponding adjoint procedure is then simply deduced by transposing the latter matrix such that

    ūk F k F k-1 F k-2     =     1 0 0 0 a k ∆t 1 0 0 b k ∆t 0 1 0 c k ∆t 0 0 1         ū * F k F k-1 F k-2     . (4.18)
Therefore, the generalized form of the discrete adjoint of the time integration scheme is given by algorithm 24.

Algorithm 23 intt_b(u k ,ū k , F k , F k-1 , F k-2 ) Initialization: F k = 0, F k-1 = 0 and F k-2 = 0 F k = F k + a k ∆tū k F k-1 = F k-1 + b k ∆tū k F k-2 = F k-2 + c k ∆tū k 4.2.

Discrete adjoint of the velocity divergence

At this point, we proceed to several operations before computing the pressure field. First and foremost, we impose the boundary conditions on each component of the velocity field and in each directions. A pre-correction treatment is applied to ensure the proper application of the boundary conditions. Finally, the velocity divergence ∇•u required to solved the Poisson equation is developped (see § 4.2.5). All these operations are gathered under the designation of velocity divergence in the structure of the code (see figure 4.3). As the velocity divergence ∇•u is computed by an elementary combination of the derivative routines, we omit here its description and refer the reader to appendix E. Thus, the present section describes the discrete adjoint of the boundary conditions and the pre-correction routine.

Boundary conditions

In terms of boundary conditions, we constructed the discrete adjoint corresponding specifically to the cylinder wake flow configuration studied within this chapter. In this work, we impose the following boundary conditions model :

-periodic condition along the z-axis, -free slip condition along the y-axis, -Dirichlet condition for the inflow and -a first order advection model for the outflow.

As discussed in § 4.1, the periodic and free-slip conditions are implicitly imposed by the numerical schemes used in the first and second order derivatives subroutines. As the inflow is a control parameter of the 4DVar assimilation technique, we refer the reader to § 2.4.5 for its adjoint. Finally, we recall the first order transport model in the stream-wise direction employed to model the outlet boundary condition,

u BC nx = u * nx -c x (u * nx -u * nx-1 ), (4.19)
where u * is the result of the time integration scheme (4.3) and c x = g k ∆t ∆x u * nx is set such as the CFL condition is verified.

The tangent of outflow model with respect to the differential variables u * can be written into the following matrix form

uBC nx uBC nx-1 = 1 -c x c x u * nx -u * nx-1 0 1 0   u * nx u * nx-1 ċx   .
The corresponding adjoint procedure is obtained by transposing the previous matrix :

  ū * nx ū * nx-1 cx   =   1 -c x 0 c x 1 u * nx -u * nx-1 0   ūBC nx ūBC nx-1
.

The adjoint procedure of the outflow boundary is given by algorithm 24.

Algorithm 24 outflow_b(ū * , v * , w * ,ū BC , vBC , wBC )

Initialization: ūnx = 0, ū * nx-1 = 0, v * nx = 0, v * nx-1 = 0, w * nx = 0, w * nx-1 = 0, ! Forward procedure c x = g k ∆t ∆x u nx ! Background procedure ū * nx = ūBC nx + c x ūBC nx v * nx = vBC nx + c x vBC nx w * nx = wBC nx + c x wBC nx cx = cx -(u nx -u nx-1 )ū nx -(v nx -v nx-1 )v nx -(w nx -w nx-1 ) wnx ū * nx-1 = ūBC nx-1 + c x ūBC nx v * nx-1 = vBC nx-1 + c x vBC nx w * nx-1 = wBC nx-1 + c x wBC nx ū * nx = ūBC nx + g k ∆t ∆x cx

Pre-correction

We perform the pre-correction on the boundaries for each component of the velocity field (4.4) recalled below,

u * * = u * + g k ∆t∇p k ,
where u * and ∇p k is the velocity field and the average of the pressure gradient respectively and computed at the previous time step. This pre-correction ensures that the boundary condition model is effectively applied after the correction of the pressure gradient deduced by solving the Poisson equation (see § 4.2.6). The discrete adjoint corresponding to this pre-correction is straightforward and reads :

ū * = ū * + ū * * , ∇p k = ∇p k + g k ∆tū * * .
A second correction on the outflow is performed to ensure that the outflow rate corresponds to the inflow rate. This correction is an elementary operation which requires the use of a global communicator MPI_AllReduce in order to compute the sum of the flow on the entire plane. The adjoint associated to routine is described in § 4.2.1.

Discrete adjoint of the Poisson solver

As described in § 4.1.5, the incompressibility condition is verified by solving a Poisson equation,

∇•∇p k+1 = ∇.u * * g k ∆t ,
in the spectral space. In this work, we choose to construct the discrete adjoint of the generic FFT developed by [START_REF] Li | 2decompfft a highly scalable 2d decomposition library and fft interface[END_REF] which is based on the derivation of Glassmans's general N Fast Fourier Transform [START_REF] Ferguson | A simple derivation of glassman general-n fast fourier transform[END_REF]. The FFT routine denoted SPCFFT can be found in http://www.jjj.de/fft/glassman-fft.f. The construction of the adjoint procedure corresponding to SPCFFT revealed that the routine is autoadjoint in the numerical sense which eases the implementation of the adjoint operator.

The poisson routine takes as input the right hand side of the Poisson equation, denoted rhs, and outputs the pressure field p. The procedure is described by algorithm 25. The wave number kxyz is determined by the auxilary procedure waves, as it doesn't depend on an active variable of the routine, there is no need to compute its adjoint. The routine poisson also depends on two auxiliary routines assemble, disassemble and FFT_postproc. The routines assemble and disassemble are conjugate routines which simply performs a redistribution of the element of the input matrix. In other words, the discrete adjoint of assemble corresponds to disassemble and vice versa. The routine FFT_postproc performs a one dimensional post processing applied in every direction. The reader can find in appendix F the construction of its corresponding discrete adjoint.

Algorithm 25 poisson(rhs,p)

Initialization: call waves(kxyz) assemble(rhs) FFT(rhs,cw,1) ! Post processing in the spectral space call FFT_postproc(cw,cwp,1) ! Solve the poisson equation tmp1 = real(kxyz) ; tmp2 = aimag(kxyz) cwp = cmplx(-cwp/tmp1, -cwp/tmp2) ! Backward post processing in the spectral space call FFT_postproc(cwp,cwpp,-1) FFT(p,cwpp,-1) assemble(p)

Given the discrete adjoint of the auxiliary routines, the discrete adjoint corresponding to the poisson routine is deduced by algorithm 26.

Algorithm 26 poisson_b(rhs,p)

Initialization: call waves(kxyz) disassemble(p) FFT(p,cwpp,-1) ! Post processing in the spectral space call FFT_postproc_b(cwpp,cwp,-1) ! Solve the poisson equation tmp1 = real(kxyz) ; tmp2 = aimag(kxyz) cw1 = cmplx(-cwp/tmp1, -cwp/tmp2) ! Backward post processing in the spectral space call FFT_postproc_b(cwp,cw,1) FFT(rhs,cw,1) assemble(rhs)

Discrete adjoint of the pressure gradient and velocity correction

The pressure gradient ∇p k+1 , computed by the subroutine gradp, is determined by sixth-order finite difference schemes described in § 4.1.5. The numerical schemes are linear combinations of the value of the pressure at different position of the staggered grid. Therefore, the corresponding adjoint procedure gradp_b can be determined straightforwardly by using TAPENADE.

We recall the expression of the final correction (4.5) of the velocity field with the pressure gradient computed previously u k+1 = u * * -g k ∆t∇p k+1 , its discrete adjoint ūk+1 , implemented to the corgp_b routine, is straightforwardly deduced by : ū * * = ū * * + ūk+1

∇p k+1 = ∇p k+1 + -g k ∆t∇ ˜uk+1 .

Discrete adjoint of the main program

All in all, the discrete adjoint of the original program Incompact3d is summarized in algorithm 25.

Code validation

We applied each validation test presented in § 2.4.7 to each adjoint routine built previously. We show in this section the results run on the discrete adjoint Incompact3d_b (see algorithm 27) corresponding to the original Incompact3d code.

Algorithm 27 Incompact3d_b (u 0 ,ū 0 )

Initialization: k = N ,u 0 ,ū N = 0, F N -1 = 0, F N -2 = 0 while k>0 do corgp_b(u k ,ū k ,∇p,∇p) gradp_b(∇p,∇p,p,p) poisson_b(p,p,u * * ,ū * * ) pre_correc_b(p,p) intt_b(u * ,ū * ,F k-1 , F k-1 ,F k-2 , F k-2 ) convdiff_b(u k-1 ,ū k-1 ,F k-1 , F k-1 ,F k-2 , F k-2 ) k = k -1 end while

Tangent process validation

We recall the tangent validation procedure for the operator f , which represents in this case the original code Incompact3D, and its corresponding tangent denoted ∇f ,

R 1 := f (x 0 + αd) -f (x 0 ) α < ∇f (x 0 ), d > . (4.20)
We evaluate the ratio R 1 for different differentiation step α over a temporal window of 10∆t then 100∆t . The results are summarized in table 4.2. The tangent code accuracy meets our expectation as the ratio R 1 is close to 1 up to 10 -α . We obtain the same accuracy for both temporal windows.

α R 1 for 10 iterations R 1 for 100 iterations 10 -3 1.0006899708371482 1.0000100833054026 10 -4 1.0000689567331724 1.0000001285470035 10 -5 1.0000068952890127 1.0000000040572325 10 -6 1.0000006894438878 1.0000000003169462 10 -7 1.0000000695308335 1.0000000000302249 10 -8 1.0000000375275666 0.99999999989057253 10 -9 1.0000000677486891 0.99999999908904569

Dot product validation

By definition, the discrete adjoint procedure must correspond exactly to the previous tangent procedure associated to original code. Therefore, as the validation test is carried out in FORTRAN double precision (O(10 -14 ), we expect the result to reach the same accuracy. We computed the following ratio R 2 for 50∆t and obtained the following result which met our expectations,

R 2 = < (∂ X M)X, (∂ X M)X > < X, (∂ X M # )Y > = < Ẏ , Ẏ > < Ẋ, Ȳ > = 0.99999999999998990. (4.21)

Adjoint process validation

Finally, we evaluated the validity of the entire 4DVar assimilation code by computing the following ratio for several differential step α

R 3 = J(x 0 + αd) -J(x 0 ) α < λ(t 0 ), d > . (4.22)
This test was at first carried out at first in FORTRAN double precision, then in FORTRAN single precision which corresponds to the accuracy of the simulation and assimilation. In opposition to the previous validation test, we don't require the adjoint procedure to be close to the finite difference as the latter is a rough approximation of the cost function gradient. The results are summarized in table 4.3.

Validation on a synthetic case with 3D observations

We consider an ideal synthetic experiment to apprehend the mechanism of the reconstruction of a 3D turbulent flow using the 4DVar code developed during this thesis. This experiment doesn't aim at reproducing any specific experimental setting, instead our goal is to assess the capacity of the code to correct a background flow with a set of 3D observations sparse in time and at the DNS spatial resolution. The experimental context and the synthetic reference flow are described in § 4.3.1. As it is unfeasible to perform the assimilation within the same domain as the reference flow, we investigate the consequence of running the assimilation within a smaller domain in § 4.3.2. The description of the synthetic observations and the initial background are given in § 4.3.3 and § 4.3.4, respectively. Finally, the results are presented in § 4.3.5.

Synthetic configuration

In this synthetic context, we consider a cylinder wake flow at Re = 300 as a reference flow. The latter is generated by running Incompact3d from an initial null velocity field within the reference domain Ω ref = 20D × 12D × 6D descritized into 360 × 217 × 108 nodes. The turbulent structures are induced by imposing an inflow velocity u in /U = 1 + at each DNS time step ∆t = ∆t•U/D = 0.005. The cylinder is modeled by an Immersed Boundary Method [START_REF] Parnaudeau | Combination of the immersed boundary method with compact schemes for dns of flows in complex geometry[END_REF][START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at reynolds number 3900[END_REF]. We impose free-slip conditions on the y-axis borders and periodic conditions along the cylinder axis. The outflow is modeled by a first-order advection model (see § 4.1). We obtain satisfying turbulent structures at T ref = 34000∆t = 170 (see figure 4.5). The temporal setup for the generation of a reference flow followed by the ap-plication of the 4DVar is illustrated in figure 4.6. The characteristics of each run performed within the present synthetic configuration are summarized in table 4.4. We then extract the three components of the velocity at t 1 = 170 and add a 15% noise and a sin(y 3 ) profile. The latter is used to initialize the 4DVar code performed over the assimilation window T a = 50∆t = 0.25, by injecting three dimensional synthetic observations every ∆t obs = 10∆t .

Initialization from a pre-existing turbulent flow

In regards to the computational resources, we clearly can't afford to run the 4DVar assimilation on the reference domain Ω ref . Unfortunately, running the DNS on a smaller domain inevitably leads to incompatible boundary conditions models of the reference domain and the smaller domain. Furthermore, most of the experimental setups can't provide the quantities required for the restart procedure described in § 4.1.6 (pressure field and convection/diffusion terms at previous time steps). We assess in the present section the inherent errors induced by the model and the initialization with the components of the velocity field only.

We ran Incompact3d within the smallest sub-domain possible Ω = 8D ×8D ×6D which contains the turbulent structure, over a period of 500∆t = 2.5. This DNS

Run Lx × Ly × Lz nx × ny × nz Re Duration CPU time (h) Reference 20D × 12D × 6D 361 × 217 × 108 300 30000∆t 50 4DVar 8D × 8D × 6D 145 × 145 × 108 300 50∆t 5 4DVar-rampe 8D × 8D × 6D 145 × 145 × 108 300 50∆t 5
Table 4.4 -Characteristics of each run performed to reconstruct the three dimensional flow from 3D synthetic observations. Each data assimilation cases considered a static background error covariance matrix B -1 = 1. The 4DVar and 4DVar-rampe runs fix the observation error covariance matrix to R -1 = 1 and to the formulation given by (4.23), (4.24) and (4.25), respectively.

is only initialized with the velocity field components directly extracted from the reference flow described in § 4.3.1, the pressure and gradient pressure are set to zero and we use of a first order temporal scheme at the start. We also input at each DNS time step ∆t the exact inflow extracted from the reference. In our context, we recall that we perform the 4DVar over a small assimilation window T a = 50∆t as illustrated in figure 4.6. Within the assimilation window, the errors are reasonably confined at the boarders and remain in a range of [-0.01; 0.01].

A solution could consist in using sponge layers [START_REF] Bogey | Numerical simulation of sound generated by vortex pairing in a mixing layer[END_REF] which will require the construction of the adjoint as suggested by [START_REF] Lemke | Adjoint based pressure determination from piv-data -validation with synthetic piv measurements[END_REF]. However, as explained in § 4.3.3, for the sake of simplicity we chose to reduce boundary zone effects by tuning the observation error covariance matrix R. Finally, we perform the statistics assessments within a smaller Ω obs , illustrated in figure 4.5.

Synthetic observations

In this ideal synthetic configuration, we assume we possess a set of 3D synthetic observations every ∆t obs = 10∆t , the observation time step. These observations are directly extracted from the reference flow within the sub-domain Ω.

In light of the analysis made in § 4.3.2, the errors induced by the incompatibility of the boundary condition models over the assimilation window T a are confined within boundary zones (see figure 4.7). Therefore, we attempt to reduce the effect of these errors by tuning the observation error covariance matrix R within the subdomain Ω. We run two assimilation cases : Case a : we assume that we trust the observations which comes to impose R -1 = I within the sub-domain Ω. This case is referred to as "4Dvar-all" in the following results. Case b : we consider the incompatibilities between the model and the observations within the boundary zones. This case is referred to as "4Dvar-rampe" in the following results. We set the observation covariance matrix as a decreasing function depending on the distance to the boundaries such that

R -1 (x, y, z) = diag exp - (x -6D) 2 2 , x ∈ [6D; 8D], (4.23) R -1 (x, y, z) = R -1 diag (x, y, z)exp - (y -1D) 2 2 , y ∈ [0D; 1D], (4.24) R -1 (x, y, z) = R -1 diag (x, y, z)exp - (y -7D) 2 2 , y ∈ [7D; 8D]. (4.25)
The resulting observation covariance matrix R -1 is illustrated in figure 4.8.

Background choice

In this purely synthetic application, we build a background initial condition which is likely to generate a divergent trajectory to the reference flow. Our main concern in this experiment is to assess the 4DVar code ability to retrieve the reference from an un-phased flow and few observations by correcting the initial state. We recall that we don't aim at reproducing a synthetic equivalent of a realistic experiment in this application.

As illustrated in figure 4.6, the background initial condition is deduced from the reference flow at t 1 by adding a 15% white noise and a sin y D 3 profile bias along the y-axis. The comparison between the synthetic reference and the background initial stream-wise velocity fields is illustrated in figure 4.9. In this configuration, we fix the background error covariance matrix to B = 1.

Results

We assess the reconstruction of the velocity fields and compare each experimental case described in § 4.4.2. Both cases were carried out on a single assimilation window T a with a single outer loop iteration. Since the observation time step is small (∆t obs = 10∆t ), there are few changes between two consecutive observations. Besides, the observed inlet corresponds to the reference, thus the inlet correction didn't bring much improvements to the analysis flow. Therefore, we only assess the correction brought by the assimilation to the initial condition at t 1 .

Optimization characteristics

First and foremost, we compare the evolution of the cost function J(δγ 0 ) and the norm of its gradient ∇J(δγ 0 ) with respect to the number of inner loop iterations. Figure 4.29 shows that both cases have a similar convergence rate of the cost function. In the one hand, the case a performs more inner loop iterations and reach a lower cost function evaluation than in case b. In the other hand, the case b found the correct descent gradients most of the time. We assume that the observation covariance matrix guided the case b in the right path.

RMSE comparison

We study the evolution of the Root Mean Square Error (RMSE) defined for a given state x with respect to a reference state x ref such that,

RM SE(x) = 1 N -1 N i=1 (x -x ref ), (4.26)
where N is the size of the state of interest x. We compute and compare the RMSE for each component of the velocity field within the sub-domain Ω. The RMSE provides a global estimation of the correction of the large scale structures. We compare in figure 4.11 the RMSE evolution of the background, observation and the analysis results in both cases, with respect to the reference flow over the data assimilation window T a = 50∆t . Each RMSE was computed within the subdomain Ω obs . We recall that case a and case b correspond to the labels "4DVar-all" and "4DVar-rampe", respectively. In both cases, the data assimilation performed few corrections of the initial state and the solution tends to get closer to the reference in time. Both cases provide very similar results, we note that the case b provides the closest solution to the reference flow for the components v and w. .11 -Assimilation with 3D synthetic observations : comparison of the RMSE evolution of the background (in blue), the observations (in green), the analysis in case a (in red) and the analysis in case b (in black) with respect to the reference state for each velocity components over the assimilation window of T a = 50∆t within the sub-domain Ω obs .

Map comparison

We compare in figure 4.12 the snapshots of each component u, v and w of the velocity field in the stream-wise plane z = 3D at t 1 = 170. The figure shows that over the short assimilation window T a = 50∆t , the 4DVar technique was able to correct the bias imposed on the initial background flow. In the downside, both cases haven't lowered the level of noise added in the initial background flow despite the observations and the inputted inlet are noise-free. We also compare in figure 4.13 the snapshots of each component u, v and w of the velocity field in the stream-wise plane z = 3D at t 2 = 170.25. We notice that the DNS code is able to fairly denoise the initial perturbed flow in time. However the imposed bias influences the behavior of the background flow. Both assimilation cases corrected the bias and ensured a better reconstruction in time. The case b provides a slightly better result than the case a in terms of noise reduction within the sub-domain Ω obs . This improvement is especially observed on the w component.

Energy spectrum comparison

Finally, we compare in figure 4.14 the energy spectrum in each case at t 1 = 170 and t 2 = 170.25. In both 4Var methods corrected the flow at the largest scale but didn't perform very well at the smallest scale. We notice that case b performed a slightly better correction of the initial smallest scales than case a. At the end of the assimilation window, the analysis results tend to stick to the background flow, nevertheless the case b still shows slight improvement with respect to the case a.

Conclusion and discussion

We assessed throughout this section the ability of the 4DVar code to correct a non-physical noise and bias over a short assimilation window T a = 50∆t .

First and foremost, we identified the inherent source of model errors that can be encountered during the assimilation. In the one hand, these errors stem from the incomplete initialization of the DNS code. While the errors remain in a small range of err/U = [-0.01; 0.01] over 500∆t , they might be much more amplified at a higher Reynolds number. A solution would consist in adding an additional control of the initial quantities in the same spirit as the control on the initial velocity components. On the other hand, the incompatibility between the boundary conditions model led to local errors on the y-axis borders and the outlet. We could have used sponge layers and construct their corresponding adjoint in the same vein as [START_REF] Lemke | Adjoint based pressure determination from piv-data -validation with synthetic piv measurements[END_REF]. However, the implementation of the parallel 4DVar code for two different domains is a delicate task so we opted for an assimilation within the domain Ω and tuned the covariance matrix R and B (see § 4.3.3).

Globally, the 4DVar code was able to reconstruct the initial conditions from a noisy and biased background and a set of ideal three dimensional observations. The tuning of the observation error covariance matrix R improves slightly the solution at the largest and the smallest scale and eases the optimization process.

Towards the reconstruction of a full 3D velocity field with orthogonal-plane stereo PIV observations

We investigated in this section the capacity of the 4DVar technique to reconstruct a fully three dimensional flow with orthogonal-plane stereo PIV observations. In this context, we consider the reconstruction of a cylinder wake at Reynolds 300 given a set of observations within the inlet and the stream-wise observation planes. We layout in § 4.4.1 the experimental context and the measurement techniques employed to provide the experimental observations. We consider at first the reconstruction of the three dimensional velocity field with synthetic observations in § 4.4.2. We then consider in § 4.4.3 the reconstruction of the three dimensional velocity field with real observations extracted from two stereo-PIV systems.

Experimental context

The experimental cylinder wake flow was generated by the wind tunnel located at Irstea regional center in Rennes. The wind tunnel, illustrated in figure 4.15, was designed to simulate experimentally a mixing layer flow and was adapted to simulate experimentally a cylinder wake flow. It is equipped with two independent aeraulic circuits, each of them equipped with a 7.5kW blower, a 42kW cold battery and a 31kW heater. The two parallel air streams are blown separately from two air intakes into a conditioning chamber composed of filters and screens, before they are blown into a pressure conditioning chamber. This configuration can generate low turbulence homogeneous flows. The air streams are then blown into a convergent entrance cone with a contraction ratio of 2.5 and are separated at the inlet by an aluminum plate. Both flows end up in a test section of 3m × 1m × 1m. The velocity and temperature of each incoming flows can be set continuously and independently in a range of 0.5 to 5m/s and in a range of 5 to 35 • C. We acquired alternatively 3000 image pairs of the inlet and the stream-wise flows with SPIV1 and SPIV2 respectively, with a time interval of 2000µs between two successive images (500Hz). The time acquisition of each stereo PIV system is illustrated in figure 4.17. As the wind tunnel was initially designed for a mixing layer configuration, a cushion at the separation plaque was added to reduce the mixing layer effect. The circular cylinder has a length L = 1000mm and a diameter D = 10mm. It is equipped with two thin rectangular end plates with the specification recommended by Stansby (1974). The distance between the end plates is 300mm providing an aspect ratio L/D = 30. The clearance between the walls and the end plates is about of 350mm which is much larger than the thickness of the boundary layer. The blockage ratio is 1. The free stream velocity is adjusted to U = 0.48m/s so that the Reynolds number Re is 300.

Synthetic orthogonal plane stereo PIV observations

General context

We consider the same cylinder wake reference flow defined in § 4.3.1 to generate the synthetic observations. We recall that the reference flow is obtained by running Incompact3d from an initial null velocity field and by inputting an inflow velocity within the reference domain Ω ref = 20D × 12D × 6D . We conserve the same DNS step ∆t = 0.005 and run the simulation until T ref = 30000∆t . We mimic the PIV observations by extracting alternatively the inlet plane x = 0D and a stream-wise plane z = 3D from the reference flow every ∆t obs = 10∆t . Thus, every observation in a given plane occurs every ∆t PIV = 20∆t . The time setup is illustrated in figure 4.18 and the characteristics of each run performed during this experiment are summarized in table 4.5.

Synthetic observations

In order to apply the 4DVar code to this problem, we have to carry out a preprocedure which generates three dimensional observations from each occurrence of Table 4.5 -Assimilation with synthetic orthogonal-plane observations : characteristics of each run performed to reconstruct the three dimensional flow from 3D synthetic observations. The empty box case corresponds to the simulation performed to obtain an acceptable initial background from an initial null velocity field and the injection of an inlet every ∆t PIV = 20∆t . The data assimilation method 4DVar considered a static background error covariance matrix B -1 = 1 and the observation error covariance matrix R -1 built by algorithm 29.

an observed plane. We draw the reader's attention on the necessity to filter the observations at the borders in order to verify the free-slip boundary conditions and ensure that the code won't explode numerically. The entire pre-procedure that transforms the 2D observations into the 3D observations is summarized in algorithm 28.

u(x, y, z) = u 1 (y, z)w 1 (x) + u 2 (x, y)w 2 (z) w 1 (x) + w 2 (z) 6:
Filter the values at the boundary conditions in the y-axis to fit the free-slip boundary conditions

In light of the previous results in the synthetic case (see § 4.3), we tune the observation error covariance matrix with respect to the level of confidence at each point. In the orthogonal-plane observations context, we impose the highest confidence on the observed planes and lower rapidly the confidence in time. The construction of the observation error covariance matrix is summarize in algorithm 29 and is illustrated

Choice of the intial background

As we don't possess any additional information aside the observed planes, the background initial solution is obtained by performing the DNS from a null velocity and input the reference inflow every ∆t . The inlet planes at every ∆t are obtained by a linear temporal interpolation of the observed inlet every ∆t PIV = 20∆t . This procedure is referred to as "Empty box" in table 4.5.

A first simulation was carried out by considering an ideal inflow which is directly Algorithm 29 Construction of the invert observation error covariance matrix R -1

1: We set the most confidence on the observed orthogonal planes x = 0D and z = 3D

R -1 (x, y, z) = diag max exp - x 2 0.2 , exp - (x -3D) 2 0.2 2:
We reduce the confidence along they-axis to take into account the incompatibility between the background free-slip boundary condition and observations' model

R -1 (x, y, z) = R -1 (x, y, z).diag exp - (y -1D) 2 0.2 , y < 1D R -1 (x, y, z) = R -1 (x, y, z).diag exp - (y -7D) 2 0.2 , y > 7D 3:
We reduce the confidence along the x-axis to take into account the incompatibility between the background and observations model at the outlet

R -1 (x, y, z) = R -1 (x, y, z).diag exp - (x -6D) 2 0.2 , x > 6D 4:
We reduce the confidence along the z-axis as the generated observations don't strictly verify the periodicity condition We induced a 10% synthetic noise into the inflow and performed once more the simulation from scratch. The synthetic noise is adjusted according to a Gaussian profile to avoid the apparition of instabilities at the free slip boundaries which can compromise the simulation. After T b = 5000∆t , the obtained background has a different structure than the reference as illustrated in figure 4.22. We also notice that the noise injected at the inlet is denoised by the DNS code.

R -1 (x, y, z) = R -1 (x, y, z).diag exp - (z -1D) 2 0.2 , y < 1D R -1 (x, y, z) = R -1 (x, y, z).diag exp - (z -5D) 2 0.2 , y > 5D y = 4D z = 3D x = 0D

Results

We show in this section the results of the assimilation run with the initial background and synthetic observations described in § 4.4.2 and § 4.4.2, respectively. The assimlation was carried out with two outer loop iterations over the data assimilation window T a = 50∆t.

We analyze the evolution of the cost function J(δγ 0 ) and the norm of its gradient ∇J(δγ 0 ) with respect to the number of inner loop iterations. Figure 4.29 shows that during the first outer loop iteration, the value of the cost function decreased rather slowly. The first correction of the second outer loop provided a better correction and led the assimilation to a better solution.

We display the evolution of the velocity norms of each component of the reference, background, observation and analysis at each observed planes x = 0 and z = 3D. Background Figure 4.24 shows that the assimilation achieved notable corrections at the inlet plane every ∆t PIV , however, in general the analysis trajectory tends to stick to the background trajectory. The snapshots of each velocity components showed very few visible corrections of the background. Therefore, the assimilation failed to provide a satisfying correction of the initial and inflow conditions with synthetic observations, however the assimilation performed better with real observations as described in § 4.4.3.

Real orthogonal-plane stereo PIV observations

General context

The PIV measurements were calibrated to the same resolution than the DNS. In this configuration, the cylinder diameter is discretized on 1D = 17 nodes instead of 18 nodes. The velocity observations are defined on 20D×16D at the inlet plane x = 0 and on 16D × 20D at a stream-wise plane z = 3D. The characteristic velocity is set as U = 0.48m/s. The temporal setup of this experiment is described in figure 4.25.

The characteristics of each run performed within this experiment are summarized in table 4.6. Table 4.6 -Assimilation with real orthogonal-plane observations : characteristics of each run performed to reconstruct the three dimensional flow from 3D piv observations. The empty box case corresponds to the simulation performed to obtain an acceptable initial background from an initial null velocity field and the injection of an inlet every ∆t PIV = 20∆t . The data assimilation method 4DVar considered a static background error covariance matrix B -1 = 1 and the observation error covariance matrix R -1 built by algorithm 29.

Similarly to the synthetic case, we define the smallest sub-domain possible which contains the turbulent structures. In this context, the definition of the sub-domain depends on the observations. The real observations structures aren't as regular as the synthetic structures and are occasionally spread along the y-axis. In order to fit as much as possible the free-slip condition along the y-axis, we have to increase the sub-domain to L y = 10D. In addition to the errors inevitably induced by the intialization of the DNS and the boundary condition model incompatibilities as described in § 4.3.2, the observed flow doesn't present a strict periodicity length as assumed in the DNS. We estimated instead the statistical periodicity length by comparing the mean and the standard deviation of each component of the velocity field in time at different position at fixed position in the inlet plane. Figure 4.26 shows a comparison between the mean and the standard deviation of each component of the velocity field in time at a punctual position. We fixed the periodic length at 6D. All in all, we run the assimilation on the sub-domain Ω = 8D × 10D × 6D.

In terms of observations, we reconstruct the three dimensional observations from the two dimensional observations by applying the algorithm 28. The first step of the algorithm is replaced by a filter that We also employ the algorithm 29 to construct the R -1 matrix which is illustrated in figure 4.27 at plane y = 4D, z = 3D and x = 0D. In this configuration, as we don't have a strictly periodic flow it is necessary to lower the confidence along the z-axis.

Choice of the intial background

Similarly to the synthetic case described in § 4.3.4, the background initial solution is obtained by performing the DNS from a null velocity and input the temporal interpolation of the reference inflow every ∆t . The latter is obtained by a linear temporal interpolation of the observed inlet every ∆t PIV . Figure 4.28 compares the snapshots the real observation with the initial background at the stream-wise plane z = 3D obtained after running 4000∆t from scratch. We notice that the background flow is less detailed with respect to the observations. We also notice that the flow seems to diverge from the observation at the outflow as the influence of the inlet faded and the outflow boundary condition is incompatible with the observation. Nevertheless, the three dimensional background flow is conform to the dynamics of the system and captures the trend of the observations on the observed plane, which is sufficient to start the 4DVar code.

Results

We analyze the evolution of the cost function J(δγ 0 ) with respect to the number of inner loop iterations. Figure 4.29 shows that the assimilation process has difficulties to find the correct descent direction.

As we only observe two orthogonal planes, it is clearly unreasonable to trust the extrapolated data within the 3D volume. Therefore, we limit our statistics comparison to the observed plane x = 0 and z = 3D. Moreover, the observations are noisy and only represent a trend of the true solution which lies somewhere in between the observations and the dynamical representation of the flow. Therefore, instead of studying the evolution of the RMSE in the observed planes, we analyze norm of each velocity component of the observation, the background and the analysis at the observed planes x = 0 and z = 3D. Figure 4.30 shows that the analysis achieves notable corrections on the inlet similarly to the the synthetic case § 4.4.2. This phenomenon could be improved by implementing a temporal average on the inflow condition as proposed in [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF]. At the stream-wise plane, the analysis trajectory generally tends to get closer to the observations, excepted for the v component.

We analyze in figure 4.31 the snapshots of each component u/U , v/U and w/U of the velocity field in the observed stream-wise plane z = 3D at the beginning of the assimilation window t 0 . We compare the snapshots of the observations, the background and the analysis obtained by the 4DVar code. We observe that the analysis retrieved similar turbulent structures to the observations for each component, while it conserved the same velocity levels as the background flow. This correction leads to a denoised and more detailed flow in the observed plane z = 3D.

In terms of inlet correction, the correction of the inlet at the beginning of the assimilation window t 1 = 20 is too small to be discerned directly as illustrated in figure 4.32.

We can assess the correction performed within the entire volume by analyzing each component in the cross-section plane y = 4D. Figure 4.33 shows that the analysis flow leveled the values of each velocity component to follow the trend imposed by the observations. Nevertheless, the analysis map remains close to the background as it doesn't possess any other information beside the observed planes.

Conclusion and discussion

We carried out the reconstruction of the three components of the velocity fields from a sequence of two orthogonal-plane synthetic and real observations. In both cases, it is crucial to perform a post processing procedure on the observations described in algorithm 28. The latter generated by Incompact3d from a null velocity field with a input of the temporal interpolation of the inflow at each ∆t . In the synthetic case, we are able to retrieve a very comparable flow to the reference, however, by adding noise to the inflow we obtain a sufficiently modified flow to justify the use of the variational assimilation technique. In opposition, the initial background obtained by the temporal interpolation of the real inflow observations isn't as detailed as the observations but provides a sufficient estimation of the flow to perform the assimilation.

In the one hand, we performed the assimilation with synthetic noisy observations. In this configuration, the assimilation had difficulties to converge and achieved negligible corrections of the background trajectory. This case could be improved by implementing a temporal average on the inflow condition similarly to [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF]. As the assimilation performed better results with real orthogonalplane observations, this case lays out the difficulty of generating synthetic realistic observations.

On the other hand, we achieved the reconstruction of the three dimensional flow from pre-processed observations acquired from a stereo-PIV setup. This experiment represented a particularly tedious task as the observations don't obey to a strict periodicity conditions imposed at the boundaries along the z-axis in addition to the model incompatibilities described in § 4.3.2. The observation error covariance matrix R played an important role in this reconstruction as it indicates the areas in which the boundary conditions model is incompatible. We also remind that these results were obtained on a short assimilation window T a = 60∆t containing 6 observations including 4 observations at the inlet x = 0. The increase of this temporal window can clearly ease the reconstruction of the flow.

Conclusions and perspective

Conclusion

The fluid mechanics field context possess a wide range of various methods to reconstruct a realistic turbulent flow. During these last decades, the trend is to couple experimental measurements with dynamical models, particularly by using data assimilation techniques. These techniques emerged from the Numerical Weather Prediction field and are well suited for large sized problems.

We covered in chapter 2 the main data assimilation approaches so far. In the framework of this thesis, we applied a variational data assimilation approach, referred to as 4DVar in the following, which consists in solving a optimal control problem under the constraint of the dynamical model. This approach is well-known for its accuracy and relies on the construction and the implementation of a so-called adjoint procedure. The latter is built upon the original dynamics code and provides the gradient descent direction for the optimization. The construction of the adjoint is eased by the use of an automatic differentiation tool but requires a systematic and rigorous validation procedure to ensure that the accuracy of the original dynamics code is conserved.

Reconstruction of a free surface flow

In chapter 3, we implemented and evaluated a classical incremental 4DVar technique [START_REF] Courtier | A strategy for operational implementation of 4d-var, using an incremental approach[END_REF] to a simplified SWE code. We compared this method with the ensemble-based method 4DEnVar, proposed by [START_REF] Liu | An ensemble-based four-dimensional variational data assimilation scheme. part i : Technical formulation and preliminary test[END_REF], and various strategies based on the latter [START_REF] Yang | Enhanced ensemble-based 4dvar scheme for data assimilation[END_REF], to reconstruct the complete characteristics of a free surface from partial and noisy observations.

In general, both methods retrieved successfully the characteristics of the free surface with synthetic and real-ground observations. The major advantage of the presented 4DEnVar strategies lies in the fact that we don't have to construct the adjoint procedure in opposition to the 4DVar approach. In general, the ensemblebased methods provided good results in terms of RMSE for partial and complete observations without the use of the adjoint procedure. In return, these techniques require a significant number of ensemble members to have comparable results with the 4DVar technique. The parallelization of the generation of the ensemble members also yields a computational cost comparable to the 4DVar. While the parallelization can be easily achieved for simple dynamical models, this task can be tedious for a DNS or even a LES model which require themselves a parallelization by domain decomposition.

Reconstruction of a turbulent three dimensional flow

The reconstruction of the three dimensional three components of a turbulent flow was at the heart of this thesis. In chapter 4, we carried out the reconstruction of the three dimensional three components velocity fields at the downstream of a cylinder wake at Reynolds 300. In this chapter, we combined the parallelized highly-accurate DNS code Incompact3d and a set of orthogonal-plane stereo PIV obserations. This study follows the work of [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF] who performed successfully the reconstruction in a two dimensional context.

First and foremost, we emphasize on the difficult task of the adjoint construction. Indeed, we had to reach the double objective of conserving the accuracy and the parallelized structure of the original Incompact3d code, as laid out in § 4.2.

Once we validated the adjoint procedure, we carried out a first reconstruction in an ideal context. In § 4.3, we assessed the ability of the 4DVar code to retrieve a reference flow from an initial noisy and biased flow. We encountered difficulties to find a compromise between the size of the domain (an implicitly the computational cost) and the accuracy of the flow with respect to a reference run in a larger domain. While the DNS code alone is capable to denoise the initial noise, the 4DVar successfully corrected the initial flow bias. The tuning of the observation error covariance matrix also improved the reconstruction and paved the way for future improvements.

Finally, we performed in § 4.4.3 the 4DVar with a set of orthogonal-plane stereo PIV observations. We draw the reader's attention on the difficulty of this challenge. Indeed, despite the quality of the DNS code and the DNS-like resolution of the PIV measurements, the passage from a 2D reconstruction to the 3D reconstruction is a delicate task. We described in § 4.4.3 the construction of three dimensional observations and their associated covariance error matrix.

Perspective

In light of the conclusions of this work, we lay out various subjects for further research.

Preconditioning model

In the first place, the incremental variational method could be enhanced by implementing a preconditioning technique similar to § 3.2.1 to avoid the ill-posed nature of the system of interest. We can also improve the variational method by updating the background error covariance matrix B by injecting the dependence of each controlled variables from the spatial correlation. These lines of research can be easily implemented for simple dynamics and be compared to hybrid data assimilation techniques in the same vein as the work carried out in chapter 3.

Control of the boundary conditions

In chapter 4 we achieved the correction of the flow by controlling the initial condition and the inflow at each time step. In the first step, a short term improvement to the presented results consists in implementing a temporal average on the inflow condition similarly to [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF]. As further study in § 4.3.2 revealed that the boundary conditions and the incomplete initialization of the DNS can lead to additional errors. While these errors remained in a reasonable range, at a higher Reynolds they are likely to be amplified. This situation could be avoided by including the control of the initial pressure field and every boundary conditions of the domain into the assimilation process.

Towards three dimensional observation

Despite the possibility to reconstruct a turbulent flow with orthogonal-plane stereo PIV observations, it is preferable to have a three dimensonal representation of the observation. During the last decade, experimental fluid dynamics community head to 3D PIV techniques. While the tomo-PIV measurements are only capable to retrieve a very small domain (40 mm × 40 mm × 10mm) [START_REF] G E Elsinga | Tomographic particle image velocimetry[END_REF], recent studies are able to retrieve cubic-meter sized volume [START_REF] Kuhn | Large-scale tomographic particle image velocimetry using helium-filled soap bubbles[END_REF][START_REF] Scarano | On the use of helium-filled soap bubbles for large-scale tomographic piv in wind tunnel experiments[END_REF]. These observations are particularly sparse in space and could be coupled with an LES model in the future.

Annexe B Time integration discretization schemes

We recall the general discretization of the time advancement of (4.2) which deduces the state u k+1 at time t k+1 from the previous state u k at time t k ,

u * -u k ∆t = a k F k + b k F k-1 + c k F k-2 -g k ∇p k , u * * -u * ∆t = g k ∇p k , u k+1 -u * * ∆t = -g k ∇p k+1 , . (B.1)
We specify the discretization scheme by fixing the coefficient triplet {a k , b k , c k } k=1,••• ,n k (where g k = a k + b k + c k for simplicity) on n k sub-time steps with t 1 = t n and t n k = t n+1 (∆t = t n+1 -t n being the full time step). The discretization schemes available in Incompact3d are the Runge Kutta methods at order 3 (RK3) and 4 (RK4), and the Adams-Bashford methods at order 2 (AB2) and 3 (AB3). The values attributed to each coefficient for each technique are given by :

AB2 n k = 1 (a 1 ,b 1 ,c 1 ) = (3/2 ,-1/2,0) AB3 n k = 1 (a 1 ,b 1 ,c 1 ) = (23/12 ,-16/12,5/12)
RK3 [START_REF] Williamson | Low-storage runge-kutta schemes[END_REF]

n k = 3 (a 1 ,b 1 ,c 1 ) = (8/15 ,0,0) (a 2 ,b 2 ,c 2 ) = (

Automatic Differentiation

Automatic Differentiation (AD) is a set of techniques that aims at evaluating derivatives of a given function defined by a computer program. Regardless of the complexity of a given computer program, it consists in a sequence of simple instructions consisting in elementary arithmetic operations (additions, multiplications etc...) and/or elementary functions (exponential, cosinus etc...). As we know beforehand the derivative of each instruction, AD consists simply in applying the chain rule of derivative calculus repeatedly. Conceptually, AD is different from symbolic and numerical differentiation and takes advantage of both approaches. Indeed, the application of the chain rule avoids round-off errors encountered by finite differentiation and provides the most accurate derivatives of the original program. At the same time, it avoids the delicate task of performing a discretization by providing directly a discretized program.

We describe in this chapter the main principles of automatic differentiation, with the different computational strategies. We then describe the implementation techniques and present the AD tool used to ease the implementation. Finally, we illustrate the AD technique on a simple example. The reader can find further information in [START_REF] Hascoët | Automatic differentiation by program transformation[END_REF][START_REF] Hascoët | Adjoints by Automatic Differentiation[END_REF].

C.1 Automatic Differentiation principles

Let's consider a continuous function f : R m → R n which can be decomposed by a sequence of k functions such as

f = f p • f p-1 • • • • • f 1 . The value of f at an initial value x := x 0 ∈ R m is given by f (x) = f p (x p-1 ) = f p-1 • f p-1 (x p-2 ) = f p-1 •f p-1 (x p-2 )•• • ••f 1 (x 0 )
. The function f is implemented as a program P composed of a sequence of p instructions such that P = {I 1 ; ... ; I p-1 ; I p }. Our aim consists in defining the differential tangent and adjoint program of P corresponding to the function f with respect to x ∈ R m . According to the chain rule, the differential tangent associated to f with respect to x is given by :

f (x) = f k (x p-1 ) • f p-1 (x p-2 ) • • • f 1 (x 0 ) (C.1)
The derivative f (x) is too large and requires expensive computation resources. We can counter this issue by introducing the differential vector ẋ such that :

ẏ = f (x) • ẋ = f k (x p-1 ) • f p-1 (x p-2 ) • • • f 1 (x 0 ) • ẋ (C.2)
Instead of computing each matrix f k of size n for k ∈ {1, • • • , p} and deduce f (x) • ẋ, we compute each vector of size n × m. This is done by computing (C.2) from right to left : we first compute the vector f

1 (x 0 ) • ẋ, then f 2 (x 1 ) • (f 1 (x 0 ) • ẋ) and finally deduce ẏ = f (x) • ẋ.
The reverse mode is deduced from the chain rule by introducing the differential adjoint vector ȳ ∈ R n and the adjoint derivatives f # p for p ∈ {1, ..., k} such that :

x = f * (x) • ȳ = f * 1 (x 0 ) • • • f * p-1 (x p-2 ) • f * p (x p-1 ) • ȳ (C.3)
As in the tangent mode, (C.3) is also computed from right to left which rises an important issue. Indeed, for the first step calculation for instance, we need the value of x p-1 which, we recall, is obtained by In the framework of this thesis, we are led to opt for the checkpointing SA strategy. This choice stems from the implementation of the AD as discussed in the following section.

x p-1 = f p-1 (x p-2 ) • • • • • f 1 (x 0 ).

C.1.1 Automatic Differentiation implementation

The numerical implementation of the AD can be achieved by applying one of the following strategies :

Source Code Transformation (SCT) consists in replacing the source code by an automatically generated from the source code. This new code includes additional variables, arrays and data structures that hold the derivatives and all the necessary new instructions to compute them. It is also easier to implement the reverse mode with this strategy. On the downside, the transformation of the source code is more complex as the complexity of the source code increases. Therefore, in practice we resort to AD tools that performs automatically and rapidly the transformation. This strategy can also be implemented for all programming languages and the resulting program can be compiled into an efficient code.

Operator Overloading (OO) consists in telling the compiler that each real number is replaced by a pair of real numbers, the second holding the differential. This strategy is obviously possible for a source code written in a language which supports operator overloading. Each elementary operation on real numbers is overloaded, in other words, it is replaced by a new one. working on pairs of reals, that computes the value and its differential. In opposition to source code transformation, the original program is virtually unchanged since everything is done at compile time. However, the resulting program will run slowly because it constantly builds and destroys pairs of real numbers. Moreover, it is very hard to implement the reverse mode with this strategy.

In the framework of this thesis, we turn to the SCT strategy and perform the code transformation with the help of the AD tool TAPENADE (Hascoët, 2004;[START_REF] Hascoët | The Tapenade Automatic Differentiation tool : Principles, Model, and Specification[END_REF]. TAPENADE, which is successor of Odyssée [START_REF] Faure | Odyssée user's guide version 1[END_REF], was developed by the Tropics team (Inria Sophia Antipolis) since 1999. This tool is specialized among others in AD reverse and tangent mode and adopts a store-all strategy with the possibility of checkpointing on calls. This tool is implemented at 90% Java and 10% C and supports code sources written in Fortran (Fortran95, Fortran177 and older). TAPENADE was previously employed in several such as optimum design problems citepHascoet2003,Dervieux2006, oceanography [START_REF] Ferron | Capacités actuelles de la différentiation automatique : l'adjoint d'opa par tapenade[END_REF]Tber et al., 2007) and was also used for the reconstruction of a 2D turbulent flow by a variational assimilation technique [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF]. In the framework of this thesis, we recall that we aim at reconstruction the state of a 2D flow chapter 3 and a 3D turbulent flow chapter 4. On the one hand, the 2D application is based on a homemade MATLAB code, which can be easily translated into a Fortran code ; and in another hand, the 3D application is based on the DNS code Incompact3d which is entirely written in Fortran. Therefore, TAPENADE is a good candidate for both applications and was naturally selected.

In practice, TAPENADE can be used directly as a web server by a distant user, without any need for a local installation at the following address : http://tapenade.inria.fr:8080/tapenade/index.jsp However, it is usually more convenient to download and install the tool on a local computer from the ftp server ftp://ftp-sop.inria.fr/tropics/tapenade and call the tool directly as a command with arguments. The reader can refer to TAPENADE's user guide [START_REF] Hascoët | Tapenade 2.1 user's guide[END_REF].

C.1.2 Automatic Differentiation example using TAPENADE

Many examples simple examples can be found in the litterature, in this section we illustrate the construction of the tangent and adjoint of a homemade short code source written in Fortran90. In this example, we consider a subroutine foo2, defined by 2, nested in a main subroutine foo1, defined by 6, and aim at differentiating the latter with respect to the active variables m1 ∈ R n and v4 ∈ R. where, These source codes are automatically transformed by TAPENADE and are verified by the validations tests § 2.4.7. The tangent mode results to 11 and 4, where the variables ending by d denotes the differential tangent variable which are always initialized at the begining of the code for security :

We recall that the reverse mode use the checkpointing SA strategy. As illustrated previously, this procedures consists in a forward computation of the source code and a backward computation to deduce the differential adjoints. The checkpointing is performed by repeatedly calling PUSH procedures which push the value at a given stage of the source algorithm during the forward procedure and POP procedures which reads the first value of the pile. The result of the reverse mode is given by 26 and 6. (v1, v1b, v2, v2b, v3, v3b)

m1b(i+1) = m1b(i+1) -SIN(m1(i+1))*v2b m1b(i) = m1b(i) + 2*m1(i)*v1b end for v4b = 0.0 Annexe D

Discrete adjoint of auxiliary procedures used in SWE

We describe in the tangent and adjoint procedures corresponding to each auxiliary routines used within the shallow water code SWE. We recall in algorthm 36 the structure of the original SWE code.

Algorithm 36 SWE(h 0 ,u 0 ,v 0 ) Initialization: n = 0, h 0 , u 0 , v 0 while n < N do call u2q(h n ,u n ,v n ,q n x ,q n y ) call RK3(h n ,q n

x ,q n y ,h n+1 ,q n+1

x ,q n+1 y ) call q2u(h n+1 ,u n+1 ,v n+1 ,q n+1 x ,q n+1 y ) save (h n+1 ,u n+1 ,v n+1 ) n = n + 1 end while return h N , u N , v N We will provide the numerical code or the equivalent matrix formulation for each routine laid out in this chapter. In terms of the notations, we inherit from the formalism imposed by TAPENADE such that :

-the tangent procedures are ended by the subscript _d, the differential tangent variable associated to a given variable x are denoted ẋ, -the adjoint procedures are ended by the subscript _b, the differential adjoint variable associated to a given variable x are denoted x.

D.1 Adjoint of the boundary conditions procedure

(i, j) = h(i, j) + hE (i, j) ū(i, j) = ū(i, j) + ūE (i, j) v(i, j) = v(i, j) + vE (i,
(i, j) = h(i, j) + hW (i + 1, j) ū(i, j) = ū(i, j) + ūW (i + 1, j) v(i, j) = v(i, j) + vW (i + 1,
h(1, j) = h(1, j) + hW (1, j) ū(1, j) = ū(1, j) -ūW (1, j) v(1, j) = v(1, j) + vW (1, j) end for h(1, ny) = h(1, ny) + hW (1, ny + 1) ū(1, ny) = ū(1, ny) -ūW (1, ny + 1) v(1, ny) = v(1, ny) -vW (1, ny + 1) return h, ū, v

D.1.3 Contribution of the adjoint variable in the north direction

Algorithm 39 BC_north_b( hN , ūN , vN )

Define value within the domain for i=1,nx do for j=1,ny do h(i, j) = h(i, j) + hN (i, j) ū(i, j) = ū(i, j) + ūN (i, j) v(i, j) = v(i, j) + vN (i, j) end for end for Sum values from the ghost cells into the boundaries of the domain for i=1,nx do h(i, ny) = h(i, ny) + hN (i, ny + 1) ū(i, ny) = u(i, ny) + ūN (i, ny + 1) v(i, ny) = v(i, ny) -vN (i, ny + 1) end for for j=1,ny do h(nx, j) = h(nx, j) + hN (nx + 1, j) ū(nx, j) = ū(nx, j) -ūN (nx + 1, j) v(nx, j) = v(nx, j) + vN (nx + 1, j) end for h(nx, ny) = h(nx, ny) + hN (nx + 1, ny + 1) ū(nx, ny) = ū(nx, ny) -ūN (nx + 1, ny + 1) v(nx, ny) = v(nx, ny) -vN (nx + 1, ny + 1) return h, ū, v

D.1.4 Contribution of the adjoint variable in the south direction

Algorithm 40 BC_south_b( hS , ūS , vS )

Define value within the domain for i=1,nx do for j=1,ny do h(i, j) = h(i, j) + hS (i, j + 1) ū(i, j) = ū(i, j) + ūS (i, j + 1) v(i, j) = v(i, j) + vS (i, j + 1) end for end for Sum values from the ghost cells into the boundaries of the domain for i=1,nx do h(i, 1) = h(i, 1) + hS (i, 1) ū(i, 1) = u(i, 1) + ūS (i, 1) v(i, 1) = v(i, 1) -vS (i, 1) end for for j=1,ny do h(nx, j) = h(nx, j) + hS (nx + 1, j + 1) ū(nx, j) = ū(nx, j) -ūS (nx + 1, j + 1) v(nx, j) = v(nx, j) + vS (nx + 1, j + 1) end for h(nx, 1) = h(nx, 1) + hS (nx + 1, 1) ū(nx, 1) = ū(nx, 1) -ūS (nx + 1, 1) v(nx, 1) = v(nx, 1) -vS (nx + 1, 1) return h, ū, v

D.1.5 Discrete adjoint of the BC routine

Given the adjoint procedures presented in this section, the discrete adjoint ( h, ū, v) of the boundary conditions are given by algorithm 41.

Algorithm 41 BC_b( h, ū, v, hW , ūW , vW , hE , ūE , vE , hN , ūN , vN , hS , ūS , vS ) The Roe routine which solves the 1D Riemann problem from the left and right state x L = (h L , u L , v L ) and x R = (h R , u R , v R ) is described by algorithm 42. We describe in this section the adjoint procedures of the auxiliary routines Priestley, eigenvalues, eigenvectors, fluxH and fluxV called from the Roe routine.

Algorithm 42 Roe(h L , u L , v L , h R , u R , v R , f lag) ! Preprocessing w L = Priestley (h L , u L , v L ) w R = Priestley (h R , u R , v R ) ∆w = w L -w R w = 1
2 (w L + w R ) ! Characteristics of the Roe matrix Ã(X R , X L ) call eigenvalues( w,e 1 , e 2 , e 3 ) call eigenvectors( w,r 1 , r 2 , r 3 ) call coefficients( w,∆w,α 1 , α 2 , α 3 ) ! Compute the left and right fluxes if f lag then

F L = fluxH (h L , u L , v L ) F R = fluxH (h R , u R , v R ) else F L = fluxV (h L , u L , v L ) F R = fluxV (h R , u R , v R ) end if return 1 2 (F L + F R ) -1 2 3 k=1 α k |e k |r k D.2.

Adjoint of the Preistley algorithm

In order to implement the Roe solver to the 2D SWE model, Priestley (1987) introduced parameter vectors which are determined by the left and right state vectors x L and x R (see § 3.3.2). Given an input x = (h, u, v), we define the Preistley routine by the following function

Preistley(x) = ( √ h, u √ h, v √ h)
The tangent procedure with respect to a linearization state x = (h, u, v) can be written into the following matrix form, Similarly to the eigenvalues, the eigenvectors r 1 , r 2 , r 3 of the linearized Jacobian Ã(x R , x L ) are computed from the input variable w := 1 2 (w L + w R ) . We recall below the expression of each eigenvectors,

r 1 =   w 1 w 2 -w 1 w 2 1 d w 3 -w 1 w 2 1 (1 -d)   , r 2 =   0 0 w 1   , r 3 =   w 1 w 2 + w 1 w 2 1 d w 3 + w 1 w 2 1 (1 -d)   ,
where d = 1 for the horizontal fluxes and d = 0 for the vertical fluxes. We provide the tangent and adjoint procedures with respect to a given linearization state w := (w 1 , w 2 , w 3 ).

The tangent and adjoint procedures corresponding to the calculation of the eigenvector r 2 are elementary, therefore we only provide the formulation for the eigenvectors r 1 and r 3 which are summarized by :

r 1,3 =   w 1 w 2 ± w 1 w 2 1 d w 3 ± w 1 w 2 1 (1 -d)   .
The tangent procedure can be written into the following matrix form

ẇ =   1 0 0 ± w 2 1 d 1 0 ± w 2 1 (1 -d) 0 1   ṙ1,

D.3 Adjoint for the conversion of non-conservative variables to conservative variables

The SWE routine requires two auxiliary routines u2q and q2u which transforms the non-conservative x = (h, u, v) T into the conservative variable state x = (h, h.u, h.v) T and vice versa. These routines are defined by elementary operations, therefore their corresponding tangent and adjoint routines are straightforwardly obtained with TA-PENADE.

In the one hand, the u2q routine is defined by the elementary operations q x = hu, (D.4) q y = hv, (D.5) its corresponding tangent procedure u2q_d is given by qx = h u + ḣu, (D.6) qy = h v + ḣv, (D.7)

and its corresponding adjoint procedure u2q_b is deduced by h = h + uq x + v qy (D.8) ū = hū (D.9) v = hv.

(D.10)

In the one hand, the q2u routine is defined by the elementary operations u = q x /h, (D.11) v = q y /h, (D.12) its corresponding tangent procedure q2u_d is given by Annexe E

u = h qx -

Discrete adjoint of the compact schemes

The first and second order derivatives are computed for each direction, after the swap to the adequate pencil. For the sake of simplicity, we consider within this section that swapped to the X-pencil configuration and we aim at computing the first and second order derivatives of a given vector f := (f i ) nx . We recall that the 6th order numerical scheme used to approximate the first order derivative f (x i ) at x i reads

αf i-1 + f i + αf i+1 = a f i+1 -f i -1 2∆x + b f i+2 -f i-2 4∆x .
and the approximation of the second order derivative f (x i ) at x i reads

αf i-1 + f i + αf i+1 = a f i+1 -2f -i + f i-1 2∆x 2 + b f i+2 -2f i + f i-2 4∆x 2 + c f i+3 -2f i + f i-3 9∆x 2 ,
The main asset of this numerical scheme is the fact that it's a linear algebraic combination, therefore the adjont procedure can be easiely deduced. It is possible to write the compact schemes in the following matrix forms :

A x f = 1 ∆x B x f, A x f = 1 ∆x B x f,
where the matrices A x , B x , A x and B x of size n x × n x depends on the selected boundary conditions. Hence, the adjoint procedure corresponding to the compact schemes are simply deduced by :

B T x f = 1 ∆x A T x f , (B x ) T f = 1 ∆x (A x ) T f .
Periodic condition conserves the order of the numerical scheme and is easy to implement. The values of f, f , f at the boundaries are simply deduced by :

f 0 = f nx , f -1 = f nx-1 , f 0 = f nx , f 0 = f nx .
Therefore the compact scheme can be formulated as : Free slip condition has the same assets than the periodic condition and corresponds whether to the symmetric condition

A x =             1 α α α 1 α . . . . . . . . . . . . α 1 α α α 1             and B x = 1 ∆x            
f 0 = f 2 , f -1 = f 3 , f 0 = f 2 , f 0 = f 2 ,
This boundary condition can be considered as a symetric or antisymetric condition according to the velocity component considered. If f is even, the matrix A x and B x reads : or the asymmetric condition

A x =             1
f 0 = -f 2 , f -1 = -f 3 , f 0 = f 2 , f 0 = -f 2 .
If f is odd, the matrix A x and B x reads : Open boundary condition doesn't make any assumptions on the flow outside the computational domain. Single sided formulations are used for the approximation of first derivatives and second derivatives for these types of boundaries using relations of the form :

A x =            
f 1 + 2f 2 = 1 2∆x
(-5f 1 + 4f 2 + f 3 ),

f 1 + 11f 2 = 1 2∆x 2 (13f 1 -27f 2 + 15f 3 -f 4 ),
that are third order accurate [START_REF] S K Lele | Compact finite difference schemes with spectral-like resolution[END_REF]. At the adjacent nodes, because a three point formulation must be used Padé schemes are employed with forth-order accurate schemes :

1 4 f 1 + f 2 + 1 4 f 3 = 3 2 f 3 -f 1 2∆x , 1 10 f 1 + f 2 + 1 10 f 3 = 6 5 f 3 -2f 2 + f 1 ∆x . A x =             1

PIV system

In the double frame configuration, we can change each system's acquisition time which corresponds to the period between two snapshots which are used to determine the velocity of each particules. Then we set a delay between two acquisition in order to perform another acquisition with system 2.In figure G.1, the acquisition of the system 1 and system 2 is 500µs and 400µs, respectively. We set a cooldown of 2000µs between each acqusition. The period between the begin of the system 1 and system 2 velocity acquisition is 1000µs. Thus, in this configuration we can set the observation time step ∆t obs = 1000µs. In the single frame configuration, each system's acquisition time is fixed, we perform a shit between each acquisition to avoid an overlap between each system's acquisition. In figure G.2, the acquisition of both is 2000µs. We set system 2's acquisition 1000µs after system 1's acquisition. Thus, in this configuration we can set the observation time step ∆t obs = 1000µs.

The stereo PIV systems were calibrated beforehand to ensure that the common area of both observed planes are equivalent. The correct spatio-temporal positioning of the common area was verified by comparing measurements of the velocity profile in the common area for each PIV system as illustrated in figure G.3. The comparison was performed by A. Guibert.

As the wind tunnel was initially designed for a mixing layer configuration, we

Résumé

L'étude des écoulements fondamentaux turbulents est au coeur de nombreux problèmes académiques et industriels. Une des caractéristiques principales des écoulements turbulents est la sensibilité par rapport aux structures turbulentes aux petites échelles. En effet, ces dernières sont capables de modifier la trajectoire de l'écoulement de manière conséquente, ainsi, il est nécessaire d'avoir une bonne estimation de l'apport de ces structures afin de reconstruire un écoulement turbulent réaliste.

D'une part, les équations de Navier-Stokes permettent de décrire l'évolution de tels écoulements dans le temps et dans l'espace. La littérature propose une multitude de méthodes numériques permettant la résolution de telles équations. La simulation directe numérique (DNS) permet la reconstruction de toutes les échelles turbulentes mais nécessite un coût de calcul considérable même pour des écoulements faiblement turbulents. Des modèles de la turbulence (LES, RANS) permettent de simuler les structures à grandes échelles des écoulements fortement turbulents à un coût de calcul raisonnable tout en incluant la contribution des structures aux petites échelles. Malgré la précision des solveurs numériques actuels, les modèles physiques sont souvent sujets à des erreurs, typiquement dues aux approximations des conditions de bords du domaine de calcul, qui peuvent nous éloigner de la solution réelle. Cette approche consiste à minimiser une fonction coût qui traduit l'écart entre l'ébauche et les observations sous contrainte du modèle dynamique du système considéré. Le problème de minimisation est résolu par une méthode itérative quasi-Newtonienne qui nécessite le calcul du gradient de la fonction coût. Une manière élégante, précise et peu coûteuse de calculer ce gradient consiste à résoudre le modèle adjoint associé au modèle dynamique du système. La construction et l'implémentation du modèle adjoint représente un effort conséquent qui peut être facilité par l'utilisation d'outils de différentiation automatique. Ainsi, cette thèse s'articule autour de la résolution des problèmes de mécanique des fluides par une méthode variationelle, dénommée par la suite 4DVar.

Enfin, ces dernières années ont été marqué par l'apparition de méthodes hybrides qui combinent les deux approches précédentes afin de tirer parti des avantages et compenser les inconvénients de chaque approche. Bien que ces approches soient prometteuses, elles requièrent un effort de parallèlisation afin de concurrencer les approches classiques en termes de coût de calcul. Malheureusement, cette parallèlisation peut se révéler difficile à implémenter dans le cas de systèmes dynamiques complexes.

Dans un premier temps, nous avons mené une première application de la 4DVar à la reconstruction d'un écoulement simple à la surface libre d'un fluide. Cet écoulement est modélisé par les équations shallow-water (SW) qui déterminent l'évolution de la hauteur et du champ de vitesse d'un écoulement peu profond. Nous avons construit le processus adjoint correspondant au modèle SW à l'aide de l'outil de différentiation automatique TAPENADE. Le code adjoint obtenu a été sujet à un processus de validation systématique afin d'assurer la correspondance au code dynamique. Au cours de cette thèse, nous avons eu l'opportunité de comparer les performances d'une méthode variationelle classique (4DVar) avec plusieurs versions d'une méthode d'assimilation hybride (4DEnVar). Les méthodes hybrides ont l'avantage de ne pas avoir recours à la construction de l'adjoint, cependant elles ont nécessité un effort de parallèlisation afin de concurrencer la méthode 4DVar.

Une première reconstruction a été effectuée à partir d'observations synthétiques complètes, puis partielles du système. Dans ce cas de figure, nous comparons les résultats des différentes techniques d'assimilations par rapport à un état de référence donné. L'ébauche a été construite en modifiant l'inclinaison initiale de la surface libre étudiée. Les observations synthétiques ont été générées en ajoutant un bruit purement synthétique à l'écoulement de référence. Dans cette application synthétique, les méthodes hybrides ont apportées une légère amélioration par rapport à la méthode 4DVar à un coût de temps de calcul réduit grâce à la parallélisation effectuée.

Par la suite, nous avons mené la reconstruction à partir d'observations de la hauteur de la surface libre d'une eau colorée contenue dans un bac rectangulaire. Ces observations ont été prélevées par un capteur de profondeur (capteur Kinect) et sont caractérisées par un bruit lié aux erreurs de mesures du capteur et à l'absence de données dans des zones situées à la paroi et localement à l'intérieur du domaine. Ces zones d'absence de données sont dues à la réflexion à la surface libre de l'eau colorée observée. Nous avons mis en oeuvre un algorithme de prétraitement des observations et ajusté les matrices de covariance aux erreurs d'observations en conséquence. En termes d'ébauche, la hauteur initiale a été interpolée spatialement à partir de l'observation initiale, et le champ de vitesse initiale a été fixé arbitrairement. Etant donné que nous ne possédons pas d'écoulement de référence, les comparaisons ont été effectuées uniquement par rapport aux observations de la hauteur. Dans ce cas de figure, la reconstruction par la 4DVar présente des discontinuités entre les zones dans lesquelles les observations sont manquantes alors que la reconstruction par les méthodes hybrides conduit à une solution plus lissée dans ces régions.

Dans un deuxième temps, nous avons mené la reconstruction volumique de l'aval d'un sillage de cylindre à Reynolds 300. L'écoulement turbulent a été simulé par un code DNS Incompact3D. Ce dernier combine des schémas numériques compacts et spectrales afin d'assurer une précision suffisamment élevée pour simuler toutes les échelles de la turbulence. La parallèlisation est implémentée selon une décomposition de domaine afin de réduire le temps de calcul de la DNS. Ces éléments ont été pris en compte dans la construction du code adjoint associé qui s'est fait à l'aide de l'outil de différentiation automatique TAPENADE. Une procédure de validation rigoureuse a été systématiquement effectuée afin d'assurer que le code adjoint conserve la précision du code dynamique originel.

Une première reconstruction a été effectuée à partir d'observations volumiques synthétiques. Afin de réduire le coût de calcul nécessaire à l'application de la 4DVar, nous avons mené une étude préliminaire sur l'impact des erreurs engendrées par le modèle dynamique dans un domaine de calcul réduit. Les erreurs liées aux conditions de bords ont été prises en compte dans la reconstruction en ajustant les matrices de covariances aux erreurs d'observations. La condition initiale de l'ébauche a été ensuite générée en ajoutant un bruit synthétique et un biais important à la condition initiale de référence. Ce choix a pour objectif d'évaluer la capacité de la 4DVar à corriger la condition initiale à partir de quelques observations sur une fenêtre d'assimilation courte. Les résultats ont montré que l'assimilation a permis de corriger le biais imposé à l'ébauche, cependant le bruit synthétique initial a été peu corrigé. De plus, l'ajustement des matrices de covariances ont apporté une contribution qui peut améliorer la condition initiale.

Par la suite, nous avons effectué la reconstruction volumique à partir d'observations de plans orthogonaux par la stéréo-PIV. Une première reconstruction a été menée à partir d'observations synthétiques qui sont déduites par l'ajout d'un bruit synthétique à l'écoulement de référence. Il a été nécessaire de construire un algorithme de prétraitement qui génère des observations volumiques à partir d'observations plans. Les matrices de covariances des erreurs d'observations ont été construites en conséquence. Afin de reproduire les conditions similaires au cas expérimental, l'ébauche a été générée à partir d'un écoulement de vitesse nul en ajoutant une condition d'entrée bruitée. Au bout d'un certain nombre de pas de temps, la comparaison entre la référence et l'ébauche ont montré des similarités très proches qui ne justifiaient pas l'utilisation de la 4DVar. En revanche, dans le cas expérimental, l'ébauche s'est avérée beaucoup plus lissée par rapport aux observations, ce qui motive l'application de la 4DVar. De la même manière que dans la première application au problème de SW, les résultats de la 4DVar ont été comparés aux plans observés puisque nous ne possédons pas d'écoulement de référence. La condition initiale d'analyse obtenue démontre que la 4DVar est capable de modifier les structures à grandes échelles de l'ébauche initiale afin de correspondre aux observations réelles.

Avec le développement récent et prometteur des techniques de mesures volumiques, les résultats de cette thèse ouvrent la voie à une application des méthodes d'assimilations de données à des problèmes en trois dimensions. L'étude des erreurs engendrées par le modèle ont également permis d'améliorer la stratégie de la 4DVar, notamment en termes d'ajustement des matrices de covariance des erreurs de modèle et d'observations.

4. 3

 3 Validation on a synthetic case with 3D observations . . . . . . . . . . 4.3.1 Synthetic configuration . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Initialization from a pre-existing turbulent flow . . . . . . . . 4.3.3 Synthetic observations . . . . . . . . . . . . . . . . . . . . . . 4.3.4 Background choice . . . . . . . . . . . . . . . . . . . . . . . . 4.3.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.6 Conclusion and discussion . . . . . . . . . . . . . . . . . . . . 4.4 Towards the reconstruction of a full 3D velocity field with orthogonalplane stereo PIV observations . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Experimental context . . . . . . . . . . . . . . . . . . . . . . . 4.4.2 Synthetic orthogonal plane stereo PIV observations . . . . . . 4.4.3 Real orthogonal-plane stereo PIV observations . . . . . . . . . 4.4.4 Conclusion and discussion . . . . . . . . . . . . . . . . . . . . Differentiation principles . . . . . . . . . . . . . . . . . . . C.1.1 Automatic Differentiation implementation . . . . . . . . . . . C.1.2 Automatic Differentiation example using TAPENADE . . . . D Discrete adjoint of auxiliary procedures used in SWE D.1 Adjoint of the boundary conditions procedure . . . . . . . . . . . . . D.1.1 Contribution of the adjoint variable in the east direction . . . D.1.2 Contribution of the adjoint variable in the west direction . . . D.1.3 Contribution of the adjoint variable in the north direction . . D.1.4 Contribution of the adjoint variable in the south direction . . D.1.5 Discrete adjoint of the BC routine . . . . . . . . . . . . . . . . D.2 Adjoint of the routines involved in the Roe solver . . . . . . . . . . . D.2.1 Adjoint of the Preistley algorithm . . . . . . . . . . . . . . . D.2.2 Adjoint of the eigenvalues of Ã(x R , x L ) . . . . . . . . . . . . . D.2.3 Adjoint of the eigenvectors of Ã(x R , x L ) . . . . . . . . . . . . D.2.4 Adjoint of the nonlinear flux . . . . . . . . . . . . . . . . . . . D.3 Adjoint for the conversion of non-conservative variables to conservative variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E Discrete adjoint of the compact schemes F Discrete adjoint of FFT_postproc G PIV system List of figures 1.1 Kolmmogorov energy spectrum representing all the scales of a turbulent flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Ludwig Prandtl in 1904 showing a test bench, known as the Prandtl channel, to visualize the flow evolution. . . . . . . . . . . . . . . . . . 1.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.2 1D representation of the piecewise partition of the solution. The finite volume solution has a constant value within each cell C i = [x i-1/2 ; x i+1/2 ]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.3 Reflecting boundary conditions : we set the value of the flow at ghost cells (in gray) with respect to the values of the height and velocity fields at the boundary cells. . . . . . . . . . . . . . . . . . . . . . . . 67 3.4 Left side : BC creates ghost cells from the input state x (black box) to x E , x W , x N and x S (red box) which are used to compute the horizontal flux F E , F W , G N and G S and G i+1/2,j respectively. Right side : BC_b sums the values at the ghost cells with the values at the boundaries of each adjoint variable state x (black box) for each adjoint xE , xW , xN and xS (green box). . . . . . . . . . . . . . . . . . . . . . 77 3.5 Representation of the initial configuration of the reference (on the left) and the background (on the right). The reference is obtained by a smooth 21% slope along the x-axis and a 10% slope along the y-axis. The background height is set as a smooth slope at 20% slope along the x-axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.6 Temporal setup for the data assimilation techniques : we perform the assimilation over a period of t f •U/L x . The observations are injected each ∆t obs = 50∆t•U/L x . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.7 RMSE comparison between an incremental 4DVar and 4DEnVar assimilation approaches with a a partial height observation through a noisy free surface height. . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.8 RMSE comparison between an incremental 4DVar and 4DEnVar assimilation approaches with a a partial velocity observation through a noisy free surface height. . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.9 Experimental setup with the Kinect depth sensor . . . . . . . . . . . 88 3.10 Observed free surface height retrieved from the Kinect depth map (a), interpolated free surface height used as the background (b), background velocity set accordingly to the height observation (c) . . . . . 89 3.11 Temporal setup for the sliding data assimilation techniques : we perform the assimilation over a period of t f •U/L x . The observations are injected each ∆t obs = 30∆t•U/L x . We then slide the assimilation window to the next observation and perform once again the data assimilation. This procedure is repeated until we reach the ninth observation. 90 3.12 Comparison of the evolution of the mean surface height of the wave crest region.a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 3.13 Height field comparison at t•U/L x ≈ 0.0652 (in the left) and at t•U/L x ≈ 0.5859 (in the right). From top to bottom : observation, background, 4DVar, En4DVar proposed by Liu et al. (2008), the 4DEnVar-OL-LC and the 4DEnVar-OL-LE. . . . . . . . . . . . . . . . 92 4.1 2D domain decomposition example of a 3D computation domain using a 4 × 4 processor grid : (a) X-pencil, (b) Y-pencil and (c) Z-pencil. . . 96 4.2 Representation of the velocity meshgrid (black points) and the pressure meshgrid (white points). The pressure meshgrid is obtained by shifting by a half-mesh, ∆x/2 and ∆y/2, along the spatial direction x and y respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.3 General structure of the DNS code Incompact3d. The pencil configuration swaps carried out in each subroutine are specified in the right column. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.4 Discrete adjoints of global MPI communicators . . . . . . . . . . . . . 102 4.5 Snapshot of the stream-wise velocity field u/U in the reference domain Ω ref . The turbulent structures are obtained from an initial null velocity field and the input of an inflow velocity over T ref = 34000∆t . The assimilation is run within the sub-domain Ω (in blue) and we perform statistics of the results in the the sub-domain Ω obs (in green). 113 4.6 Temporal setup for the reconstruction of a complete three dimensional flow with the 4DVar technique. In a first step we build the reference flow of this study by running a DNS over T ref = 34000∆t = 170 within the reference domain Ω ref . We then extract the three components of the velocity at t 1 = 170 and add a 15% noise and a sin(y 3 ) profile. The latter is used to initialize the 4DVar code performed over the assimilation window T a = 50∆t = 0.25, by injecting three dimensional synthetic observations every ∆t obs = 10∆t . . . . . . . . . . . . . . . . 114 4.7 Error between the reference and the simulated stream-wise velocity field at T ref + 50∆t = 170.25 and T ref + 500∆t = 172.5 at the streamwise plane z = 3D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 4.8 Representation of the inverted observation covariance matrix R -1 at the stream-wise plane z = 3D. We reduce gradually the confidence within the boundary zones by the application of several 1D Gaussian profiles described by (4.23), (4.24) and (4.25). . . . . . . . . . . . . . 117 4.9 Assimilation with 3D synthetic observations : comparison between the stream-wise velocity field of the reference flow (left) with the background flow (right) obtained by adding a 15% white noise and a sin(y 3 ) profile bias along the y-axis at t 1 = 170. . . . . . . . . . . . . 117 4.10 Assimilation with 3D synthetic observations : comparison of the evolution of the cost function and cost function gradient for the case a (on the left) and the case b (on the right) with respect to the inner loop optimization iterations. . . . . . . . . . . . . . . . . . . . . . . . 118

  123 4.15 General view of the mixing layer wind tunnel . . . . . . . . . . . . . . 124 4.16 Representation of the 2D3C PIV systems : SPIVS1 (in red) measures the three components velocity at the plane x = 19D and PIVS2 (in green) measures the three components velocity at the plane z = 0D. . 125 4.17 Single frame configuration : each stereo PIV system requires a frequency dt to acquire two consecutive frames. In the single frame configuration, dt = 1000µs. . . . . . . . . . . . . . . . . . . . . . . . . 125 4.18 Assimilation with synthetic orthogonal-plane observations : temporal setup for the reconstruction of a complete three dimensional flow with the 4DVar technique. In a first step we build the reference flow of this study by running a DNS over T ref = 30000∆t within the reference domain Ω ref . In a second step, we generate an initial background from a null velocity by injecting the three components of the velocity field every ∆t PIV = 20∆t at the inlet x = 0D within the sub-domain Ω. After a T b = 5000∆t time period, the initial background is used to initialize the 4DVar code performed over the assimilation window T a = 60∆t = 0.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 4.19 Assimilation with synthetic orthogonal-plane observations : comparison of the initial reference and synthetic observation stream-wise velocity u/U at the observed planes x = 0D (top) and z = 3D (bottom) at t 2 = 175. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 4.20 Assimilation with synthetic orthogonal-plane observations : representation of the inverted observation covariance matrix R -1 in the plane y = 4D, z = 3D and x = 0D . . . . . . . . . . . . . . . . . . . . . . . 129 4.21 Assimilation with synthetic orthogonal-plane observations : comparison between the streamwise velocity of the reference flow (left) with the background flow (right) after T b = 5000∆t . The background is obtained by injecting the temporal integration of the exact inflow at each ∆t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Assimilation with synthetic orthogonal-plane observations : comparison between the strea-mwise velocity of the reference flow (left) with the background flow (right) after T b = 5000∆t . The background is obtained by injecting the temporal integration of the a noisy inflow at each ∆t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.23 Assimilation with synthetic orthogonal-plane observations : evolution of the cost function (left) and cost function gradient (right) with respect to the inner loop iterations for two outer loop iterations. . . . . 4.24 Assimilation with synthetic orthogonal-plane observations : comparison of the norm of each component of the velocity field from t 2 = 175 to t 3 = 175.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.25 Assimilation with real orthogonal-plane observations : temporal setup for the reconstruction of a complete three dimensional flow with the 4DVar technique. In a first step we build the reference flow of this study by running a DNS over T ref = 30000∆t = 150 within the reference domain Ω ref . In a second step, we generate an initial background from a null velocity and by injecting alternatively the three components of the velocity field every ∆t obs = 10∆t at the plane x = 0D and z = 3D. After a T b = 5000∆t time period, initial background is used to initialize the 4DVar code performed over the assimilation window T a = 50∆t = 0.25. . . . . . . . . . . . . . . . . . . . . . . . 4.26 Assimilation with real orthogonal-plane observations : comparison of the mean and the standard deviation of each component of the velocity field in time at a punctual position on the inlet plane at y = 4D . 4.27 Assimilation with real orthogonal-plane observations : representation of the observation covariance matrix R -1 in the plane y = 4D, z = 3D and x = 0D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.28 Assimilation with real orthogonal-plane observations : comparison of the initial observation (left) and background (right) stream-wise velocity u/U at the observed stream-wise plane z = 3D. The background is obtained after running the DNS over a 4000∆t period from a null velocity field and the input of the temporal interpolation of the observation inlet at each ∆t . . . . . . . . . . . . . . . . . . . . . . . . . 4.29 Assimilation with real orthogonal-plane observations : cost function evolution with respect to the inner loop iterations. . . . . . . . . . . . 4.30 Assimilation with real orthogonal-plane observations : comparison of of the norm of each component of the velocity field from t 1 = 20 to t 2 = 20.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.31 Assimilation with real orthogonal-plane observations : snapshots of the velocity field for each component u, v and w at the beginning of the assimilation window at t 1 = 20 in the stream-wise plane z = 3D. Each row corresponds to the observation, the background and the analysis, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.32 Assimilation with real orthogonal-plane observations : snapshots of the velocity field for each component u, v and w at the beginning of the assimilation window at t 1 = 20 in the stream-wise plane z = 3D. Each row corresponds to the observation, the background and the analysis, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.33 Assimilation with real orthogonal-plane observations : snapshots of the velocity field for each component u, v and w at the beginning of the assimilation window t 1 = 20 in the stream-wise plane y = 4D. Each row corresponds to the observation, the background and the analysis, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . C.1 SA strategy (top) stores all the intermediate values (black points) at the end of each instructions {I k ; k = 1, • • • , p -1}. In the reverse mode, each values are invoked and used. Checkpointing SA strategy (bottom) consists in applying the SA strategy on small blocs of the program. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.2 RA strategy (top) performs all the intermediate instructions from the initial value x (black points) to obtain the required value. Checkpointing RA strategy (bottom) recomputes the intermediate values from each snapshot (black points) . . . . . . . . . . . . . . . . . . . . . . . G.1 Double frame configuration : the acquisition of the system 1 and system 2 is dt1 and dt2 respectively. . . . . . . . . . . . . . . . . . . . . G.2 Single frame configuration . . . . . . . . . . . . . . . . . . . . . . . . G.3 Comparison of the velocity profile u at the inlet . . . . . . . . . . . . G.4 Dual plane position in the cylinder wake configuration at Reynolds List of tables 3.1 Evaluation of 1 -R 1 results for different values of the differentiation step α over 10∆t and 200∆t . These results are calculated in MAT-LAB's default accuracy. . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 4DVar validation result for different values of the differentiation step α over an assimilation window T = 200∆t . . . . . . . . . . . . . . .

4. 1

 1 List of the clusters used in this thesis . . . . . . . . . . . . . . . . . . 4.2 Tangent validation results for different values of the differentiation step α over 10∆t and 100∆t . . . . . . . . . . . . . . . . . . . . . . . 4.3 4DVar validation result for different values of the differentiation step α in FORTRAN single and double precision. . . . . . . . . . . . . . .

Figure 1

 1 Figure 1.1 -Kolmmogorov energy spectrum representing all the scales of a turbulent flow

Figure 1

 1 Figure 1.2 -Ludwig Prandtl in 1904 showing a test bench, known as the Prandtl channel, to visualize the flow evolution.

Figure 2

 2 Figure 2.1 -Sequential data assimilation representation : we seek for the analysis state x a k+1 that corrects the forecast state x f k+1 generated by the dynamics model given the observation Y k+1 at the instant t k+1 . The dynamics model start from the analysis state x a k of a previous assimilation process.

Figure 2

 2 Figure 2.2 -Ensemble Kalman Filter representation : we consider a set, called ensemble, of possible state samples, each referred to as an ensemble member, noted {x

Figure 2

 2 Figure 2.3 -Standard 4DVar assimilation representation : we aim at determining the optimal assimilation trajectory (red curve) from the background trajectory (blue curve) and a set of observations (green). We minimize the distance between the background trajectory and the observations by a gradient-descent algorithm. The cost function gradient is deduced by performing an adjoint procedure backwards in time (orange curve) which takes into account the observations correction.

  x

Figure 2

 2 Figure 2.4 -Incremental variational assimilation representation : we seek for the optimal increment δx (1) 0 of the linearized convex cost function around the initial state x (1) 0 = x b 0 (blue curve) and unpdate to x

repeat

  Compute δx t with the tangent linear dynamics (2.52) and (2.53) Compute the cost function (2.62) Compute the adjoint variable λ(t 0 ) with the backward integration of relation (2.47)

Figure 2

 2 Figure 2.5 -Construction of the differentiation code df(x) associated to the analytic function y = f (x) by the discretize-then-differentiate (solid line) and the differentiate-then-discretize (dashed line) strategies.

  . De Pando et al. (2012) evaluated the direct and adjoint linearized dynamics from compressible flow solvers. In another context, Ponçot et al. (2013) constructed the discrete adjoint to implement a variational technique for xenon dynamical forecasts in neutronic. Gronskis et al. (2013) employed for the first time the discrete adjoint in the implementation of a variational technique which combines the 2D version of Incompact3d and PIV measurements. This strategy has also the great advantage to guarantee that the constructed adjoint is the adjoint of the discrete code.

  3)The equations (3.1), (3.2) and (3.3) can be written in a conservative form by considering the conserved variables as x =

Figure 3

 3 Figure 3.2 -1D representation of the piecewise partition of the solution. The finite volume solution has a constant value within each cell C i = [x i-1/2 ; x i+1/2 ].

Figure 3

 3 Figure 3.3 -Reflecting boundary conditions : we set the value of the flow at ghost cells (in gray) with respect to the values of the height and velocity fields at the boundary cells.

Algorithm 9 Derive

 9 Ensemble-based variational data assimilation algorithm : Localize covariance approach Set an initial state x(t 0 ) = x b (t 0 ) and the ensemble x j,b 0 as an arbitrary choice (for the 1st cycle) or as the forecast state an the ensemble forecast derived from the previous assimilation cycle respectively loop Define the matrixX(t) = [x(t 0 ), x 1,b 0 , ..., x N,b 0 ]concatenating the initial ensemble and compute in parallel X(t) with the forward integration of the nonlinear dynamics (3.14) Generate ensemble observations (3.29) and innovations (3.30) Derive the background perturbation matrix localization technique from (3.39) Initialize the increment matrix δX 0 = [δx(t 0 ), δx 1,b 0 , ..., δx N,b 0 ] Do an inverse control variable transformation δR 0 = (P b ) -1 δX 0 if necessary where (P b ) -1 is calculated by SVD Optimize in parallel δR 0 in the inner loop, the cost function and the gradient are calculated based on the modified versions of (3.23) Update control initial space R 0 and calculate δ 0 X 0 by transforming δR 0 to the state space with (3.22) Update the initial ensemble (3.33) and the ensemble perturbation matrix (3.34) end loop Evolve the analysis state x a (t 0 ) to the beginning of the next cycle through the nonlinear dynamics (3.14). The forecast state and the forecast enesemble are used to initialize the next assimilation cycle. Algorithm 10 Ensemble-based variational data assimilation algorithm : Localize ensemble approach Set an initial state x(t 0 ) = x b (t 0 ) and the ensemble x b 0 as an arbitrary choice (for the 1st cycle) or as the forecast state an the ensemble forecast derived from the previous assimilation cycle respectively loop Define the matrix X(t) = [x(t 0 ), x 1,b 0 , ..., x N,b 0 ] concatenating the initial ensemble and compute in parallel X(t) with the forward integration of the nonlinear dynamics (3.14) Define local space and the transformation from state space to local space (3.40) Parallelizing minimization computation at each grid point (p, q) Initialize the increment vector δx 0,l = M l,s δx 0 Apply an inverse control variable transformation in the local space

Figure 3

 3 Figure 3.4 -Left side : BC creates ghost cells from the input state x (black box) to x E , x W , x N and x S (red box) which are used to compute the horizontal flux F E , F W , G N and G S and G i+1/2,j respectively. Right side : BC_b sums the values at the ghost cells with the values at the boundaries of each adjoint variable state x (black box) for each adjoint xE , xW , xN and xS (green box).

  4DVar validation result for different values of the differentiation step α over an assimilation window T = 200∆t .

Figure 3

 3 Figure 3.6 -Temporal setup for the data assimilation techniques : we perform the assimilation over a period of t f •U/L x . The observations are injected each ∆t obs = 50∆t•U/L x .

Figure 3

 3 Figure3.8 -RMSE comparison between an incremental 4DVar and 4DEnVar assimilation approaches with a a partial velocity observation through a noisy free surface height.

Figure 3

 3 Figure 3.9 -Experimental setup with the Kinect depth sensor

Figure 3

 3 Figure 3.10 -Observed free surface height retrieved from the Kinect depth map (a), interpolated free surface height used as the background (b), background velocity set accordingly to the height observation (c)

Figure 3 .

 3 Figure 3.12 -Comparison of the evolution of the mean surface height of the wave crest region.a

  Figure 3.13 -Height field comparison at t•U/L x ≈ 0.0652 (in the left) and at t•U/L x ≈ 0.5859 (in the right). From top to bottom : observation, background, 4DVar, En4DVar proposed by Liu et al. (2008), the 4DEnVar-OL-LC and the 4DEnVar-OL-LE.

Figure 4

 4 Figure 4.1 -2D domain decomposition example of a 3D computation domain using a 4 × 4 processor grid : (a) X-pencil, (b) Y-pencil and (c) Z-pencil.

Figure 4

 4 Figure 4.2 -Representation of the velocity meshgrid (black points) and the pressure meshgrid (white points). The pressure meshgrid is obtained by shifting by a halfmesh, ∆x/2 and ∆y/2, along the spatial direction x and y respectively.

Figure 4

 4 Figure 4.3 -General structure of the DNS code Incompact3d. The pencil configuration swaps carried out in each subroutine are specified in the right column.

Figure 4 . 4 -

 44 Figure 4.4 -Discrete adjoints of global MPI communicators

Figure 4

 4 Figure 4.5 -Snapshot of the stream-wise velocity field u/U in the reference domain Ω ref . The turbulent structures are obtained from an initial null velocity field and the input of an inflow velocity over T ref = 34000∆t . The assimilation is run within the sub-domain Ω (in blue) and we perform statistics of the results in the the sub-domain Ω obs (in green).

Figure 4

 4 Figure 4.6 -Temporal setup for the reconstruction of a complete three dimensional flow with the 4DVar technique. In a first step we build the reference flow of this study by running a DNS over T ref = 34000∆t = 170 within the reference domain Ω ref .We then extract the three components of the velocity at t 1 = 170 and add a 15% noise and a sin(y 3 ) profile. The latter is used to initialize the 4DVar code performed over the assimilation window T a = 50∆t = 0.25, by injecting three dimensional synthetic observations every ∆t obs = 10∆t .

  Figure 4.7 -Error between the reference and the simulated stream-wise velocity field at T ref + 50∆t = 170.25 and T ref + 500∆t = 172.5 at the stream-wise plane z = 3D.

Figure 4

 4 Figure 4.8 -Representation of the inverted observation covariance matrix R -1 at the stream-wise plane z = 3D. We reduce gradually the confidence within the boundary zones by the application of several 1D Gaussian profiles described by (4.23), (4.24) and (4.25).

Figure 4

 4 Figure 4.10 -Assimilation with 3D synthetic observations : comparison of the evolution of the cost function and cost function gradient for the case a (on the left) and the case b (on the right) with respect to the inner loop optimization iterations.

Figure 4

 4 Figure4.11 -Assimilation with 3D synthetic observations : comparison of the RMSE evolution of the background (in blue), the observations (in green), the analysis in case a (in red) and the analysis in case b (in black) with respect to the reference state for each velocity components over the assimilation window of T a = 50∆t within the sub-domain Ω obs .

Figure 4

 4 Figure 4.12 -Assimilation with 3D synthetic observations : snapshots of the velocity field for each component u, v and w at the beginning of the assimilation window at t 1 = 170 in the stream-wise plane z = 3D. The row correspond to the reference, background, 4DVar-all and 4DVar-rampe flow, respectively.

Figure 4

 4 Figure 4.13 -Assimilation with 3D synthetic observations : snapshots of the velocity field for each component u, v and w at the end of the assimilation window at t 2 = 170.25 in the stream-wise plane z = 3D. The row correspond to the reference, background, 4DVar-all and 4DVar-rampe flow, respectively.

Figure 4 .

 4 Figure 4.14 -Assimilation with 3D synthetic observations : comparison of the energy spectrum evolution of the reference flow (in blue), the background flow (in green), the 4DVar-all flow (in red) and the 4DVar-rampe black (in black) at t 1 = 170 (on the left) and at t 2 = 170.25 (on the right) in a cross-section plane x = 4D. The comparison of the spectra of the error for the same data with respect to the reference are shown in inset.

Figure 4 .

 4 Figure 4.15 -General view of the mixing layer wind tunnel

Figure 4 .

 4 Figure 4.16 -Representation of the 2D3C PIV systems : SPIVS1 (in red) measures the three components velocity at the plane x = 19D and PIVS2 (in green) measures the three components velocity at the plane z = 0D.

Figure 4 .

 4 Figure 4.17 -Single frame configuration : each stereo PIV system requires a frequency dt to acquire two consecutive frames. In the single frame configuration, dt = 1000µs.

Figure 4 .

 4 Figure 4.18 -Assimilation with synthetic orthogonal-plane observations : temporal setup for the reconstruction of a complete three dimensional flow with the 4DVar technique. In a first step we build the reference flow of this study by running a DNS over T ref = 30000∆t within the reference domain Ω ref . In a second step, we generate an initial background from a null velocity by injecting the three components of the velocity field every ∆t PIV = 20∆t at the inlet x = 0D within the sub-domain Ω. After a T b = 5000∆t time period, the initial background is used to initialize the 4DVar code performed over the assimilation window T a = 60∆t = 0.3.

Figure 4 . 2 1 4 :

 424 19 shows two snapshots of the reference and the observed flows at each observed plane at t 2 .Algorithm 28 Synthetic observations pre-processing in the orthogonal-plane PIV observations configuration 1: Add a 10% white noise to each component of the observed velocity field at each occurence ∆t PIV and at each point of the grid 2: At each plane x = 0D and z = 3D perform a temporal interpolation at every ∆t obs 3: Set a weight for the inlet plane w 1 (x) = exp -x 2 2σ Set a weight for the stream-wise plane w 2 (z) = exp -(z-3D) Compute the value of each component of the velocity flow u(x, y, z) at each non-observed point (x, y, z) from the values of the observed velocities at the inlet u 1 (y, z) and at the stream-wise plane u 2 (x, y)

Figure 4

 4 Figure 4.19 -Assimilation with synthetic orthogonal-plane observations : comparison of the initial reference and synthetic observation stream-wise velocity u/U at the observed planes x = 0D (top) and z = 3D (bottom) at t 2 = 175.

Figure 4 .

 4 Figure 4.20 -Assimilation with synthetic orthogonal-plane observations : representation of the inverted observation covariance matrix R -1 in the plane y = 4D, z = 3D and x = 0D

Figure 4 .Figure 4

 44 Figure 4.22 -Assimilation with synthetic orthogonal-plane observations : comparison between the strea-mwise velocity of the reference flow (left) with the background flow (right) after T b = 5000∆t . The background is obtained by injecting the temporal integration of the a noisy inflow at each ∆t .

Figure 4 .

 4 Figure 4.24 -Assimilation with synthetic orthogonal-plane observations : comparison of the norm of each component of the velocity field from t 2 = 175 to t 3 = 175.3.

Figure 4 .

 4 Figure 4.25 -Assimilation with real orthogonal-plane observations : temporal setup for the reconstruction of a complete three dimensional flow with the 4DVar technique. In a first step we build the reference flow of this study by running a DNS over T ref = 30000∆t = 150 within the reference domain Ω ref . In a second step, we generate an initial background from a null velocity and by injecting alternatively the three components of the velocity field every ∆t obs = 10∆t at the plane x = 0D and z = 3D. After a T b = 5000∆t time period, initial background is used to initialize the 4DVar code performed over the assimilation window T a = 50∆t = 0.25.

Figure 4 .

 4 Figure 4.26 -Assimilation with real orthogonal-plane observations : comparison of the mean and the standard deviation of each component of the velocity field in time at a punctual position on the inlet plane at y = 4D

Figure 4

 4 Figure 4.27 -Assimilation with real orthogonal-plane observations : representation of the observation covariance matrix R -1 in the plane y = 4D, z = 3D and x = 0D.

Figure 4

 4 Figure 4.29 -Assimilation with real orthogonal-plane observations : cost function evolution with respect to the inner loop iterations.

Figure 4 .Figure 4 Figure 4 Figure 4

 4444 Figure 4.30 -Assimilation with real orthogonal-plane observations : comparison of of the norm of each component of the velocity field from t 1 = 20 to t 2 = 20.3.

  5/12 ,-17/60,0) (a 3 ,b 3 ,c 3 ) = (3/4 ,-5/12,0) RK4 (Kennedy et al., 2000) n k = 5 (a 1 ,b 1 ,c 1 ) = (0,0.1496590219993,0) and g 1 = 0.1496590219993 (a 2 ,b 2 ,c 2 ) = (-0.4178904745 ,0.1496590219993,0) and g 2 = 0.220741935365 (a 3 ,b 3 ,c 3 ) = (-1.192151694643 ,0.3792103129999,0) and g 3 = 0.25185480577 (a 4 ,b 4 ,c 4 ) = (-1.697784692471 ,0.8229550293869,0) and g 4 = 0.33602636754 (a 5 ,b 5 ,c 5 ) = (-1.514183444257 ,0.1530572479681,0) and g 5 = 0.041717869325 Annexe C

  Thus, at each step of the reverse mode computation, we require all the values of x p for p ∈ {1, ..., k}. There are two main strategies to tackle the reverse mode issue : Store all (SA) consists in saving all the values of {x p ; p = 1, • • • , k -1} and invoke each intermediate value when it is required. This technique is illustrated by figure C.1. While this strategy costs twice the computation time of the source program P, it clearly requires a tremendous amount of memory which can't be provided when we consider very large size problems. We can deal with this memory problem by a checkpoint strategy which consists in applying the SA strategy on some parts of the program and recompute from the beginning of the subprograms as illustrated by figure C.1.

Figure C. 1 -

 1 Figure C.1 -SA strategy (top) stores all the intermediate values (black points) at the end of each instructions {I k ; k = 1, • • • , p -1}. In the reverse mode, each values are invoked and used. Checkpointing SA strategy (bottom) consists in applying the SA strategy on small blocs of the program.

Figure

  Figure C.2 -RA strategy (top) performs all the intermediate instructions from the initial value x (black points) to obtain the required value. Checkpointing RA strategy (bottom) recomputes the intermediate values from each snapshot (black points)

[

  h, ū, v]= BC_east_b( hE , ūE , vE ) [ h, ū, v]= BC_west_b( hW , ūW , vW ) [ h, ū, v]= BC_north_b( hN , ūN , vN ) [ h, ū, v]= BC_south_b( hS , ūS , vS )D.2 Adjoint of the routines involved in the Roe solver

  of the eigenvalues of Ã(x R , x L )The eigenvalues e 1 , e 2 , e 3 of the linearized Jacobian Ã(x R , x L ) are computed from the input variable w := 1 2 (w L + w R ) such that d) + w1 , (D.3)where d = 1 for the horizontal fluxes and d = 0 for the vertical fluxes. In the following, we provide the tangent and adjoint procedures of the eigenvalues routine with respect to a given linearization state w := (w 1 , w 2 , w 3 ).The tangent procedure can be written into the following matrix form d -w 3 (1 -d) + w 1 w 1 w 1 -w 2 d -w 3 (1 -d) w 1 w 1 -w 2 d -w 3 (1 -d) + w 1 w 1 w 1   ẇ,and the corresponding adjoint procedure simply reads, d -w3 (1 -d) + w 1 -w 2 d -w 3 (1 -d) -w 2 d -w 3 (1 -d) + w 1 of the eigenvectors of Ã(x R , x L )

F

  of the nonlinear fluxWe recall that the nonlinear flux of a given state x = (h, u, v) T in the horizontal direction is given by F H := F (xrecall that the nonlinear flux in the vertical direction is given by

Figure G. 1 -

 1 Figure G.1 -Double frame configuration : the acquisition of the system 1 and system 2 is dt1 and dt2 respectively.

Figure

  Figure G.2 -Single frame configuration

Figure G. 4 -

 4 Figure G.4 -Dual plane position in the cylinder wake configuration at Reynolds 300

  Coupling CFD with EFD . . . . . . . . . . . . . . . . . . . . . . . . . Bayesian approach . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3 Best Linear Unibiased Estimator (BLUE) . . . . . . . . . . . 2.3 Sequential approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Kalman filter . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Extended Kalman Filter . . . . . . . . . . . . . . . . . . . . . 2.3.3 Ensemble Kalman Filter . . . . . . . . . . . . . . . . . . . . . 2.4 Variational approach . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Variational formulation . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Standard 4DVar technique . . . . . . . . . . . . . . . . . . . . 2.4.3 Incremental Variational Assimilation . . . . . . . . . . . . . .

Introduction Introduction 1 Turbulent 3D flow reconstruction 1.1 Computational Fluid Dynamics . . . . . . . . . . . . . . . . . . . . . 1.1.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Introduction to turbulence modeling . . . . . . . . . . . . . . 1.1.3 Numerical methods for CFD . . . . . . . . . . . . . . . . . . . 1.1.4 Limitations of computational fluid dynamics . . . . . . . . . . 1.2 Experimental fluid dynamics . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Particle Image Velocimetry for 2D observations . . . . . . . . 1.2.2 3D Particle Image Velocimetry . . . . . . . . . . . . . . . . . . 1.2.3 Limitation of Experimental Fluid Dynamics . . . . . . . . . . 1.3 1.3.1 Early attempts to combine measurements with a model . . . . 1.3.2 Towards data assimilation techniques . . . . . . . . . . . . . . 2 Data assimilation techniques 2.1 General context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Linear estimation theory . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Linear least squares estimation . . . . . . . . . . . . . . . . . 2.2.2 2.4.4 Inflow and initial conditions control . . . . . . . . . . . . . . . 2.4.5 Incremental variational assimilation with inflow and initial condition control . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.6 Numerical adjoint construction strategies . . . . . . . . . . . .

  Table 4.1 -List of the clusters used in this thesis framework of this work are summarized in table 4.1. The computational time and memory resources strongly depend on the size of the domain of study, therefore they will be precised in § 4.3 and § 4.4.

	Name CPU	REFERENCE RAM DISK	Network
	Calda Panda Gouda	2 x 8 cores Sandy Bridge 2 x 4 cores Clovertown 2 x 4 cores Clovertown	Intel(R) Xeon(R) CPU E5-2450 0 @ 2.10GHz Intel(R) Xeon(R) CPU L5310 0 @ 1.60GHz Intel(R) Xeon(R) @ 2.33GHz CPU E5345 0	48GB 2 x 600GB Infiniband + 1GB/s 8GB 1 x 73GB 1 GB/s 8GB 73GB 1 GB/s SAS 15k

Table 4 .

 4 2 -Tangent validation results for different values of the differentiation step α over 10∆t and 100∆t .

Table 4 .

 4 3 -4DVar validation result for different values of the differentiation step α in FORTRAN single and double precision.

  D.1.1 Contribution of the adjoint variable in the east direction

	Algorithm 37 BC_east_b( hE , ūE , vE )
	Define value within the domain
	for i=1,nx do
	for j=1,ny do
	h

  D'autre part, nous sommes capables de mesurer de manière non-intrusive le champ de vitesse d'un écoulement turbulent donné. En général, les techniques de mesures de Particle Image Velocimetry (PIV) sont utilisées pour mesurer les champs de vitesse en deux dimensions. Ces dernières années ont été marqué par un progrès technologique remarquable qui a permis la mesure des champs de vitesse en trois dimensions grâces à des techniques de mesure de tomographie-PIV. Cependant, les observations obtenues sont souvent bruitées et limitées par une résolution spatio-temporelle qui pourrait fausser une représentation réaliste de l'écoulement turbulent observé.Compte tenu des limitations de chaque approche, plusieurs méthodes ont été proposées pour combiner ces dernières afin de bénéficier des avantages et compenser les inconvénients de chacune d'entre elle. Dans le cadre de cette thèse, nous nous intéressons aux techniques d'assimilation de données qui combinent les modèles numériques avec les observations afin de déterminer une meilleure approximation du système d'intérêt. Dans le contexte de la mécanique des fluides, de telles méthodes ont été éprouvées dans des applications deux dimensionnelles, nous proposons dans cette thèse une première application à la reconstruction d'un écoulement turbulent en trois dimensions.Les techniques d'assimilation de données ont émergé depuis le domaine de la météorologie et se sont répandues dans de nombreuses applications géophysiques (océanographie, glaciologie etc…). Dans le contexte de l'assimilation de données, nous disposons d'un modèle décrivant le comportement du système considéré et permettant de calculer l'état de ce système à un temps donné. Nous possédons une première trajectoire, appelé l'ébauche, déduite d'une condition initiale a priori qui peut être le résultat d'une assimilation passée. Nous disposons d'autre part d'observations complètes ou partielles des caractéristiques du système d'intérêt. L'objectif des techniques d'assimilation de données est d'estimer l'état, appelé l'analyse, de ce système qui minimise à la fois son écart par rapport à l'ébauche et aux observations. Les techniques d'assimilation de données sont généralement regroupées en deux familles principales : les techniques séquentielles et les techniques variationelles. D'une part, comme son nom l'indique, l'approche séquentielle consiste à corriger l'ébauche à chaque fois que nous rencontrons une observation dans le temps. Cette approche est basée sur le filtre de Kalman et consiste à déterminer l'état d'analyse par le calcul d'une matrice de gain qui prend en compte les contributions passées. Cette approche a l'avantage d'évaluer l'incertitude sur l'analyse par le calcul de la matrice de covariance d'erreur d'analyse. En contrepartie, le calcul de la matrice de gain et la mise à jour de la matrice de covariance d'erreur d'analyse engendrent un coût en mémoire et temps de calcul considérables. En conséquence, de telles méthodes ne peuvent être appliquées dans le contexte de la mécanique des fluides. D'autre part, l'approche variationelle est basée sur les résultats de la théorie du contrôle optimal.

Acknowledgments

Annexe A BLUE

We seek for the gain matrix K involved in the analysis state x a given by (2.9) and recalled below

The background error is defined by η = x tx b associated with the background error covariance matrix B = E(ηη T ), (A.1) and the observation error ε = Y -Hx t associated with the observation error covariance matrix R = E(εε T ).

(A.

2)

The gain matrix minimizes the variance E(||ε a || 2 ) by using the information provided by the covariance matrices (A.1) and (A.2). The cost function of interest reads

and we cancel ints gradient with respect to K :

We then deduce the analysis covariance matrix v1, v1d, v2, v2d, v3, v3d)

Discrete adjoint of FFT_postproc

We describe the construction of the FFT_postproc routine which intervenes in the poisson routine (see algorithm 25. This routine perfoms a one dimensional post processing applied in every direction. In this section, we consider the construction in an arbitrary direction for an input a vector in of size n. The input variable sign indicates if we consider the forward post processing (sign = 1) or the backward post processing (sign = -1). The coefficients used in this routine are determined by the routine abxyz which doesn't involve any active variables, therefore we don't need to compute its discrete adjoint.

The original FFT_postproc routine is described by algorithm 43.

Algorithm 43 FFT_postproc(in,out,sign)

Initialization: call abxyz(ax,bx,ay,by,az,bz)

The discrete adjoint corresponding to the FFT_postproc routine is obtained straightforwardly by TAPENADE and is described by algorithm 44.

Annexe F. Discrete adjoint of FFT_postproc

Algorithm 44 FFT_postproc_b(in,in,out,out,sign) Initialization: call abxyz(ax,bx,ay,by,az,bz) Abstract : In the one hand, flow dynamics are usually described by the Navier-Stokes equations and the literature provides a wide range of techniques to solve such equations. On the other hand, we can nowadays measure different characteristics of a flow (velocity, pressure, temperature etc...) with non-intrusive Particle Image Velocimetry techniques. Within this thesis, we take interest in the data assimilation techniques, that combine a dynamics model with measurements to determine a better approximation of the system. This thesis focus on the classic variational assimilation technique (4DVar) which ensures a high accuracy of the solution by construction.

We carry out a first application of the 4DVar technique to reconstruct the characteristics (height and velocity field) of a uni directional wave at its free surface. The fluid evolution is simulated by the shallow water equations and solved numerically. We use a simple experimental setup envolving a depth sensor (Kinect sensor) to extract the free surface height. We compared the results of the 4DVar reconstruction with different versions of the hybrid data assimilation technique 4DEnVar.

Finally, we apply the 4DVar technique to reconstruct the downstream of a three dimensional cylinder wake at Reynolds 300. The turbulent flow is simulated by the high-performance multi-threading DNS code Incompact3d. This dynamics model is first combined with synthetic three dimensional observations, then with real orthogonal-plane stereo PIV observations to reconstruct the full three dimensional flow. Key words : variational assimilation, fluid dynamics, particle image velocimetry Résumé : D'une part, les équations de Navier-Stokes permettent de décrire les écoulements fluides, la littérature est riche de méthodes numériques permettant la résolution de celle-ci. D'autre part, nous sommes capables de mesurer de manière non-intrusive différentes caractéristique d'un écoulement (champ de vitesse et pression, etc.). Dans le cadre de cette thèse, nous nous intéressons aux techniques d'assimilation de données qui combinent les modèles numériques avec les observations afin de déterminer une meilleure approximation du système. Cette thèse s'articule autour de l'assimilation de donnée variationelle (4DVar) qui est plus précise par construction.

Nous avons mené une première application sur la reconstruction de la hauteur et vitesse de la surface libre d'un fluide contenu dans un récipient rectangulaire à fond plat. L'écoulement est modélisé par les équations de shallow water et résolues numériquement. Les observations de l'évolution de la hauteur de la surface libre ont été prélevées par un capteur de profondeur (Kinect). Nous avons comparé les résultats de la reconstruction par ls 4DVar avec plusieurs version de la méthode d'assimilation hybride 4DEnVar.

Enfin, nous avons appliqué la technique 4DVar à la reconstruction volumique de l'aval d'un sillage de cylindre à Reynolds 300. L'écoulement turbulent a été simulé par un code DNS parallèle Incompact3D. La reconstruction a été éffectué en combinant tout d'abord des observations synthètiques en trois dimension, puis en combinant des observations de plans orthogonales en stéréo PIV. Mots clefs : assimilation variationelle de données, dynamique des fluides, particle image velocimetry