Linh Van Nguyen

Jean-Philippe Laval

Pierre Chainais

A Bayesian fusion model for space-time reconstruction of finely resolved velocities in turbulent flows from low resolution measurements

The study of turbulent flows calls for measurements with high resolution both in space and time. We propose a new approach to reconstruct High-Temporal-High-Spatial resolution velocity fields by combining two sources of information that are well-resolved either in space or in time, the Low-Temporal-High-Spatial (LTHS) and the High-Temporal-Low-Spatial (HTLS) resolution measurements. In the framework of co-conception between sensing and data post-processing, this work extensively investigates a Bayesian reconstruction approach using a simulated database. A Bayesian fusion model is developed to solve the inverse problem of data reconstruction. The model uses a maximuma posteriori estimate, which yields the most probable field knowing the measurements. The DNS of a wall-bounded turbulent flow at moderate Reynolds number is used to validate and assess the performances of the present approach. Low resolution measurements are subsampled in time and space from the fully resolved data. Reconstructed velocities are compared to the reference DNS to estimate the reconstruction errors. The model is compared to other conventional methods such as Linear Stochastic Estimation and cubic spline interpolation. Results show the superior accuracy of the proposed method in all configurations. Further investigations of model performances on various range of scales demonstrate its robustness. Numerical experiments also permit to estimate the expected maximum information level corresponding to limitations of experimental instruments.

Introduction

Turbulence, though governed by Navier-Stokes equations, is extremely hard to predict due to its spatiotemporally intermittency as well as three-dimensional and irregular properties. It is also a multi-scale phenomenon where a very wide range of scales from the largest eddies to Kolmogorov micro-scales co-exist and interact. Since the ratio between the largest and the smallest scales increases with Reynolds number as Re 3/4 , flows with high Reynolds are the most challenging. Wall bounded flows are particularly difficult to model due to the overlap of several scaling regions as a function of distance to the wall. Coherent structures in such flows can extend up to several boundary layer thickness. The modeling of such structures and scales therefore requires extremely detailed flow information in both space and time.

Despite a constant progress, none of the experimental techniques, even in academic researches, is capable of providing spatiotemporally resolved information in sufficiently wide spatial domains and for diverse flow conditions. Particle Image Velocimetry (PIV), the most advanced turbulence measurement technique, cannot measure space-time resolved velocities. Stereoscopic PIV measures three-component velocities at high spatial resolution and large field-of-view, but limited to a low acquisition rate compared to the flow dynamics. High repetition tomographic PIV and Time-resolved PIV (TrPIV) are improving but still limited to small volumes and low speed flows. Other point-measurement techniques such as Hot Wire Anemometry (HWA) measure the full temporal dynamics. However, the combination of these devices to get a better spatial resolution is not straightforward and remains intrusive.

Direct Numerical Simulation (DNS) can provide reliable and fully resolved velocities of turbulent flows. It simulates the flows by directly solving Navier-Stokes equations. The computational cost of such a numerical Figure 1: Sketch of the inverse problem, with the two sources of measurements: the LTHS (color images) and a coarse grid of HTLS (red dots among black ones of LTHS). The inverse problem of HTHS data reconstruction is to fill in the space-time data-cube.

approach is very high since the number of simulated grid points increases as Re 9/4 . DNS therefore can simulate flows with low to moderate Reynolds and simple geometries only.

To have fully resolved velocities, one idea is to measure and combine two types of complementary measurements in space and time: the Low-Temporal-High-Spatial resolution (LTHS) and the High-Temporal-Low-Spatial resolution (HTLS) measurements. One particular example of such an idea is presented in [START_REF] Coudert | Double large field stereoscopic piv in a high reynolds number turbulent boundary layer[END_REF]. This joint experiment provides a database of high Reynolds boundary layer flows. The data are provided by a stereoscopic PIV synchronized with a rake of HWA probes. PIV is with a large field-of-view and at a high spatial resolution but at low acquisition frequency. HWA measurements are at an extremely high temporal resolution, but the spatial discretization of the rake of probes is very coarse compared to Kolmogorov scales.

Various methods have been proposed to combine such measured data of turbulent flows to recover the maximum information level. Linear Stochastic Estimation (LSE) is the most common one. Its introduction into turbulence community dates back to the works by Adrian [START_REF] Adrian | On the role of conditional averages in turbulence theory[END_REF][START_REF] Adrian | Conditional eddies in isotropic turbulence[END_REF] and has been further investigated later [START_REF] Guezennec | Stochastic estimation of coherent structures in turbulent boundary layers[END_REF][START_REF] Adrian | Stochastic estimation of conditional structure: a review[END_REF][START_REF] Ewing | Examination of a lse/pod complementary technique using single and multitime information in the axisymmetric shear layer[END_REF]. These works use LSE as a tool to extract coherent structures from the measurements. Later works proposed various extensions such as multi-time, nonlinear or higher-order LSE [START_REF] Mokhasi | Predictive flow-field estimation[END_REF][START_REF] Durgesh | Multi-time-delay lse-pod complementary approach applied to unsteady high-reynolds-number near wake flow[END_REF][START_REF] Nguyen | Proper orthogonal decomposition-based estimations of the flow field from particle image velocimetry wall-gradient measurements in the backwardfacing step flow[END_REF][START_REF] Meyer | Provide detailed and real-time indoor environmental information using pod-lse and limited measurements[END_REF]. In these works, unknown velocities are reconstructed from measurements of other quantities such as pressure or shear-stress. LSE can be also linked to Proper Orthogonal Decomposition (POD) to reduce the order of reconstruction problems [START_REF] Bonnet | Stochastic estimation and proper orthogonal decomposition: complementary techniques for identifying structure[END_REF].

The idea of combining sparse velocity measurements to obtain fully-resolved fields has not been addressed until recently [START_REF] Melnick | Experimental investigation of a turbulent boundary layer using simultaneous 3-d flow visualization and 2-d piv[END_REF][START_REF] Tu | Integration of nontime-resolved piv and time-resolved velocity point sensors for dynamic estimation of velocity fields[END_REF]. In [START_REF] Melnick | Experimental investigation of a turbulent boundary layer using simultaneous 3-d flow visualization and 2-d piv[END_REF], 3D smoke intensity and 2D PIV measurements are combined using a POD-LSE model to get fully resolved 3D velocities of a flow over a flat plate. POD-LSE estimation model has been developed further [START_REF] Tu | Integration of nontime-resolved piv and time-resolved velocity point sensors for dynamic estimation of velocity fields[END_REF] with a reconstruction scheme based on a multi-time LSE. Either a Kalman filter or a Kalman smoother is used depending on the problem as real time estimation or data post processing. The model is tested using TrPIV measurements of a bluff-body wake at a low Reynolds number. Sparse velocity measurements are virtually extracted from the high resolution ones, while original data are used to estimate reconstruction errors.

LSE suffers from critical limitations though extensively used. First, as a conditional average, LSE estimates a set of coefficients that associate the so-called conditional eddies to one flow pattern [START_REF] Adrian | Conditional eddies in isotropic turbulence[END_REF]. Using these coefficients to reconstruct all velocity fields, LSE fails to capture coherent structures and misleads physical interpretations when particular patterns exist. Second, the reconstructed structures are independent of event magnitudes [START_REF] Moser | Statistical analysis of near-wall structures in turbulent channel flow[END_REF]. Reconstructed flows are associated with weak fluctuations only. Last, LSE as a low pass filter reconstructs large scales only and lose flow details even at measured positions.

u 2 Σ = u T Σ -1 u u 2 2 square of a Euclidean distance: u 2 2 = u T u n ∼ N(0, Σ n ) Gaussian noise vector of zero-mean and covariance matrix Σ n N(u|µ u , Σ n )
Gaussian distribution of variable u that takes the mean µ u and fluctuates due to n of covariance Σ n N(u|v) distribution of u knowing (or conditioning on) v

The present work proposes a novel model to reconstruct the fully resolved HTHS velocities from HTLS and LTHS measurements. This model is based on a Bayesian inference framework using a Maximum A Posteriori (MAP) estimate [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]. It is inspired by the multispectral image fusion problem with the limited resolution of image measurements in space-wavelength domains [START_REF] Hardie | Map estimation for hyperspectral image resolution enhancement using an auxiliary sensor[END_REF]. This framework has been discussed early in communication problems [START_REF] Van Trees | Detection, estimation, and modulation theory, optimum array processing[END_REF][START_REF] Sage | Estimation theory with applications to communications and control[END_REF] and is used more extensively in image processing, remote-sensing, and data fusion [START_REF] Levitan | A maximum a posteriori probability expectation maximization algorithm for image reconstruction in emission tomography[END_REF][START_REF] Ma | Simultaneous detection of lane and pavement boundaries using model-based multisensor fusion[END_REF][START_REF] Challa | Bayesian and dempster-shafer fusion[END_REF][START_REF] Koks | An introduction to bayesian and dempster-shafer data fusion[END_REF][START_REF] Durrant-Whyte | Multisensor data fusion[END_REF][START_REF] Zhang | Noise-resistant wavelet-based bayesian fusion of multispectral and hyperspectral images[END_REF][START_REF] Joshi | Map estimation for multiresolution fusion in remotely sensed images using an igmrf prior model[END_REF]. The Bayesian fusion model takes benefit from both sources of information in space and time simultaneously by searching for the most probable flow for given measurements. Better performances are expected since space and time correlations are equally important. The model also recovers flow details inaccessible from single interpolations. By integrating directly the measurements, it proposes a compromise estimate such that detailed flow information close to the sensor positions are well preserved. This approach also overcomes the limitations of LSE, which acts as a low pass filter due to the mean square error minimization. To test the model, the DNS database of a turbulent wall-bounded flow is used. These space-time fully resolved data allow the model optimization and validation. Sparse measurements of HTLS and LTHS are extracted from the full dataset, while the reference DNS data are used in the end to evaluate reconstruction errors. Performances are evaluated for various configurations with different subsampling ratios.

The paper is organized as follows. Section 2 presents the Bayesian model using a MAP estimate. Model simplification and statistical parameters estimation are also discussed. Section 3 describes the DNS database used to test the model and also other reconstruction methods for comparison. Results for various configurations are presented. Conclusions and future works are in Section 4.

Bayesian fusion model 2.1 Bayesian model

Let x and y denote LTHS and HTLS measurements, and z denote HTHS data to reconstruct. z, x and y are random, zero-mean vectors of size NP × 1, NQ × 1 and MP × 1 respectively. N and M are numbers of spatial points in each snapshot, while P and Q are numbers of snapshots. The present work is a challenging inverse problem since we consider M N and Q P. Let the subscript "s" denote operators performing in space, and "t" be those in time; is an interpolator; is for subsampling; is a Low-Pass Filter (LPF). The cubic spline interpolation either 1D or 2D is used as for its state-of-the-art interpolation results and finite support [START_REF] Unser | Splines: A perfect fit for signal and image processing[END_REF][START_REF] Thevenaz | Interpolation revisited [medical images application[END_REF]. is a 5th-order least-square spline filter [START_REF] Reinsch | Smoothing by spline functions[END_REF][START_REF] Cook | The smoothing spline: a new approach to standardizing forest interior treering width series for dendroclimatic studies[END_REF] for its sharp cutoff response to better separate large scales from small scales. Table 1 lists all notations used in this paper.

The direct model of the measurement system involves the subsampling operator and some measurement noises:

x = t z + b t (1)
y = s z + b s (2)
where b t and b s are the (typically white Gaussian) measurement noises. Therefore, given sparse measurements of either y in space or x in time, two estimators ẑ1 and ẑ2 of the fully resolved vector z can be reconstructed by single interpolations. The 1D time interpolation goes from NQ to NP dimensional space, i.e. x -→ ẑ1 = t x , while the 2D space interpolation goes from MP to NP dimensional space, i.e. y -→ ẑ2 = s y . Let NP dimensional vectors h s and h t denote the information that cannot be recovered by simple interpolations; z can be modeled in two ways from these separate estimates:

z = ẑ1 + h t = t x + h t (3)
z = ẑ2 + h s = s y + h s (4)
Missing information h t and h s essentially feature small scales. Using either x or y, it is not possible to estimate h t and h s . The idea of Bayesian fusion is to combine the two models by using t x in (3) to estimate the unknown h s in (4) and vice-versa.

Let N(u|µ u , Σ u ) denote the multivariate Gaussian distribution of a NP dimensional random vector u with mean value µ u and covariance matrix Σ u . The NP × NP matrix is the expectation of (u -µ u )(u -µ u ) T . The probability density function (pdf) of u with a multivariate Gaussian distribution N(u|µ u , Σ u ) is:

p(u) = 1 (2π) NP/2 |Σ u | 1/2 exp - 1 2 u -µ u 2 Σ u (5)
where |.| denotes the matrix determinant, and uµ u 2 Σ u is the Mahalanobis distance:

u -µ u 2 Σ n = u -µ u T Σ -1 n u -µ u (6)
Let assume that t x and h t are approximately independent; t x captures temporal large scales of x. Similarly, s y and h s are assumed to be approximately independent. Due to subsampling, aliasing terms are also present in each pairs of ( t x, h t ) and ( s y, h s ). Assume also that h t and h s are zero mean Gaussian noises, i.e. h t ∼ N(0, Σ h t ) and h s ∼ N(0, Σ h t ). Pdfs of these unknowns are modeled as:

p(h t ) = 1 (2π) NP/2 |Σ h t | 1/2 exp - 1 2 h t 2 Σ h t (7)
and similarly for p(h s ). Posterior distributions of z knowing either x or y are then modeled as:

N(z|x) ∼ N(z| t x, Σ h t ) (8) N(z|y) ∼ N(z| s y, Σ h s ) (9)
where N(z|x), resp. N(z|x), is the posterior distribution of z knowing x, resp. knowing y.

MAP estimation

The present Bayesian model aims to build an estimate of z given x and y using the probability models ( 8) and [START_REF] Nguyen | Proper orthogonal decomposition-based estimations of the flow field from particle image velocimetry wall-gradient measurements in the backwardfacing step flow[END_REF]. The model uses a MAP estimate to search the most probable ẑ given x and y such that ẑ maximizes the posterior pdf p(z|x, y):

ẑ = argmax z p (z|x, y) (10)
Using Bayesian rules [START_REF] Kendall | Kendall's Advanced Theory of Statistics[END_REF], one has:

p (z|x, y) ∝ p (x, y|z) p(z) (11)
Assuming that x and y are independent conditioned on z, equation ( 11) becomes:

p (z|x, y) ∝ p(x|z)p(y|z)p(z) (12)
In equation ( 12), the likelihood functions p(x|z), p(y|z) and the prior pdf p(z) appear, while only the posterior probabilities p(z|x) and p(z|y) are available in ( 8) and ( 9).

To complete the model, the likelihood functions can be expressed in term of the posterior pdfs and prior of z using Bayesian rules. Section 2.3.3 in [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF] introduces an alternative way to estimate these functions from posterior pdfs using a linear Gaussian model. Various tests of different Gaussian priors p(z) lead to the use of a noninformative prior. This prior, referred also as vague or flat prior, assumes that all the values of z are equally likely [START_REF] Gelman | Bayesian data analysis[END_REF]. The estimation of ẑ is now solely based on the measurements and not influenced by external information. The prior distribution therefore has no influence on the posterior pdfs.

With the assumption of a noninformative prior, p(z) is constant. Using Bayes rules, the relation between the likelihood function and the posterior pdf is:

p (z|x) ∝ p (x|z) p(z) (13)
Since p(z) is replaced by a constant, one gets p (x|z) ∝ p (z|x). Similarly, p (y|z) ∝ p (z|y). Eq [START_REF] Melnick | Experimental investigation of a turbulent boundary layer using simultaneous 3-d flow visualization and 2-d piv[END_REF] becomes:

p (z|x, y) ∝ p(z|x)p(z|y) (14)
The MAP estimation is:

ẑ = argmax z p(z|x)p(z|y) (15)
Logarithms of p (z|x) and p (z|y) are:

-ln p(z|x) = 1 2 z -t x 2 Σ h t + C 1 ( 16
)
-ln p(z|y) = 1 2 z -s y 2 Σ hs + C 2 (17)
where C 1 and C 2 are independent of x, y and z. Solving ( 15) is equivalent to minimize the cost function:

C(z) = 1 2 z -t x 2 Σ h t + 1 2 z -s y 2 Σ hs (18)
Computing the gradient of C(z) and setting to zero:

∂C(z) ∂z = Σ -1 h s (z -s y) + Σ -1 h t (z -t x) = 0 ( 19
)
the solution to the optimization problem [START_REF] Meyer | Provide detailed and real-time indoor environmental information using pod-lse and limited measurements[END_REF] is:

ẑ = Σ -1 h t + Σ -1 h s -1 Σ -1 h s s y + Σ -1 h t t x (20)
Applying the matrix inversion lemma [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]:

(A + BD -1 C) -1 = A -1 -A -1 B(D + CA -1 B) -1 CA -1 (21)
Equation ( 20) can be rewritten as:

ẑ = Σ h t + Σ h s -1 Σ h t s y + Σ h s t x (22)
Equation ( 22) is the final full form of the proposed Bayesian fusion model using a MAP estimate and assuming a noninformative prior of z. Variance matrices Σ h t and Σ h s are parameters to be estimated.

Model simplification

Though providing the full theoretical estimate of ẑ, equation ( 22) is impractical to use as is for several reasons. The full covariance matrices Σ h t and Σ h s , representing all the sources of correlations in space and time, cannot be estimated from only the measurements x and y. This is because the unknown h t and h s are only accessible at the measured positions in space and time. Also, the covariance matrices of size NP × NP are very large, making them very difficult to accurately estimate and to inverse. Additional assumptions on the shape of Σ h t and Σ h s are necessary.

A common and simple approach is to assume diagonal covariance matrices. This implies the independence between all elements of h t and h s . The simplified version of equation ( 22) becomes a point-wise formula:

ẑ(i) = σ 2 h s (i) σ 2 h s (i) + σ 2 h t (i) t x(i) + σ 2 h t (i) σ 2 h s (i) + σ 2 h t (i) s y(i) (23)
where i(t, s) is the index of each point in time (t) and space (s). The variances σ 2 h s and σ 2 h t are functions of each position in space and time. Their estimation is detailed in next section. Then equation [START_REF] Durrant-Whyte | Multisensor data fusion[END_REF] will be used to reconstruct HTHS data. As a weighted average, it proposes a compromise estimate from the measurements. With a symmetrical form in space and time, the model uses information from both measurements to correct large scales reconstruction and recover certain information at smaller scales.

Statistical parameters estimation

Let Z, X, Y,H t , H s , Γ h t and Γ h s be the (time,space) matrix forms of z, x, y, h t , h s , σ 2 h t and σ 2 h s respectively. Z, H t , H s , Γ h t and Γ h s are of size P × N, while X and Y are of size Q × N and P × M. Γ h t and Γ h s are matrices of empirical variances, which are functions of time and space (t, s).

Variance matrices are estimated from H t and H s , which are available at the measurements positions only. We use:

t H s = X -s s X ( 24
)
s H t = Y -t t Y (25)
where t subsamples in time from P to Q time steps, and s subsamples in space from N to M points. These Q instants and M positions are the same as for LTHS and HTLS measurements. Since the flow is approximately stationary and spatial interpolation is independent of time, Γ h s (t, s) becomes Γ h s (s), a function of spatial locations only. These variances are estimated by averaging over all time steps:

Γ h s (s) = 1 Q Q t=1 (X(t, s) -s s X(t, s)) 2 (26)
Variance in Γ h t is a function of distances τ to the previous LTHS time step only, where τ/δt = 0, 1, 2, ..., P/Q, and δt is the time lag between two consecutive HTHS time steps. Γ h t becomes a function of space and τ only, i.e. Γ h t (τ, s). It is estimated by averaging over Q blocks (of P/Q snapshots) bounded by two consecutive LTHS instants:

Γ h t (τ, s) = 1 Q t s (Y(t s , s) -t t Y(t s , s)) 2 (27)
where t s /δt = τ/δt, τ/δt+P/Q, τ/δt+2P/Q, ..., τ/δt+(Q-1)P/Q. Since the flow is approximately homogeneous in spanwise direction, Γ h s (s) and Γ h t (τ, s) are also averaged over all blocks defined by the four neighboring HTLS measurements, see Fig. 2. The variances are then functions of only vertical positions and relative distances to the four closest HTLS sensors. These estimated variances are rearranged into a vector form σ 2 h t (i) and σ 2 h s (i) to complete the fusion model using the simplified formula in equation ( 23).

Numerical experiments

Section 3.1 describes the DNS database used to test the model. Section 3.2 discusses other reconstruction methods for comparison. Section 3.3 presents results of the fusion model in various cases.

DNS database

DNS database of a turbulent wall-bounded flow is used to test the model. This simulation uses the numerical procedure described in [START_REF] Marquillie | Direct numerical simulation of a separated channel flow with a smooth profile[END_REF]. The flow is at a Reynolds number Re τ = 550 based on the friction velocity. Cartesian coordinates of the simulation in space are (x, y, z) for streamwise, vertical and spanwise directions respectively. The domain size L x × L y × L z normalized by half the channel height H is 2π × 2 × π. Fully resolved fluctuating streamwise velocities in a plane normal to the flow direction are considered as HTHS data. This data includes P = 10000 snapshots at spatial resolution of N = 288 × 257 and at sampling frequency of 40 Hz. Sparse LTHS and HTLS measurements are subsampled from HTHS data to learn the fusion model. HTHS is used as the ground truth to estimate reconstruction errors. The extension to spanwise and vertical velocity components follows the same procedure.

Various cases are investigated. The subsampling ratios N/M applied in each direction of space are 5, 10 and 20. These ratios correspond to a number M of HTLS sensors of 51 × 57, 26 × 29 and 13 × 15 respectively. Each ratio has a spacing between two successive HTLS points in spanwise and vertical directions of ∆z and ∆y. Subsampling ratios P/Q in time are 4 (Q = 2500), 10 (Q = 1000) and 20 (Q = 500). Each ratio, both in space and time, corresponds to a certain amount of energy loss. This is essentially the energy of small scales separated from large scales by a low pass filter . Here is the 5 th -order least square spline filter, either temporal 1D ( t ) or spatial 2D ( s ), using measurements as knots. This spline filter has the advantages of a sharp cutoff response and finite support. The energy loss is defined by comparing the filtered field z and the original field z:

∆κ = j∈ z 2 j - j∈ [ z] 2 j j∈ z 2 j ( 28
)
where is the considered set of points. Table 2 gathers the energy loss in time (∆κ t ) and in space (∆κ s ) estimated with t and s respectively. The set contains all points at y/H = 1.

Other methods for comparison

Other reconstruction methods are used for comparison with the present model. Cubic spline interpolation: Interpolation techniques reconstruct HTHS velocities from either LTHS or HTLS measurements independently, i.e. x -→ ẑ = t x or y -→ ẑ = s y. The cubic spline interpolations [START_REF] Unser | Splines: A perfect fit for signal and image processing[END_REF], either 1D in time or 2D in space, are used. These interpolations are by Matlab built-in functions, which follow the algorithm in [START_REF] Boor | A practical guide to splines[END_REF].

Linear Stochastic Estimation: LSE estimates ẑ as a linear combination of measurements. Coefficients are estimated from the measurements by solving a system of linear equations to minimize the mean square errors of reconstructed fields. Refs. [START_REF] Adrian | Conditional eddies in isotropic turbulence[END_REF][START_REF] Adrian | Stochastic estimation of conditional structure: a review[END_REF] describe the physical interpretations of this procedure. This section derives the model differently [START_REF] Hastie | The elements of statistical learning[END_REF][START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF] but in accordance with turbulence literature.

Matrix forms X, Y and Z described in Section 2.4 are used to build the LSE model. Let Y s = t Y of size Q × M denote a part of Y subsampled at the same instants as X. LSE model finds the optimal matrix B of size N × M that minimizes the residual sum of squared errors:

B = argmin B Y s B -X 2 2 (29)
Let set the gradient of this residual sum to zero:

∂ Y s B -X 2 2 ∂B = Y T s (Y s B -X) = 0 ( 30
)
the optimal B is obtained as:

B = Y T s Y s -1 Y T s X (31)
Equation ( 31) requires the inversion of (Y T s Y s ) that can be singular, leading to a high variance model with large coefficients. A small change of predictors Y then can lead to a very different reconstruction of Z, causing model's instability. Tikhonov regularization [START_REF] Cordier | Calibration of pod reduced-order models using tikhonov regularization[END_REF], well-known in machine learning problems as L2 penalty or Ridge Regression [START_REF] Hastie | The elements of statistical learning[END_REF][START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF], can be used as a remedy. It aims to solve this ill-posed problem by imposing a L2 penalty term on the residual sum of errors. The optimization problem [START_REF] Cook | The smoothing spline: a new approach to standardizing forest interior treering width series for dendroclimatic studies[END_REF] becomes:

B = argmin B Y s B -X 2 2 + λ B 2 2 ( 32
)
Setting the gradient of the cost function (for λ > 0) to zero:

∂ Y s B -X 2 2 + λ B 2 2 ∂B = Y T s (Y s B -X) + λB = 0 ( 33
)
the closed form of B is:

B = Y T s Y s + λI -1 Y T s X (34)
The regularization parameter λ can be optimized by ten-fold cross-validation [START_REF] Geisser | Predictive inference[END_REF]. The fully resolved field of Z is then estimated using these coefficients:

Z = Y B ( 35
)
Matrix B encodes the predictor of Z knowing Y learnt from the joint observation of X and Y. A completely analogous procedure can be used switching the roles of X and Y.

Table 2: NRMSEs of all scales reconstruction errors for various cases. The subsampling ratios of HTLS measurements are N/M and equal in both spatial directions. The ratios of LTHS measurements in time are P/Q. The equivalent spacing in spanwise direction is normalized by half channel height as ∆z/H and the spacing in time is ∆t. The normalized energy losses in space ∆κ s and in time ∆κ t are defined in Eq. [START_REF] Reinsch | Smoothing by spline functions[END_REF]. and max are the mean and max NRMSE defined in Eq. [START_REF] Cordier | Calibration of pod reduced-order models using tikhonov regularization[END_REF]. is averaged over all space-time positions in the outer region y/H ∈ [0.25, 1.75], while max is computed for one of the most difficult position in space and time (the most remote from all nearby measurements). The smallest errors in each cases are boldfaced.

Subsampling ratios

Impact of subsampling ratios

The fusion model uses equation [START_REF] Durrant-Whyte | Multisensor data fusion[END_REF] to reconstruct fully resolved velocities ẑ in various cases. Reconstructed fields are compared with the original DNS via the Normalized Root Mean Square Error (NRMSE):

NRMS E =             j∈ (ẑ j -z j ) 2 j∈ z 2 j             1/2 (36)
where is the considered set of points used to estimate the error. The field z is more or less difficult to estimate depending on the considered instant and position with respect to available measurements. To qualify, two types of NRMSE, the mean NRMSE and the maximum NRMSE max , are estimated. is estimated over including all space-time positions in the outer region of y/H ∈ [0.25, 1.75], where the flow is approximately homogeneous. It represents how far the reconstructed field departs from ground truth in order to evaluate reconstruction accuracy. max is estimated using all blocks (in time and in spanwise directions) bounded by HTLS sensors at y/H = 0.94 and y/H = 1.06, see Fig. 2. The set includes centers at local coordinates (∆y/2, ∆z/2, Pδt/2Q) of all blocks. and max of s y, t x and LSE reconstruction are also estimated for comparison.

Table 2 describes 7 cases with their settings and reconstruction errors. In cases 1 and 2, the energy losses due to subsampling in time are much higher than in space, and vice-versa in cases 3 and 4. The model gives similar errors compared to the best interpolation, with smaller and comparable max . In cases 5 to 7, the losses are due to both the subsamplings in space and time in a balanced manner. The proposed model reduces by 15% to 30% and max by 10% to 20% compared to the best of other methods.

Improvements are expected from the weighted average in equation [START_REF] Durrant-Whyte | Multisensor data fusion[END_REF]. The present model uses variances σ 2 s (i) and σ 2 t (i) as parameters of the flow's physics, and t x and s y as the specific flow information. It imposes the reconstruction to be consistent with measurements at nearby positions and proposes compromise estimates elsewhere. Simple interpolations use either HTLS or LTHS measurements only, losing information from the other source. LSE learns its coefficients from both measurements but inherits the limitations of the conditional averaging.

Large and small scales reconstruction

In cases 1 to 4, the fusion model performs as the best interpolation with small improvements. This is expected since one measurement of HTLS or LTHS is much better resolved than the other. Cases 5 to 7 are the most interesting since energy losses due to subsampling in space and time are comparable. The model brings complementary information from both measurements and improves the reconstruction.

We study reconstructions of large and small scales in details for these three cases. Spatial 2D filters s (see Section 3.1) are used to separate large scales from small scales. These filters take HTLS points as knots to have a cutoff close to the Nyquist frequency. The reconstructed large scales by all methods are compared to the reference s z. Small scales are estimated using Is where I is the identity matrix. Table 3 shows NRMSEs estimated using equation [START_REF] Cordier | Calibration of pod reduced-order models using tikhonov regularization[END_REF] but normalized by the RMS of either s z or (Is )z.

The fusion model recovers part of small scales from complementary measurements. It gives the lowest and max of small scales reconstruction in all cases. It also better reconstructs large scales than other methods. For large scales, max remains the same in case 5 of small subsampling ratios and improves significantly in cases 6 and 7 of high ratios, with max reduced by 5 % and 25 % respectively, and by 20% to 40% compared to the best of other methods.

Model performance analysis

We focus on case 6 for a model performance analysis. This case has about 5 % energy losses due to both time and space subsamplings, which are critical to highlight interests of the present approach. The model reduces and max by 25 % and 35 % respectively for all scales reconstruction, 10 % and 5 % for large scales reconstruction. To analyze reconstructions in space, Fig. 3(a) shows spatial NRMSE maps by all methods as functions of local coordinates (α, β). For each (α, β), NRMSE is estimated using equation [START_REF] Cordier | Calibration of pod reduced-order models using tikhonov regularization[END_REF], where includes points at (α, β, Pδt/2Q) of all blocks used to estimate max (see Section 3.3.1). For all methods, NRMSEs are small close to the four HTLS positions in the corners and increase when approaching the center. Time interpolation behaves differently since its errors are independent of spatial coordinates. The fusion model yields the smallest errors at all positions. It improves significantly near the center compared to spatial interpolation, the best of other methods.

To analyze reconstructions in time, Fig. 3(b) shows the NRMSE curves by all methods as functions of distances τ from the previous LTHS time step. For each τ, NRMSE is estimated using including points at local coordinates (∆y/2, ∆z/2, τ) of all blocks used to estimate max . NRMSEs are small close to the LTHS measurements (τ/δt = 0 and τ/δt = P/Q) and increase when moving toward the middle (τ/δt = P/2Q). Spatial interpolation are different with NRMSEs independent of time. The fusion model yields the minimum errors at all time steps. Even in the middle of two LTHS instants, the maximum fusion error remains significantly lower than that of all other methods. Fig. 4 shows a time evolution of the point at y/H = 1 and z/H = 0 (α = ∆y/2 and β = ∆z/2 in local coordinates), the most remote from its neighboring HTLS sensors. A good agreement between fused and reference velocity is still obtained. A zoom-in period is shown also for detailed comparisons with other methods. While time interpolation captures only low frequencies, spatial interpolation generates high frequencies but weakly correlated with the truth. The fusion model proposes a good compromise to improve both large and small scales reconstruction. It also captures detailed peaks much better than LSE, since LSE smooths these small scales out by minimizing the mean square errors.

Fig. 5 compares temporal spectra of above evolutions. Time interpolation fails to estimate the signal at higher frequencies than a certain cutoff. LSE keeps both large and small scales, but the loss of large scale energy is critical. This loss is highlighted in the zoom-in picture of low frequencies spectral. The present model improves the estimation at both low and high frequencies.

Fig. 6 shows estimates of the probability density function estimates of time increments u(x, t+τ)-u(x, t) for the original DNS field as well as for the reconstructed field. As expected, the original field displays intermittent non Gaussian distributions. More importantly, the reconstructed field, while less intermittent, still clearly exhibits non Gaussian increments at small scales. Note that the reconstruction error is essentially due to the difficulty to accurately reconstruct these small scales. It is expected that any reconstruction method will lead to fields that are less intermittent than the original one. Fig. 7 compares reconstructed snapshots by different methods. This snapshot is at the most remote instant from its two neighboring LTHS time steps. The model reconstructs correctly the velocity field with more flow details than spatial interpolation. It also recovers better large scales than LSE and time interpolation methods.

Conclusions

This work proposes a Bayesian fusion model using a MAP estimate to reconstruct high resolution velocities of a turbulent channel flow from low resolution measurements in space and time. It searches for the most probable field given available measurements. This approach yields a simple but efficient weighted average formula in equation [START_REF] Durrant-Whyte | Multisensor data fusion[END_REF]. Weighting coefficients are learnt from measurements and encode the physics of the flow. The informed fusion of information from available measurements improves the interpolation of large scales and recovers details at small scales.

Numerical experiments using a DNS database of a turbulent wall-bounded flow at a moderate Reynolds number illustrate the efficiency and robustness of the proposed method. Low resolution measurements are extracted to learn model parameters, while original data are used as the ground truth to estimate reconstruction errors. The model is tested in various cases with different subsampling ratios. Results are compared to more standard methods such as cubic spline interpolation and penalized LSE. Bayesian fusion always produces the most accurate reconstruction. The best results are obtained when missing spatial and temporal information are of the same order of magnitude. In these cases, it provides a better large scale reconstruction while a certain amount of small scale details are also recovered. The search for an even more accurate fusion and super-resolution method is the subject of ongoing work.

Figure 2 :

 2 Figure 2: Sketch of an element block with local coordinates (α, β, τ). LTHS time steps are at τ/δt = 0 and τ/δt = P/Q. HTLS measurements are represented by red dots and LTHS measurements by black ones.

Figure 3 :

 3 Figure 3: NRMSEs between reference and reconstructed streamwise velocities by all methods as: (a) functions of spatial coordinates in an element block at the most difficult instant, i.e. at (α, β, Pδt/2Q); (b) functions of time distances from the previous LTHS instant at the most difficult spatial location, i.e. at (∆y/2, ∆z/2, τ).

Figure 4 :

 4 Figure 4: A time evolution of fluctuating streamwise velocity at y/H = 1 and z/H = 0, the centers of all such (α, β) planes in Fig. (2).

Figure 5 :

 5 Figure 5: Spectra of the fluctuating velocity in Fig. 4.

Figure 6 :

 6 Figure 6: Probability distribution functions of velocity increments in (left) the original DNS, (right) the reconstructed field.

Figure 7 :

 7 Figure 7: A sample snapshot of fluctuating streamwise velocity at one of the most difficult instant to estimate (in the middle of two LTHS time steps): Reconstruction of all scales (left) and large scales only (right). The figure is better viewed on screen.

Table 1 :

 1 Summary of notations.

	z	NP-dimensional vector of HTHS DNS data
	y	MP-dimensional vector of HTLS measurements
	x	NQ-dimensional vector of LTHS measurements
	N	number of spatial points in each HTHS/LTHS snapshot
	M	number of spatial points in each HTLS snapshot
	P	number of HTHS or HTLS snapshots
	Q	number of LTHS snapshots
	s	2D 5th-order least-square spline filter in space
	t	1D 5th-order least-square spline filter in time
	∆κ	loss of kinetic energy computed using Eq. (28)
		normalized Root Mean Square Error (NRMSE) defined in Eq. (36)
	|.|	determinant of a matrix
	u T	transpose of vector u
	argmax (.)	argument of the maximum, which is z for which the function attains its
	z	maximum
	u 2 Σ	square of a Mahalanobis distance:

t 1D cubic spline interpolator in time s 2D cubic spline interpolator in space t subsampling in time from P to Q snapshots, t z = x s subsampling in space from N to M points, s z = y

Table 3 :

 3 NRMSEs of large and small scales. Notations are explained in Table.2.

	Spacings	Energy loss


Damien Do Couto
email: damien.docouto@unige.ch

Speranta-Maria Popescu
email: speranta.popescu@gmail.com

Jean-Pierre Suc
email: jeanpierre.suc@gmail.com

Mihaela Carmen Melinte-Dobrinescu

Nadia Barhoun
email: nbarhoun@yahoo.fr

Christian Gorini
email: christian.gorini@gmail.com

Laurent Jolivet
email: laurent.jolivet@univ-orleans.fr

Jeffrey Poort
email: jeffrey.poort@upmc.fr

Gwénaël Jouannic
email: gwenael.jouannic@cerema.fr

Jean-Luc Auxietre
email: jean-luc.auxietre@total.com

Reply to the comment on the paper "Lago Mare and the Messinian Salinity Crisis: Evidence from the Alboran Sea (S. Spain) by Do Couto et al. (2014) Marine and Petroleum Geology 52 (57-76)" authored by Serrano and Guerra-Merchán

Keywords: Alboran Sea, Post-Messinian Crisis reflooding, Lago Mare, Mediterranean -Paratethys high sea-level connection

. This deposit is not a "transitional unit towards normal marine conditions" but followed the marine reflooding of the Mediterranean Basin, allowing the connection with the Dacic Basin (Eastern Paratethys). We show that this Lago Mare event did not result from a climatic change leading to a dilution episode but from a high sea-level connection. Different interpretations on the respective effects of tectonics and eustatism on the sedimentary archives linked to the Messinian Salinity Crisis (MSC) are once more debated. They depend on (1) the duration assigned to the MSC and its subdivisions, (2) the amplitude of the successive sea level variations, and at last (3) the dimensional (local or regional) view of the MSC. Such discrepancies concern subsidiary matters which cannot mask first order facts and progresses in their knowledge.

We welcome the opportunity to answer this comment that mixes up first order and secondary aspects of the Messinian Salinity Crisis which are worth distinguishing. Accordingly, our reply first concerns the new context of the Malaga Lago Mare which is partly accepted by Serrano and Guerra-Merchán. Then, we return to the old debate on the respective impact of tectonic activity and sea-level changes during and just after the Messinian Salinity Crisis (MSC) which, at last, affect the interpretation of the Messinian -Zanclean archives.

Marine context of the Malaga Lago Mare deposition

As a major conclusion of the paper by [START_REF] Couto | Lago Mare and the Messinian Salinity Crisis: Evidences from the Alboran Sea (S. Spain)[END_REF], we showed that the Lago Mare (LM) of the Malaga area has been deposited during a marine context of high sea level that is agreed by Serrano and Guerra-Merchán. It corresponds to the LM 3 event [START_REF] Clauzon | Influence of the Mediterranean sea-level changes over the Dacic Basin (Eastern Paratethys) in the Late Neogene. The Mediterranean Lago Mare facies deciphered[END_REF] that is not questioned by these authors. We interpret this event as the consequence of the fast marine reflooding of the Mediterranean Basin which resulted from the collapse of the Gibraltar sill [START_REF] Garcia-Castellanos | Catastrophic flood of the Mediterranean after the Messinian salinity crisis[END_REF] and rapidly allowed the connection with the Paratethyan Dacic Basin [START_REF] Bache | A two-step for the reflooding of the Mediterranean after the Messinian Salinity Crisis[END_REF]. This is the only way to explain the large amounts of mixed Atlantic and Paratethyan fossils (calcareous nannofossils, foraminifers, molluscs, and dinoflagellate cysts on the one hand, ostracods, molluscs and dinoflagellate cysts on the other hand; [START_REF] Guerra-Merchán | Messinian Lago-Mare deposits near the Strait of Gibraltar (Malaga Basin, S Spain)[END_REF][START_REF] Couto | Lago Mare and the Messinian Salinity Crisis: Evidences from the Alboran Sea (S. Spain)[END_REF]. Serrano and Guerra-Merchán interpret the marine character of the LM 3 episode as resulting from initiated brief and repeated oceanic overflows into a brackish Mediterranean Basin. On contrary, our palaeobiological data show that this LM 3 episode results from only one event, the massive entrance of oceanic waters into an almost desiccated Mediterranean Basin.

The Serrano and Guerra-Merchán's hypothesis also invokes peculiar climatic conditions at the time of LM 3 "with a positive hydrological budget (rainfall + runoff > evaporation), contrary to the present-day one", plus restricted water exchanges through the Gibraltar Strait that together might have led to water stratification [START_REF] Guerra-Merchán | Tectono-sedimentary evolution of the peripheral basins of the Alboran Sea in the arc of Gibraltar during the latest Messinian-Pliocene[END_REF]. However, an abundant pollen documentation shows that the southwestern Mediterranean climate was significantly warmer and drier than today prior to the MSC, during its first step (as defined by [START_REF] Clauzon | Alternate interpretation of the Messinian salinity crisis: Controversy resolved?[END_REF] and immediately after the MSC that includes the LM 3 event [START_REF] Fauquette | How much did climate force the Messinian salinity crisis? Quantified climatic conditions from pollen records in the Mediterranean region[END_REF]Jiménez-Moreno et al., 2010;[START_REF] Feddi | Histoire plio-pléistocène des écosystèmes végétaux de Méditerranée sud-occidentale : apport de l'analyse pollinique de deux sondages en mer d'Alboran[END_REF]). The regional documentation is now completed in deep sea cores (DSDP -ODP Sites 134, 976, and 978) by the record of the three successive LM events [START_REF] Clauzon | Influence of the Mediterranean sea-level changes over the Dacic Basin (Eastern Paratethys) in the Late Neogene. The Mediterranean Lago Mare facies deciphered[END_REF] which occurred in the same warm and dry climatic conditions (Popescu et al., submitted). The pollen content of our sample 6 from the Río Mendelín Lago Mare (Do [START_REF] Couto | Lago Mare and the Messinian Salinity Crisis: Evidences from the Alboran Sea (S. Spain)[END_REF] figure 3), characterized by 56% of herbs including subdesertic plants, locally confirms these climatic conditions. Increased runoff and subsequent water stratification (i.e., anoxic conditions on the sea floor) must result in abundant amorphous organic matter in the palynological slides (Suc et al., 2010), a phenomenon that we did not observe in the Alboran Sea. These data contradict the Serrano and Guerra-Merchán's climatic hypothesis. Thus, the assumption of the Lago Mare being a dilution event as concluded by Rouchy and Caruso, 2006 is not supported.

Tectonic versus eustatism

Serrano and Guerra-Merchán argue that tectonics may have played a major role in the stratigraphic succession observed at the Río Mendelín section. We agree that tectonics was active between the deposition of the Lago Mare and the subsequent yellowish clays (Do Couto et al., 2013: figure 3) as attested by the occurrence of syn-sedimentary normal faults. However the authors assumed that "tectonics could have played a more radical role at the end of the Messinian" [START_REF] Guerra-Merchán | Messinian Lago-Mare deposits near the Strait of Gibraltar (Malaga Basin, S Spain)[END_REF]. This is an old debate to which we are often confronted regarding the post-MSC sedimentation. Actually, the debate concerns the respective influence of the sea-level change, with outstanding amplitude in our opinion, and the tectonic activity. Regarding this, two main points must be considered:

-a local view of the MSC, particularly in a very active tectonic area such as the Gibraltar arc, leads to consider the effects of tectonics prevailing over those of eustatism. However, at the Mediterranean scale, the widespread occurrence of post-MSC Gilbert-type fan deltas observed on all types of margins, whether active or passive [START_REF] Bache | A two-step for the reflooding of the Mediterranean after the Messinian Salinity Crisis[END_REF] forces to incline towards a generalized response to this unique event, without neglecting the local synor post-sedimentary tectonics;

-the sea-level fall and rise opening and closing the MSC are estimated between and 1500 m [START_REF] Bache | A two-step for the reflooding of the Mediterranean after the Messinian Salinity Crisis[END_REF][START_REF] Gargani | Slope destabilization during the Messinian Salinity Crisis[END_REF] which, in a very short time (140 kyrs; [START_REF] Bache | A two-step for the reflooding of the Mediterranean after the Messinian Salinity Crisis[END_REF], is more than one order of magnitude larger than the maximum tectonic uplift estimated in the studied area for a significantly longer time (e.g. 100 m; Guerra-Merchán et al., 2014).

Post-MSC Gilbert-type fan deltas

Without repeating the arguments in favour of two coalescing Gilbert-type fan deltas at El Túnel and Río Mendelín already detailed in Do [START_REF] Couto | Lago Mare and the Messinian Salinity Crisis: Evidences from the Alboran Sea (S. Spain)[END_REF], we now focus on four items questioned by Serrano and Guerra-Merchán:

-according to the ages indicated by the calcareous nannoflora, there is no reason to envisage a polyphased erosional surface in which the Gilbert-type fan deltas are nested. This is actually the Messinian Erosional Surface as now accepted by Serrano and Guerra-Merchán; -the record of the calcareous nannofossil species Triquetrorhabdulus rugosus without Ceratolithus acutus at El Túnel suggests that these deltaic sediments are older than those at Río Mendelín where these species have been recorded together (Raffi et al., 2006), i.e. sample 5 (Do Couto et al., 2014: table 1); -the thin iron-rich horizons in the Río Mendelín section, considered as "hard-ground surfaces" in relation with tectonic movements, are very common concentrated horizons in the bottomset beds of the post-MSC Gilbert-type fan deltas, more unsurprising with respect to the sierras as a source in the hinterland; -whereas Serrano and Guerra-Merchán attribute the conglomerates topping the Río Mendelín section to "post-Zanclean probably Upper Pleistocene (Tyrrhenian)" age, there is in fact no datation contradicting our ascription to the topset beds of the Gilbert-type fan delta.

Conclusion

We mainly reiterate the conclusions of the Do Couto et al. ( 2014)'s paper, particularly concerning the full marine context of the Lago Mare deposit at Río Mendelín which immediately overlies the Messinian Erosional Surface and belongs to a Gilbert-type fan delta. This Lago Mare episode is the LM 3 defined by [START_REF] Clauzon | Influence of the Mediterranean sea-level changes over the Dacic Basin (Eastern Paratethys) in the Late Neogene. The Mediterranean Lago Mare facies deciphered[END_REF] which immediately followed the brutal marine reflooding of the Mediterranean Basin, that momentarily allowed connection between the Mediterranean Sea and Dacic Bassin ( Popescu et al., 2009;[START_REF] Bache | A two-step for the reflooding of the Mediterranean after the Messinian Salinity Crisis[END_REF]. We appreciate the agreement of Serrano and Guerra-Merchán on some of these major points and we thank them to have offered us the opportunity to reassert our interpretations such as the LM 3 is not a "transitional unit towards normal marine conditions" but was deposited in a full marine context. Some minor items depend on interpretation, such as (1) the two Lago Mare sequences for Serrano and Guerra-Merchán that we consider the prograding passage from foreset to bottomset beds, or (2) the origin of the erosional contact between the Lago Mare deposits and the Zanclean ones referred by Serrano and Guerra-Merchán to a sea-level fall and to submarine gravity instability by us.

Jiménez-Moreno, G., Fauquette, S., Suc, J.-P., 2010

  . Miocene to Pliocene vegetation reconstruction and climate estimates in the Iberian Peninsula from pollen data. Review of Palaeobotany and Palynology 162, 403-415. Popescu, S.-M., Dalesme, F., Jouannic, G., Escarguel, G., Head, M.J., Melinte-Dobrinescu, M.C., Sütő-Szentai, M., Bakrac, K., Clauzon, G., Suc, J.-P., 2009. Galeacysta etrusca complex, dinoflagellate cyst marker of Paratethyan influxes into the Mediterranean Sea before and after the peak of the Messinian Salinity Crisis. Palynology 33(2), 105-134. Popescu, S.-M., Dalibard, M., Suc, J.-P., Barhoun, N., Melinte-Dobrinescu, M.C., Bassetti, M.A., Deaconu, F., Head, M.J., Gorini, C., Do Couto, D., Rubino, J.-L., Auxietre, J.-L., Floodpage, J., submitted. Messinian -Zanclean boundary and Lago Mare events in the deep southwestern Mediterranean. Marine and Petroleum Geology. Raffi, I., Backman, J., Fornaciari, E., Pälike, H., Rio, D., Lourens, L., Hilgen, F., 2006. A review of calcareous nannofossil astrobiochronology encompassing the past 25 million years. Quaternary Science Reviews 25, 3113-3137. Rouchy, J.M., Caruso, A., 2006. The Messinian salinity crisis in the Mediterranean basin: a reassessment of the data and integrated scenario. Sedimentary Geology 188-189, 35-67. Suc, J.-P., Combourieu-Nebout, N., Seret, G., Popescu, S.-M., Klotz, S., Gautier, F., Clauzon, G., Westgate, J., Sandhu, A.S., 2010. The Crotone series: a synthesis and new data. Quaternary International 219, 121-133.

Acknowledgements

Field trips and micropaleontological analyses have been supported by the CIFRE PhD grant No 584/2010 (TOTAL/UPMC) and by the CNRS/INSU "Actions Marges" Program.