Frédéric Bernardin

Mireille Bossy

Claire Chauvin

Philippe Drobinski

Antoine Rousseau

Tamara Salameh

Stochastic Downscaling Method: Application to Wind Refinement

Keywords: stochastic differential equations, downscaling method, wind simulation

In this article, we propose a new stochastic downscaling method: provided a numerical prediction of wind at large scale, we aim to improve the approximation at small scales thanks to a local stochastic model. We first recall the framework of a Lagrangian stochastic model borrowed from S.B. Pope. Then, we adapt it to our meteorological framework, both from the theoretical and numerical viewpoints. Finally, we present some promising numerical results corresponding to the simulation of wind over the Mediterranean Sea.

Introduction

Near-surface wind speeds are acknowledged as having particular importance on society (e.g., the insurance industry, coastal erosion, forest and infrastructure damage, storm surges, and air-sea exchange). Estimating near-surface wind speeds also has relevance for applications such as pollutant diffusion evaluation, wind energy resource estimation and construction issues. Surface wind speeds however exhibit variability at much smaller spatial scales than that resolved by atmospheric global circulation models and hence there is a need to develop tools for downscaling methods to generate finer scale near-surface wind climatologies. Two major methods are used to generate fine scale wind: dynamical downscaling and statistical downscaling. The dynamical downscaling consists in resolving with a numerical model the conservation equations of the physics over a specific geographical region with spatial scales ranging from few kilometers to tens of kilometers [START_REF] Salameh | Aerosol distribution over the western mediterranean basin during a tramontane/mistral event[END_REF], Žagar et al. (2006)). Statistical downscaling of wind, empirically estimates fine scale wind speed [START_REF] Pryor | Empirical downscaling of wind speed probability distributions[END_REF]) or wind components [START_REF] Salameh | Statistical downscaling of near-surface wind over complex terrain in southern france[END_REF]) by deriving statistical relationships between observed small-scale (often station level) variables and larger scale variables from atmospheric global circulation models.

In this paper, we propose a new method for the simulation of wind at small scales. Based on an existing numerical weather prediction (NWP) model, we introduce a system of stochastic differential equations (SDEs) as a local model used to implement a downscaling method: given a coarse prediction provided by the NWP model, the stochastic downscaling method (SDM) is aimed to estimate the wind distribution at small scales.

The development of SDM is motivated by its capacity to refine mean atmospheric winds, which are solved at rather coarse resolution (at best about 10 km horizontal resolution and 50 m vertical resolution in the first 1000 m of the atmosphere) by NWP models, and estimate their associated uncertainties.

The NWP model solves the equations for the atmosphere (momentum and mass conservation, equation of state, thermodynamics law, and water vapor continuity). The evolution of the model circulation is computed by time integration of those equations starting from an initial condition and forced at its boundaries by atmospheric re-analyses which are projections of the state of the atmosphere as known from the finite set of imperfect, irregularly distributed observations onto a regular grid. The NWP model used in this study is the fifth generation Penn State-National Center for Atmospheric Research MM5 model, version 3.6 (Dudhia (1993)).

Given the coarse prediction (provided here by MM5), SDM consists in simulating numerous fluid particles in a chosen domain of computation D, in order to refine the numerical prediction over it. The associated Lagrangian model introduces some randomness in order to fit the turbulent behavior of the fluid particle with the turbulence closure of the NWP model. We borrow and adapt models first introduced by S.B. Pope [START_REF] Pope | Lagrangian PDF methods for turbulent flows[END_REF]) that have been widely used in the framework of multiphasic flows.

The article is organized as follows: in Section 2, we give a general description of Pope models, recall some recent mathematical studies related to these models, and set a numerical framework that will be used in the sequel. In Section 3, we give a description of the SDM algorithm. We adapt the Lagrangian stochastic models to our specific case: issues of turbulence closure and boundary conditions are addressed. Finally, we end with an application to wind refinement, for which we present some encouraging numerical result in Section 4.

2 Lagrangian Models and Associated Algorithm

General Description

We now give a general description of Lagrangian stochastic models for turbulent flows. For more details, the reader is referred to [START_REF] Pope | Lagrangian PDF methods for turbulent flows[END_REF][START_REF] Pope | Turbulent flows[END_REF] and the related bibliography.

Consider an open set D of R 3 . We are interested in the behavior of an incompressible fluid in D, with constant mass density ρ. The case where ρ could eventually vary is left to further studies (see e.g. [START_REF] Pope | Lagrangian PDF methods for turbulent flows[END_REF]). At any time t > 0 and position x in D, the fluid has a velocity U in R 3 , and a pressure P in R, both ruled by the Navier-Stokes equations:

∂ t U + (U • ∇) U = - 1 ρ ∇P + ν△U, (1a)
∇ • U = 0, (1b) Boundary Conditions, (1c) Initial Data. ( 1d
)
To take into account the turbulent behavior of the fluid when the Reynolds number is high (i.e. the viscosity ν is tiny), the physical fields are not only supposed to depend on position x and time t, but they are also considered as random variables. The fields can thus be written, following the so-called Reynolds decomposition, as the sum of a deterministic part and a ω-dependent part, namely:

U (t, x, ω) = U (t, x) + u(t, x, ω), (2) P(t, x, ω) = P (t, x) + p(t, x, ω), ( 3
)
where ω is a random state living in a set of possible realizations Ω, equipped with a set of events F, measured by a probability P.

Reynolds averaging is described as a linear operator applied to random fields, which is assumed to commute with spatial and times derivatives. Applying the Reynolds operator to the Navier Stokes equations (1) leads to the so-called Reynolds Averaged Navier Stokes (RANS) equations:

∂ t U (i) (t, x) + ( U (t, x) • ∇) U (i) (t, x) =ν△ U (i) (t, x) - 1 ρ ∂ x i P (t, x) - 3 k=1 ∂ x k u (i) u (k) (t, x), (4a) ∇ • U (t, x) =0, ( 4b
)
where U (i) stands for the ith component of U for i = 1, 2, 3, and P is solution of the following Poisson equation:

- 1 ρ ∆ x P = 3 i,j=1 ∂ x j U (i) ∂ x i U (j) + ∂ 2 x i x j u (i) u (j) , (5)
and where the matrix { u (i) u (j) } 1≤i,j≤3 , called the Reynolds tensor, is defined as:

u (i) u (j) = U (i) U (j) -U (i) U (j) , ∀ 1 ≤ i, j, ≤ 3.
The system (4)-( 5) is not closed. One needs to provide additional parametric models for the Reynolds tensor, based both on theory and experimental observations. As examples, we mention the turbulent-viscosity models and related Smagorinsky closures or k -ε closures, or Reynolds-stress models like Rotta's model. The interested reader is referred to [START_REF] Mohammadi | Analysis of the k-epsilon turbulence model[END_REF], [START_REF] Pope | Turbulent flows[END_REF] for further details on the closure of the RANS equations.

Starting from 1985, S.B. Pope proposes an alternative approach to the RANS equations. His idea consists in describing the fluid as an ensemble of Lagrangian particles, and computing by a particle method the averaged fields solution to the RANS equations. In the case of an incompressible fluid, with constant density ρ, we assume that there exists a probability density f E (t, x; V, θ) that governs the Eulerian dynamic. In that case, the Reynolds operator writes:

U (t, x) = R 3 ×R V f E (t, x; V, θ)dV dθ, U (i) U (j) (t, x) = R 3 ×R V (i) V (j) f E (t, x; V, θ)dV dθ.
Pope describes f E by modeling a Lagrangian fluid particle, introducing some randomness in what becomes a stochastic dynamic whose state variable is (X t , U t , φ t ). These three stochastic processes respectively correspond to the position, the velocity and, possibly, some other physical scalar characteristic of the fluid particle at time t. The fundamental idea of Pope is to construct processes whose laws allow the reconstruction of the moments of the physical quantities, up to a given order. To this aim, the dynamic of the fluid particle is ruled by a system of stochastic differential equations (SDEs) of the form:

dX t = U t dt, ( 6a
)
dU t = - 1 ρ ∇ x P (t, X t )dt + D U (t, X t , U t , φ t )dt + B U (t, X t , U t , φ t )dW t , ( 6b
)
dφ t = D φ (t, X t , U t , φ t )dt + B φ (t, X t , U t , φ t )d W t , (6c)
where (W t , t ≥ 0) is a three dimensional Brownian motion and ( W t , t ≥ 0) is a one dimensional Brownian motion independent of W (see e.g. Øksendal (1995) for details on SDEs driven by Brownian motion). The averaged pressure P (t, x) is closed by the Poisson equation ( 5).

Now we precise the generic relation between the Lagrangian and Eulerian

quantities involved by the previous Lagrangian model ( 6) and the original Eulerian model ( 4).

We assume that there exists a Lagrangian density f L (t; x, V, θ), such that, at every time t, the measure f L (t; x, V, θ)dxdV dθ is the law of the random variables (X t , U t , φ t ) solution of (6).

In the case of an incompressible fluid with constant density ρ, we must have, for any x in D,

R 3 ×R f L (t; x, V, θ)dV dθ = ρ, (7)
and the link between f E (t, x; V, θ) and f L (t; x, V, θ) is the conditional relation (see e.g. [START_REF] Pope | On the relationship between stochastic Lagrangian models of turbulence and second-moment closures[END_REF]):

f E (t, x; V, θ) = f L (t; x, V, θ) R 3 ×R f L (t; x, V, θ)dV dθ = 1 ρ f L (t; x, V, θ). (8)
In other words, f E is the conditional density of (X t , U t , φ t ), knowing that X t = x. Equivalently, for any Lagrangian fields Q depending on the velocity and the scalar characteristic

Q(U, φ) (t, x) = E Q(U t , φ t ) X t = x ,
where E denotes the mathematical expectation, related to the underlying probability space (Ω, F, P) associated to the Lagrangian model ( 6). The expected velocity is evaluated conditionally to the position of the Lagrangian particle. As a consequence, the coefficients of the Lagrangian SDEs (6), D U , B U , B φ and D φ , must be designed such that the Lagrangian law (in particular the first and second conditional momentum of the velocity) is consistent with the RANS equations submitted to a given closure. This methodology, referenced as pdf method in papers, selects a family equations ( 6) and proposes a one-to-one construction between a Lagrangian stochastic model and closed RANS equations. Here, one-to-one correspondence means almost-equivalent description of momentums up to a given closure order; the SDEs coefficients are chosen according to the closure, with the help of Kolmogorov theory and experimental observations for calibration consideration (see [START_REF] Pope | On the relationship between stochastic Lagrangian models of turbulence and second-moment closures[END_REF], [START_REF] Pope | Turbulent flows[END_REF] for more details). The choice is not unique. For instance, for the velocity equation (6b), the selected family equations is the following Mean-Reverting SDE

dU (i) t = -1 ρ ∂ x (i) P (t, X t )dt + ν∂ 2 x (j) ,x (j) U (i) (t, X t )dt +G ij U (j) t -U (j) (t, X t ) dt + C 0 ε(t, X t )dW (i) t , (9)
for a specified model constant C 0 and relaxation tensor G ij , and where ε is the turbulent kinetic energy dissipation rate. A straightforward computation shows that first order momentum U (i) of ((6a), ( 9)) solves exactly (4b). For the SDM methodology, our particular choice is discussed in section 3.1.

Mathematical studies of stochastic Lagrangian models

The study of SDEs such as (6) raises many strong mathematical difficulties, as for the RANS equations in the framework of turbulent flows.

A study of the well-posedness of a simplified Lagrangian model can be found in Bossy et al. (2008). Though simple, the model under study has a Langevin structure and the spatial dependency in the coefficients of the velocity equation is a conditional expectation with respect to the position:

X t = X 0 + t 0 U s ds U t = U 0 + t 0 E(Q(U s )/X s )ds + σW t .
(10) for a given bounded and smooth function Q and a fixed σ > 0. Equation (10) has a unique solution constructed as the limit of a smoothed system, working with a smoothed version of the conditional expectation such that

E (Q(U t )φ δ (X t -x)) E (φ δ (X t -x)) tends to E(Q(U t )/X t = x) = Q(U ) (t, x),
as δ tends to 0, where (φ δ (x), δ > 0) is any regular approximation of the Dirac mass at the origin. The smoothed system is approximated itself as the limit equation of a system of weakly interacting particles ((X

i,Np t , U i,Np t ), N p ∈ N): X i,Np t = X i,Np 0 + t 0 U i,Np s ds, (11a)
U i,Np t = U i,Np 0 + t 0 1 Np Np j=1 Q(U j,Np s )φ δ (X i,Np s -X j,Np s ) 1 Np Np j=1 φ δ (X i,Np s -X j,Np s ) ds + σW i t , (11b)
where (W i • , i ∈ N) is a family of independent 3D-Brownian motions. Equation (11) can be viewed as a kind of spatial discretization of Equation ( 10) and the particle summations in (11) are Monte Carlo approximations of the conditional expectation in (10), excepted that here, a propagation of chaos result (see Bossy et al. (2008)) replaces the law of large number that induces the Monte Carlo convergence and

1 Np Np j=1 Q(U j,Np t )φ δ (X i,Np t -X j,Np t ) 1 Np Np j=1 φ δ (X i,Np t -X j,Np t ) (12)
is a converging estimator of E(Q(U t )/X t = x), as N tends to ∞ and δ tends to 0. We refer the interested reader to [START_REF] Bossy | Some stochastic particle methods for nonlinear parabolic pdes[END_REF] for a review on convergence results and rate of convergence results for mean field particles approximations.

For our downscaling method in a given bounded domain D, we need to construct a Lagrangian model confined in D and that admits a prescribed velocity at the boundary ∂D . As an example, Bossy and Jabir (2008) show that a Lagrangian model like (10) satisfies the no-permeability condition at the boundary ∂D

E (U t • n ∂D /X t = x) = 0, for any (t, x) ∈ [0, T ] × ∂D (13)
if its dynamic is modified at the boundary as

   X t = X 0 + t 0 U s ds, U t = U 0 + t 0 β(U s )ds + σW t - 0<s≤t 2(U s -• n ∂D (X s ))n ∂D (X s ) l l {Xs∈∂D} . (14)
The vector n ∂D above denotes outward normal at the boundary. The term with jumps 0<s≤t 2(U s -• n ∂D (X s ))n ∂D (X s ) l l {Xs∈∂D} correspond to a symmetrization of the normal velocity, when the particle reaches the boundary. Showing that ( 14) satisfies ( 13) requires some a priori existence result for the trace problem associated to the probability density function of the particle. Such result is obtained in Bossy and Jabir (2008) for some smooth drift function β(u).

Numerical methodology

In this section we present the generic numerical algorithm that will be specifically implemented for SDM (see Section 3 below). The objective is thus to discretize Equations ( 6), ( 7) and ( 16). To this aim, according to the initial data, we drop N p fluid particles in the domain D and mesh D into cells (C i , i = 1, . . . , N c ). The approximated Eulerian quantities Q(U, φ) (where Q(U, φ) stands for any Lagrangian field) are computed thanks to the following average: for all x ∈ D,

Q(U, φ) (t, x) ≃ 1 N p Np k=1 Q U k,Np t , φ k,Np t G(X k,Np t , x) 1 Np Np j=1 G(X k,Np t , X j,Np t ) ,
where G is a function that quantifies the contribution of the particles to the Eulerian local information. For the Nearest Grid Point method, G(•, x) is the indicator function 1 C ix (•) where i x is such that x ∈ C ix and

Q(U, φ) (t, x) ≃ Np k=1 Q U k,Np t , φ k,Np t 1 C ix (X k,Np t ) #{X j,Np t ∈ C ix , j = 1, . . . , N p } .
Note that this smoothing function G is an efficient numerical choice for the function φ δ in the estimator (12). The above notation Y k,Np , for Y = X, U or φ, stands for the property Y of the particle k among the N p particles.

Adapted from the ideas of [START_REF] Pope | PDF methods for turbulent reactive flows[END_REF], the main steps of our numerical algorithm are the following (from time t n = n ∆t to time t n+1 = t n + ∆t).

Given the values of the processes (X (n) , U (n) , φ (n) ) at time t n , we:

1. compute new values of velocity and physical quantities φ. To this aim, we compute ( U (n+1) , φ (n+1) ) thanks to Equation (6c) and the following modified Equation (6b)

d U t = D U (t, X t , U t , φ t ) dt + B U (t, X t , U t , φ t ) dW t , (15)
which does not take the pressure gradient into account. This will be done in the correction step for the velocities (see item 4 below).

The simulation of such stochastic differential equations inside D can be done thanks to an Euler scheme (e.g. [START_REF] Talay | Probabilistic numerical methods for partial differential equations: elements of analysis[END_REF]). Since we are considering the long time behavior of a potentially stiff problem, we have to implement a so-called exponential scheme (based on a local linearization of the drift and diffusion terms). We refer the interested reader to [START_REF] Mora | Weak exponential schemes for stochastic differential equations with additive noise[END_REF] and references therein.

In the particular case of SDM, where the domain D is bounded, a particular care must be taken for the boundary conditions, namely when the particles reach the boundary of D (see Section 3.2 below).

2. move particles: compute X (n+1) according to Equations (6a).

3. modify the positions from X (n+1) to X (n+1) , in order to satisfy (7). Note that (7) implies that the law of the particle position at any time must be uniform in the bounded domain D. This step is not classical at all and concerns the particle positions. The aim is to correct the particle positions in the domain such that:

• their probability density is uniform on D, which is insured numerically by maintaining the same number N pc = N p /N c of particles in each cell);

• the "transport cost" is minimum. This is a problem of discrete optimal transportation, which is known to be nonlinear and numerically very difficult in dimension 3 (see [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF], [START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Mccann | Existence and uniqueness of monotone measurepreserving maps[END_REF]), whereas the one-dimensional corresponding problem is linear, and easy to handle from the numerical viewpoint.

In [START_REF] Chauvin | Solving the uniform density constraint in a downscaling stochastic model[END_REF], we also studied another method for the discrete optimal transportation problem, based on the Auction Algorithm introduced by [START_REF] Bertsekas | Linear network optimization: algorithms and codes[END_REF][START_REF] Bertsekas | Auction algorithms for network flow problems: A tutorial introduction[END_REF]. One forthcoming goal would be to plug the Auction Algorithm routines into SDM.

By then, the present implemented method is based on a kill/build procedure. Namely, each cell is investigated independently of the others, in order to point out the cells that contain more (resp. less) than N pc particles. In each supernumerary cell, we kill particles so that the resulting number of particles in this cell is exactly N pc . Calling N killed the total number of killed particles, we then drop N killed new particles in the other (subnumerary) cells. Finally all the cells contain exactly N pc particles. This procedure is performed so that the local Eulerian quantities are not modified.

4. modify the velocities from U (n+1) to U (n+1) in order to satisfy the divergence free constraint

∇ x • U (t, x) = 0. ( 16
)
This correction seems a standard issue in computational fluid dynamics [START_REF] Guermond | Calculation of incompressible viscous flows by an unconditionally stable projection FEM[END_REF]): the aim is to obtain a divergence free velocity field. This could be done thanks to the resolution of the following Poisson problem:

∆P = 1 ∆t ∇ • U (n+1) . ( 17
)
For the SDM methodology, additional compatibility constraint with the MM5 boundary values is added insuring well-posedness of our Poisson problem.

Note that the problem above only concerns Eulerian unknowns, and hence does not contribute significantly to the numerical cost of the algorithm.

We then compute U (n+1) thanks to

U (n+1) = U (n+1) -∆t ∇P. ( 18
)
The Eulerian velocity field is then such that div U (n+1) = 0.

Beside PDF methods have been widely used in the context of multiphasic turbulent flows, few analytical works exist from our knowledge on the convergence of the method. In [START_REF] Xu | Assessement of numerical accuracy of PDF/Monte Carlo methods for tubulent reacting flows[END_REF], numerical experiments show that the statistical error converges as

N -1/2 pc
, where we recall that N pc denotes the number of particles per cell. For the convergence rate analysis of stochastic particle methods, see [START_REF] Bossy | Some stochastic particle methods for nonlinear parabolic pdes[END_REF] and references therein.

Description of SDM

In this section, we adapt the generic Lagrangian fluid particle model ( 6) to the framework of wind simulations at small scales.

Turbulence Closure

In the work of [START_REF] Dreeben | Wall-function treatment in PDF methods for turbulent flows[END_REF], the mean turbulent frequency is determined after solving an equation on the mean dissipation (see Equations ( 11) and ( 12) page 2 in the quoted reference). This would lead, in our generic equation ( 6), to substitute φ with the turbulent frequency ω. This has been done in [START_REF] Rousseau | Stochastic particle method applied to local wind simulation[END_REF], but it is not well adapted to the physical configuration that we are facing, since the involved time scales are far too tiny. Here, we proceed differently, and only consider equations for the position and velocity fields:

dX t = U t dt, ( 19a
)
dU t = - 1 ρ ∇ x P (t, X t )dt - 1 2 + 3 4 C 0 ω (t, X t ) (U t -U (t, X t )) dt (19b) + C 0 ε(t, X t )dW t .
where ε is the turbulent kinetic energy dissipation rate. Compared to equation ( 9), the kinematic viscosity ν is set to zero, since the Reynolds number is very high.

Naturally, it remains to provide an equation for ω , which is defined in [START_REF] Dreeben | Wall-function treatment in PDF methods for turbulent flows[END_REF] as:

ω (t, x) = ε(t, x) k(t, x) . ( 20
)
We start with the definition of the turbulent kinetic energy k, which can be calculated from the first an second moments of the velocity:

k(t, x) = 1 2 3 j=1 U (j) 2 (t, x) -U (j) (t, x) 2 . ( 21
)
A direct dissipation calculation, as suggested in [START_REF] Dreeben | Wall-function treatment in PDF methods for turbulent flows[END_REF], requires that all scales experiencing dissipation must be resolved. These scales include eddies at the Kolmogorov microscale which are typically less than a millimeter. The time increment ∆t required to resolve eddies on that scale is less than 10 -3 seconds (e.g. [START_REF] Piper | Surface layer turbulence measurements during a frontal passage[END_REF]). This value for ∆t is far too small in the context of atmospheric wind refinement on typical integration periods of few hours to few days. We thus use a closure relation which is classically used in meteorology [START_REF] Cuxart | A multiscale turbulence scheme apt for LES and mesoscale modelling[END_REF]), and assumes that turbulent eddies in the inertial range are resolved. This closure links the turbulent kinetic energy to the dissipation ε:

ε(t, x) = C ε L k 3/2 (t, x), (22)
where C ε is to be chosen and L is a characteristic length scale. For freestream turbulence, using spectral arguments, it can be shown that L = (∆x∆y∆z) 1/3 (L is the mesh grid size with ∆x, ∆y, and ∆z the mesh sizes in the x, y and z directions, respectively) [START_REF] Schmidt | Coherent structure of the convective boundary layer derived from large-eddy simulation[END_REF]). [START_REF] Redelsperger | Méthode de représentation de la turbulence l'échelle inférieure à la maille pour un modèle tridimensionnel de convection nuageuse[END_REF] use C ε = 0.7 and [START_REF] Schmidt | Coherent structure of the convective boundary layer derived from large-eddy simulation[END_REF], C ε = 0.845. In terms of a commonly used closure based on turbulent- [START_REF] Schmidt | Coherent structure of the convective boundary layer derived from large-eddy simulation[END_REF]). The values of C k found in the literature vary but [START_REF] Krettenauer | Numerical simulation of turbulent convection over wavy terrain[END_REF] found a weak sensitivity of numerical results to the values of C k in the commonly used range. [START_REF] Redelsperger | Méthode de représentation de la turbulence l'échelle inférieure à la maille pour un modèle tridimensionnel de convection nuageuse[END_REF] use C k = 0.067, [START_REF] Schmidt | Coherent structure of the convective boundary layer derived from large-eddy simulation[END_REF], C k = 0.086 and [START_REF] Deardorff | Stratocumulus-capped mixed layers derived from a three-dimensional model[END_REF], C k = 0.1. These values are valid far from the ground, but in the absence of stratification, [START_REF] Redelsperger | A simple and general subgrid model suitable both for surface layer and free-stream turbulence[END_REF] and [START_REF] Drobinski | Evaluation of a planetary boundary layer subgrid-scale model that accounts for near-surface turbulence anisotropy[END_REF] show that standard closure scheme, using the same constants as for free-stream turbulence, matches the Monin-Obukhov similarity laws (i.e. the wall law) providing that the mixing length is taken as L = Az, where A = κC

viscosity model ν = C k L √ k, C ε /C k = π 2 (
1/4 ε C -3/4 k = 2.8
(with κ = 0.4 is the Von Karman constant). This length scale is thus taken much larger than the Prandtl mixing length L = κz, which is physically explained by the anisotropic nature of near-surface turbulence [START_REF] Carlotti | Length-scales in wall-bounded high reynolds number turbulence[END_REF]).

The parameter C ε , which we set to 0.7 in the following, is the only degree of freedom of our model, and L = 2.8 z.

The MM5 boundary conditions

We denote by ∂D the boundary of D. To take into account the (Eulerian) boundary conditions given by the NWP model MM5, it is necessary to add some terms in the equation (19b). To lighten the notations, we rewrite it under the form (6b) and then modify it by:

U t =U 0 + t 0 D U ds + t 0 B U dW s -2 0≤s≤t U s-l l {Xs∈∂D} + 0≤s≤t V MM5 (s, X s ) l l {Xs∈∂D} , (23)
where V MM5 (t, x) denotes the field of MM5 velocities at ∂D. The process U is no more a time continuous process, and U t-denotes the left limit at time t. These two new terms, involving sums indexed by the particle's hitting times of the boundary ∂D, lead to the following Eulerian property satisfied by the Lagrangian velocity (see Bossy and Jabir (2008)):

U (t, x) = V MM5 (t, x), x ∈ ∂D. ( 24
)
The time-discretization of ( 23) is sketched between t n and t n + ∆t as follows:

1. a prediction step is performed, computing ( X (n+1) , U (n+1) ) as in the free-case (see the points 1. and 2. of Section 2.3) by a discretization of the stochastic differential equation

d U t = D U dt + B U dW t ;
2. a correction step is needed if the predicted position X (n+1) is outside D. In such a case, we denote by x out the point of ∂D which is the intersection of the segment X (n) , X (n+1) and ∂D. The velocity U (n+1) is then updated by: is outcoming (its inner product with the outward normal vector of D at x out is nonnegative) we do not apply the update (25). Instead, we kill the particle and reinjects it in a cell belonging to C incom , where C incom is the set of boundary cells for which the MM5 velocity is incoming. The choice of such a cell is determined by a random sampling according to the following discrete probability law. Its support is C incom and the weight for each cell C k ∈ C incom is proportional to the (incoming) normal part of the MM5 velocity associated to C k .

U (n+1) new = V MM5 (t n , x out ) -U (n+1) . (25)
4 Numerical experiments: application to wind refinement in a realistic case

In this section, we compare simulations of stochastic and deterministic downscaling methods to the experimental buoy. First, we detail the deterministic downscaling with the NWP model MM5. This model is used both to provide coarse resolution boundary data to SDM, and to run the small scales simulations for comparison. Next, we give some criteria for the numerical convergence of SDM. We end this section with a discussion on the results of the performed numerical simulations.

Deterministic downscaling with MM5

The NWP model MM5 is run for three days between March 23rd and 25th 1998, over the area shown in Figure 1, with three embedded domains whose horizontal resolution is respectively 27, 9 and 3 km. The initial and boundary conditions are taken from reanalyses provided by the European Centre for Medium Range Weather Forecast (ECMWF). An instrumented buoy (ASIS) measuring wind speed and direction at 7 m above the sea-surface, is located at 4.25 • E/43.0 • N , as shown in the bottom right picture of Figure 1. Experiments done on the buoy's performance showed its capability of measuring momentum fluxes, wind profiles, and wave directions of high quality, in open seas as well as under severe oceanic and wind conditions (greater than 18 m•s -1 ), [START_REF] Graber | Asis-a new air-sea interaction spar buoy: design and performance at sea[END_REF]). A detailed description of the NWP model and the meteorological data available for the NWP model validation at this period can be found in [START_REF] Salameh | Aerosol distribution over the western mediterranean basin during a tramontane/mistral event[END_REF].

Stochastic Downscaling with SDM

We split the computational domain into N x × N y × N z cells. The initial and boundary conditions of SDM are provided by the NWP model at the 8 corners of the computational domain. The boundary data are updated every hour.

The SDM time step is 1 second and is run over 24 hours on March 24th 1998. A spin-up of a few iterations is necessary for SDM to reach an equilibrium state, this spin-up is thus not computationally expensive. In each cell, the Eulerian fields (mean wind components, variances, turbulent kinetic energy) are evaluated thanks to an averaging over N pc particles. Recall that for mean field particle approximations, we have a propagation of chaos result (see Section 2.2). Nevertheless, the convergence speed of the first and second moment estimators should behave (at least asymptotically) as a Monte Carlo method (independent sampling). Thus, we expect the 95% confidence interval for the mean wind components to be driven by the local turbulent variance σ2 = ( u 2 , v 2 , w 2 ) in each cell, where u, v and w are the random part of the velocity Eq. ( 2). The error between the numerical velocity U SDM and the reference solution U in each cell can be quantified with the help of the law of large numbers and the central limit theorem as:

P U -U SDM ∈ - 2σ N pc , 2σ N pc ≥ 95%.
In our experiments, the value of σ tends to stabilize for N pc 400, the maximum value of each component of σ being smaller than 1. In addition, let us recall that SDM requires no stability (CFL) condition, and that it can be run over any domain of any size (even a domain corresponding to one single cell of a NWP model).

Discussion

We force SDM with boundary data provided by the lowermost NWP model cell of dimensions 9 km×9 km×24 m, on March 24th, 1998. For the computation, we take N x = N y = N z = 3, such that the horizontal resolution in SDM is 3 km, i.e. the smallest MM5 horizontal resolution in Figure 1. To obtain a confident interval, we run SDM with N pc = 800.

We compare the zonal and meridional components obtained for this SDM configuration (solid line), to the components of the wind from the instrumental buoy (dashed line) in Figure 2. We also add the results from the NWP model with resolution of 9km (stars) and 3km (crosses).

We first emphasize the ability to evaluate, inside SDM, some statistics such as the standard deviation, as shown in Figure 2. SDM is thus a powerful tool that can improve the realism of the results, by modeling the wind variability [START_REF] Mass | Does increasing horizontal resolution produce more skillful forecast?[END_REF]). The variability shown in Figure 2 is physically relevant, and validates the SDM model, in particular the turbulent closure approach discussed in Section 3.1. Moreover, this variability seems to be linked to the fluid dynamics, and not to the magnitude of the wind components. However, and not surprisingly, Figure 2 shows that the comparison between the deterministic (MM5 3 km) and the stochastic downscaling (SDM) in this present version is premature and needs further developments discussed in the next section.

Perspectives and future works

This paper describes a new method based on the Lagrangian stochastic model by Pope and adapted to meteorological applications, which allows the refinement of wind field and gives access to the different statistical moments of the wind components. This preliminary development needs to be carried on in order to implement blocking by the surface or by complex orography and effect of vertical thermal stability as suggested by [START_REF] Das | A Lagrangian stochastic model for dispersion in stratified turbulence[END_REF]. The stochastic downscaling method could also be implemented for other types of applications, such as oceanographic problems, for which the involved phenomena are smoother and possibly better fitted to Lagrangian simulations. This will be the object of future work.

  With x out , we compute the hitting time and we integrate the new velocity U

Figure 1 :

 1 Figure 1: MM5 coarse domain (horizontal spatial resolution 27 km) limited by the external dashed lines on figure (a). Figures (b) and (c) respectively represent middle (9 km) and fine (3 km) domains. Diamond on Fig. (c) represents the location of the buoy ASIS.

Figure 2 :

 2 Figure 2: Wind components ( U , V ) measured with the buoy (dashed lines), modeled with MM5 with resolution 3 km (red stars) and modeled with SDM guided by MM5 with resolution 9 km (solid line)); every 30 minutes, we plot the SDM components with symmetric error bars of two standard deviations long.

Acknowledgments

This work was partially supported by ADEME.