Joël Puibasset
email: puibasset@cnrs-orleans.fr

Pseudo-critical or hysteresis temperature versus pore size for simple fluids confined in cylindrical nanopores

The adsorption / desorption isotherms measured in nanoporous materials generally present a hysteresis. The hysteresis shrinks upon increasing the temperature (for a given pore size) or decreasing the pore size (for a given temperature), until it finally disappears at the so called hysteresis (or pseudo-critical) temperature T h or hysteresis (or pseudo-critical) pore size R h , not to be confused with a true critical point. In this paper, a Monte Carlo approach allowed calculating the surface free energy of confined fluid along the adsorption / desorption isotherms for various cylindrical pore sizes and temperatures. A simple phenomenological model then allowed exploiting these results to determine the relation between T h and R h . The prediction is compared to various literature models and experimental data, showing agreement within uncertainties. On the other hand, the simulations cannot be used directly to predict T h and R h since they significantly overestimate the hysteresis width. The model predicts a nonlinear relation between the reduced hysteresis temperature and the inverse pore radius 12

INTRODUCTION

Confinement in nanopores strongly distorts the thermodynamic properties of simple and complex fluids. Generally, a hysteresis is observed in the adsorption / desorption isotherms, which disappears upon increasing the temperature or decreasing the pore size. This point is sometimes referred to as a critical point for fluids confined in tubular and / or disordered nanoporous substrates, by analogy with the true confinement induced critical point shift in infinite slit pores. 1 It is however emphasized that this point is not a true critical point. 2 The expressions "pseudo-critical" or simply "hysteresis disappearance point" are preferable. In the following, the hysteresis disappearance temperature and pore radius are denoted T h and R h Density Functional Theory calculations performed in slit and cylindrical pores in the slab approximation allowed Evans and coworkers 3,4 to find a linear asymptotic relation

h c h c R a T T T σ = - (1)
where T c is the bulk critical temperature, σ the molecule size, R h the pore radius, and a is around unity. This behavior is in agreement with experimental observations in MCM-41, a mesoporous silica having straight, non-interconnected pores. 5,6 Based on numerous experimental data including various fluids, Morishige and coworkers 6 proposed a phenomenological relationship between T h and R h given by Eq. 1 with a = 2. More recently, Sonwane and Bhatia 7 addressed the issue of hysteresis disappearance in cylindrical pores using the Saam and Cole 8 approach. They include the effect of curvature on the adsorbate chemical potential and surface tension, and show that the relation between T c -T h and 1/R h is not simply linear. Note that this discussion on T h refers to simple ordered cylindrical pores.

The effect of pore morphology and porous network topology is discussed elsewhere. [START_REF] Pellenq | 7th Int. Symp. on the Characterization of Porous Solids[END_REF]10 Molecular simulation is a powerful tool to study such confinement situations due to its ability to take into account the fluid / fluid and fluid / substrate interactions at a molecular level.

These interactions are believed to be the main features which determine the confined fluid thermodynamic properties. Refs [11-20] give some recent examples. However, considering the hysteresis phenomenon, metastable states are involved, as well as the transition probabilities between them. In this case, the porous substrate structure (surface roughness and pore morphology and topology) become important in determining the energetic barriers and hence the probability transitions. Unfortunately, it has long been recognized that the Molecular Simulations give poor predictions on the limit of stability of the metastable states, and generally considerably overestimate the hysteresis width.

In this paper, it is proposed to take advantage of the ability of Monte Carlo simulations in predicting the thermodynamic properties of fluid (like the free energy) in the stable as well as the metastable states. This allows the calculation of surface free energies for various states obtained at different temperatures and for several pore sizes in a simple phenomenological approach in terms of bulk and surface contributions. It is emphasized that the systematic study of interfacial phenomena is not simple, especially if fluid / substrate interfaces are involved, due to the large number of parameters. References [21-30] give some examples of such studies using molecular simulation. As proposed in a recent paper, 31 a surface quantity may be introduced which greatly simplifies the characterization of the fluid / substrate interactions since this quantity is constant along the adsorption or desorption branches of the isotherms, while the usual excess free energy varies. 31,32 The relation between both quantities being known, the fluid properties can be calculated. The dependence of this new quantity on pore size and temperature is shown to be rather simple, which allows simple interpolation. The phenomenological analysis then allows the prediction of the relation between T h and R h . The results are then compared to literature models and experimental data, and discussed.

SIMULATION RESULTS

PORE MODELS

The aim of this paper is to give a comprehensive phenomenological interpretation, in terms of bulk pressure and fluid / substrate surface tension, of the relation between T h and R h corresponding to hysteresis disappearance in adsorption / desorption isotherms. The model is then chosen to be simple. The pores are cylindrical capillaries or slit pores of radius R or width H ranging from mesoscopic (few nanometers) down to microscopic scale (<1nm). All intermolecular interactions are supposed to be Lennard-Jones like. The interaction parameters entering the potential are chosen to be those of argon for the fluid, and argon / solid carbon dioxide to model the fluid / substrate interactions. This model is frequently used in literature for fundamental issues concerning fluids at an interface. [33][34][35][36][37][38] . For more details see also Ref radius. The pore width is taken large enough (equal to 8σ Ar-Ar ) so that the two plane interfaces can be considered as independent. 31 The calculations are performed at three temperatures for the three pores (T* = kT/ε Ar-Ar = 0.77, 0.90, and 1.20), providing nine adsorption / desorption isotherms.

SIMULATION TECHNIQUES

The Monte Carlo simulation in the Grand Canonical ensemble (GCMC) [START_REF] Nicholson | Computer Simulation and the Statistical Mechanics of Adsorption[END_REF] is in principle able to mimic a real adsorption / desorption experiment. The grand potential of the fluid is fixed by the gas reservoir, and the number of adsorbed molecules fluctuates. The Monte Carlo procedure (Metropolis algorithm) generates a Markov chain of molecular configurations in agreement with the grand canonical statistical ensemble, i.e. those which minimize the grand potential. As a consequence, the most relevant configurations necessary to calculate the thermodynamic properties of a given fluid state (stable or metastable, i.e. corresponding to a local minimum) are generated. This is why GCMC give good results for the thermodynamic properties (like the local fluid density, the free energy, etc.) of a given fluid state (stable or metastable). However, if one is interested in hysteresis width, the transition probabilities between the various states have to be evaluated. The GCMC procedure, which does not explore the unfavorable configurations like the critical nucleus state, is not suited for such studies. However, a phenomenological approach allows determining trends in hysteresis width as a function of the temperature and pore size, to be explained in details later. The requirement for this approach is to determine the confined fluid grand potential as a function of its bulk properties and fluid / substrate excess free energy or surface tension. This quantity is in principle accessible by general integration procedure. 35,[START_REF] Frenkel | Understanding Molecular Simulation[END_REF] However, in order to avoid lengthy calculation of complete sets of adsorption / desorption isotherms, direct algorithms have been used which exploit a natural symmetry of the system. In the case of a slit pore, the system presents translational invariance. For tubular pores presenting physico-chemical heterogeneities, the translational invariance is lost. However, in the case where the physicochemical heterogeneities are along the pore axis, the axial symmetry may be preserved. In such situations, the usual isobaric-isothermal algorithm [START_REF] Nicholson | Computer Simulation and the Statistical Mechanics of Adsorption[END_REF][START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF] or its generalization 40,[START_REF] Puibasset | [END_REF]45 may be used, where the grand potential density (thermodynamic pressure) is imposed and the system is allowed to equilibrate in the corresponding statistical ensemble. These algorithms are shortly presented now.

For a fluid characterized by its (imposed) volume V, temperature T, and chemical potential µ (grand canonical ensemble), the differential of the corresponding grand canonical potential ( )

N TS E V T µ µ - - = Ω , ,
where E is the internal energy, S the entropy, and N the number of molecules, is given by

µ Nd SdT dV d - - Π - = Ω (2)
where the quantity dV Π is the mechanical work associated to a volume variation dV. If the system presents translational invariance along at least one direction (for instance parallel to the substrate surface for a fluid confined in a slit pore) then its volume is proportional to the extension of the system in this direction. For pores presenting an axial symmetry, the same remark remains valid. 40,[START_REF] Puibasset | [END_REF]45 As a consequence of this invariance, partial derivatives of extensive thermodynamic quantities with respect to the volume V of the system are independent of its size (homogeneous function of degree zero). The extensive quantity

( ) µ , ,V T Ω
is then a homogeneous function of the first degree of the volume (the only extensive parameter), and Euler's theorem gives:

V Π - = Ω (3)
The grand potential in the pores is then calculable by Monte Carlo simulations in the usual (for translation invariance) [START_REF] Nicholson | Computer Simulation and the Statistical Mechanics of Adsorption[END_REF][START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF] and generalized (for axisymmetric pores) 40,[START_REF] Puibasset | [END_REF]45 isobaricisothermal ensemble (constant N, Π, and T). The corresponding Monte-Carlo algorithms are made of translation trials in order to thermalize the system, and volume variations trials in order to equilibrate the system thermodynamic pressure with the external reservoir. For further details concerning the algorithm implementation, see Ref [40]. It is emphasized that the thermodynamic pressure Π (scalar quantity) defined by Eq. 2 differs from the mechanical pressure, which is an anisotropic tensor for confined fluids. However, for a homogeneous fluid, Π reduces to the bulk scalar mechanical pressure denoted p bulk , and, for a fluid at an interface, the thermodynamic pressure also contains the fluid / substrate contributions, or surface tension contributions in a classical capillarity theory point of view (see later).

The simulations were performed in simulation boxes large enough to avoid finite system size effects. For cylindrical pores, the reduced length is equal to 12. For slit pores, the reduced size is 12 in each direction parallel to the pore walls. For each run, a minimum of 5. After averaging, one can introduce the global reduced density Γ of the confined fluid, defined as the total number of particles N divided by the average reduced volume V. By convention, this volume is defined as the region of the pore where the potential is not infinite, which is equivalent to taking the pore radius or width equal to R* or H* as previously introduced. It will be emphasized later that another more physical convention will emerge, taking into account the molecular size.

ISOTHERMS AND FLUID SUBSTRATE EFFECTIVE SURFACE TENSION

The nine adsorption / desorption isotherms corresponding to the three pore sizes (cylindrical pores of radius R* = 2.5 and 4.0, and the slit pore) and the three temperatures (T* = 0.77, 0.90, and 1.20) are given in Fig. 1 as a function of the reduced thermodynamic pressure (or minus grand potential density). The scale (linear for negative pressures and logarithmic for positive pressures) has been chosen to enhance the hysteresis loop, especially in the low thermodynamic pressure initial adsorption region. As can be seen, for the highest temperatures the isotherms are all reversible. Note the presence of a vertical rise at the capillary condensation point in the slit pore. For the intermediate temperature, a hysteresis appears in the slit pore as well as in the largest pore of size R* = 4.0. The adsorption in the smallest pore (R* = 2.5) remains reversible (supercritical pore). For the lowest temperature, all isotherms present hysteresis.

The previous calculations are qualitatively in agreement with experimental observations of a pseudo-critical temperature (disappearance of the hysteresis loop) which increases with the pore size. See for instance the work of Morishige and coworkers. 5,6 However, no quantitative agreement is obtained. This is actually not surprising, since it is a well-known fact that Monte-Carlo algorithms overestimate the hysteresis loop width due to the lack of ability to sample the whole phase space (the simulation run remains trapped in local minima and is unable to jump over nucleation barriers except if they are very small). However, far from the transition regions, the problem of the phase space sampling is minimized since only one global minimum or a deep local minimum is expected. If the fluid is not too dense, the Monte-Carlo algorithm is able to sample this minimum, and good agreement is generally observed between simulation and experiments. In the following, we focus on these portions of the adsorption / desorption branches far from the transition regions (steep adsorption and desorption regions of the isotherm), in order to characterize the adsorbed fluid properties.

Let us now consider a phenomenological description of the isotherms. In the classical capillarity theory, the changes in the free energy induced by the presence of an interface are entirely contained in the definition of the effective surface tension associated to this interface (thin wall approximation). In this context, the thermodynamic pressure of a fluid confined in a capillary reads:

V A p V / / bulk γ - = Ω - = Π (4)
where bulk p is the pressure of the bulk fluid, V is the volume occupied by the fluid, γ is the fluid / substrate surface tension and A is the fluid / substrate interfacial area. The grand potential and the pressure are calculated for the same chemical potential µ. As a consequence, the interfacial free energy γ is expected to be dependent on µ, and is therefore not a constant along the isotherm. In a recent paper, 31 another interfacial quantity ζ was introduced such that:

( ) ( ) Γ Π - Γ = bulk p V A ζ (5)
with Γ the confined fluid density. The difference with the definition of γ is that the bulk pressure is taken for the same density as the confined fluid. It can easily be shown that this quantity is always equal to zero for any fluid in a uniform external potential, or for an ideal fluid confined in any given external potential. 31 It is however a priori not equal to zero for an interacting fluid in a nonuniform external potential (especially close to the substrate wall).

The interest for introducing this quantity is that it was empirically shown to be almost constant along the adsorption or desorption branches (i.e. does not depend on µ ), which is of practical interest. However, the values obtained along the gaslike and the liquidlike branches procedure on this free parameter was then applied to minimize the variations of the right hand side of Eq. 5. The quantity ζ is then determined from this constant value, taking into account that the accessible surface area A is also affected by the redefinition of the pore volume, in a way given by the pore geometry.

The optimization procedure has been applied for the nine gaslike and nine liquidlike portions of isotherms previously calculated in the isobaric-isothermal ensemble for various pore sizes and temperatures (see Fig. 1), and the Lennard-Jones bulk properties are taken from Ref [46].

Let us first focus on the new pore sizes given by the optimization. The first point to be noted is that the new values obtained for the gaslike and the liquidlike parts of a given isotherm are identical within error bars (few percent). The second point is that the values are rather independent of the temperature. This is illustrated in Table 1, where the optimized pore sizes are given for the cylindrical and slit pores at various temperatures. The last point concerns a possible interpretation of the results. In all cases, the optimized pore size corresponds to the radius where the external potential reaches 10 3 in reduced units (2.04 and 3.54 for cylindrical pores). This corresponds to the pore region accessible due to thermal agitation. Let us now focus on the results obtained for ζ. The gaslike branches give values of ζ compatible with zero within error bars (of the order of 0.02). This is expected since, for a diluted adsorbed fluid, the thermodynamic pressure (minus the grand potential density) is given by the adsorbed fluid density times kT, which is the same as the bulk pressure. For the liquidlike branches, the values are given in Table 2. The Fig. 2 shows the quantity ( )

V A / ζ ρ + Π
as a function of ρ for the nine liquidlike parts of the isotherms (symbols) and the corresponding bulk fluid equation of state for the three temperatures (lines). As can be seen, an excellent agreement is obtained over the whole liquidlike branches of the isotherms, which validates the analysis. However, a disagreement is observed for the low density part of the branches corresponding to the spinodal limit of the bulk fluid. As can be seen, the fluid density is consistently overestimated by the simulation results. One reason for this disagreement is that the Monte-Carlo algorithm may be unable to sample correctly the phase space, especially close to the spinodal, i.e. the point where the liquidlike state is expected to become unstable.

Desorption is expected to occur by bubble nucleation in the system, which is favored by local fluctuations in fluid density. However, the Isobaric-Isothermal Monte-Carlo algorithm performs volume change trials by homogeneous dilatation or contraction of the fluid, which greatly disfavors such heterogeneities. This probably explains why the Monte-Carlo algorithm remains trapped in the liquidlike state slightly beyond the expected spinodal. On the other hand, a good agreement is observed far from the spinodal where such sampling difficulties are not expected.

DISCUSSION

INFLUENCE OF TEMPERATURE AND PORE DIAMETER

The interfacial quantity ζ is given in Fig. 3 as a function of the inverse of the optimized pore radius (symbols). The slit pore model with H* = 8 was shown 31 to be wide enough to be considered as two independent plane walls, and hence as a limiting case of cylindrical pore of infinite radius. Figure 3 shows that ζ varies significantly with the pore size. The smaller the pore size, the smaller the magnitude of the effective surface tension ζ . It was shown elsewhere 31 that this dependence is linear with the inverse optimized pore radius for pores larger than 4σ Ar-Ar (in diameter), in agreement with the Gibbs-Tolman-Koenig-Buff equation [47][48][49] giving the curvature dependence of surface tension. The application of this relation to the effective surface tension ζ gives:

( ) r r δ ζ ζ - = ∞ 1 ( 6
)
where ∞ ζ denotes the effective surface tension for the plane fluid / substrate interface, and δ is the displacement of the surface of zero mass density relative to the surface of tension (Tolman's length). The values obtained for δ are around σ Ar-Ar , with a small dependence on temperature. This dependence may be modeled in its simplest form by a linear function:

* * * * T b a δ δ δ + = . ( 7
)
The quantity ∞ ζ also depends on temperature. This is shown in Fig. 4 which gives the values of ζ as a function of the temperature for the three pores ( ∞ ζ is given by the slit pore, circles).

As can be seen, a linear dependence is compatible with the data:

* * * * T b a ζ ζ ζ + = ∞ . ( 8
)
The best fit of Eqs 6-8 with data is obtained for = 4. As can be seen, the agreement is quite good except for the cylindrical pore of nominal radius 4.0 at T* = 0.90.

HYSTERESIS TEMPERATURE VERSUS PORE SIZE

The hysteresis temperature is usually defined as the temperature of hysteresis disappearance for adsorption/desorption isotherms in the grand canonical ensemble (imposed chemical potential). It is however expected that the hysteresis observed in the isobaric-isothermal ensemble disappears for the same temperature. As a consequence, the knowledge of the surface term ζ, which allows to predict the adsorption / desorption isotherm as a function of the thermodynamic pressure, gives enough information about the existence or not of the hysteresis, and hence its disappearance temperature.

The adsorbed fluid isotherm Π(ρ) is deduced from the bulk one by a translation -ζ A/V. The surface term ζ is close to zero for the gaslike branch (in other words the gaslike branch is unaffected), and negative for the liquidlike branch. As a consequence, the hysteresis loop for the confined fluid is narrower than the bulk one. This is illustrated on Fig. 5 where the confined fluid isotherms are given for three pore sizes (optimized reduced radius r* = 8.0, 4.3 and 3.0) and T* = 0.80. The bulk isotherm is also given. For this choice of parameters, the translation term -ζ A/V equals to +0.50, +0.81 and +1.0 respectively (in reduced units) for the liquidlike branches, and is negligible for the gaslike branches. The vertical lines are guides to the eyes (hysteresis). As can be seen for r* = 8.0 the confinement reduces the hysteresis of approximately ζ A/V. For r* = 4.3, this quantity exactly cancels the bulk hysteresis, which results in reversible isotherm (pseudo-critical pore). For smaller pore size, the quantity ζ A/V is larger than the bulk hysteresis, which results in an unphysical gap between gaslike and liquidlike branches. In this case, the isotherm is actually reversible, but our approach is unable to describe correctly the region of intermediate density, where the surface term ζ is not expected to remain constant, as previously discussed.

Figure 4 shows that ζ is strictly negative for the bulk critical temperature (T*=1.36). 46 The hysteresis is then expected to disappear for any pore size for a temperature lower than the bulk critical temperature, in qualitative agreement with experimental data. Note that the surface term depends on the pore size due to both the surface to volume ratio and the curvature effect. As a general trend, the larger the pore size, the smaller the surface contribution, as illustrated in Fig. 5. As a consequence, for a large pore, the hysteresis loop is close to the bulk one, and the pseudo-critical temperature at which the hysteresis disappears is close to the bulk critical temperature. Conversely, for small pores, the surface term is enhanced, and significant reduction of the pseudo-critical temperature is expected.

A quantitative analysis of the pore size dependence of the pseudo-critical temperature is now possible. The previous pore size and temperature dependence of the effective surface tension ζ is used in addition to the well-known Lennard-Jones bulk properties 46 to predict the adsorption / desorption isotherms Π(ρ) of the confined fluid for a continuum of pore sizes and temperatures, within the ranges defined by our previous analysis: pore radius larger than 2σ Ar- Ar and reduced temperature between 0.7 and 1.3. For a given pore size, h T is the temperature at which the confined fluid hysteresis disappears, i.e. the temperature for which the surface contribution ζ A/V exactly cancels the bulk hysteresis. The problem is solved numerically.

The results are given in Fig 6 (grayed area, taking into account the uncertainties). The general trend is that the larger the pore size, the larger T h . However, the relation is not simple (linear for instance), and is a priori expected to be dependent on the fluid/substrate interactions.

Experimental data could be found in the literature, mainly the work of Morishige and is mentioned that the exact value of the parameter a entering Eq. 1 is probably slightly dependent on the exact nature of the system (fluid-fluid and fluid-substrate interactions). On the other hand, the predictions proposed by Sonwane et al 7 and this work significantly deviate from a linear behavior, which suggests that Eq. 1 is probably an approximation. An interesting point is however emphasized: both models include the curvature effect on surface tension, which is believed to be relevant for the determination of the relation between the hysteresis temperature shift and the inverse pore radius. In effect, the hysteresis disappearance is closely related to the nucleation barriers, which are themselves dependent on fluid-fluid and fluidsubstrate interactions. The work by Sonwane et al 7 and this work then suggest that a concavity should be expected in the curves. However, the general trends differ in both models (opposite concavity). Further experimental and theoretical investigations are then necessary to elucidate whether the concavity is expected to be always of the same sign or if it is system-dependent, and to discriminate between the various models.

CONCLUSION

A simple molecular model is used to study Lennard-Jones fluid adsorption in cylindrical nanopores. The effects of pore size and temperature are examined. The Monte Carlo results

show large hysteresis in adsorption / desorption isotherms at low temperature and for the largest pores. On the other hand, reversible isotherms are observed for small pores and for the highest temperatures. These results are in qualitative agreement with experimental observations. However, quantitative agreement cannot be reached because the Monte Carlo simulations significantly overestimate the hysteresis width (the simulations are trapped in local minima corresponding to metastable states, and are unable to jump over the energetic barriers). The simulation data are however expected to be reliable far from the transitions where the local minima are deep. A simple phenomenological model was then used to fit the data, in terms of bulk plus surface contributions in the thermodynamic pressure versus density curves. This model then allows the determination of the temperature of hysteresis disappearance for any pore size within the range of values used in the simulations. As expected, the hysteresis temperature is lower than the bulk critical temperature, and in agreement with experimental data within uncertainties. The prediction of the pseudo-critical temperature shift as a function of the inverse reduced pore size is shown to deviate from the simple linear relations proposed by Evans et al. 3,4 (Density Functional Theory calculations)

and Morishige et al. 5,6 (phenomenological approach). This non-linearity was previously observed by Sonwane and Bhatia 7 , who predicted the opposite concavity. It is emphasized that both models include the pore size dependence of the surface contributions to the free energy (surface tension in the work of Sonwane et al. 7 and effective surface term in this work). This is believed to be an important feature to explain the relation between the hysteresis temperature shift and the inverse reduced pore size. The experimental and theoretical predictions are however insufficient to validate one approach, and one may wonder whether the predicted concavity is confirmed experimentally, and if its sign is universal or system-dependent. Further experimental and theoretical investigations are then necessary to address these issues. Lines and vertical branches are guides to the eye. The bulk equation of state (Ref [46]) is also given (line without symbols).

[

  39]. All thermodynamic quantities are normalized to the Ar-Ar parameters (reduced quantities, denoted by an asterisk). The interaction sites of the solid matrix are supposed to be uniformly distributed (smooth wall approximation), and the fluid / adsorbent interactions give rise to an external potential, denoted cylindrical pore of reduced radius R*. The obtained external potential fields for various pores sizes are given elsewhere.31,40 In this work, two pore sizes (R* = 2.5 and 4.0) are considered. The slit pore is used to model plane surfaces, which are considered to be the limiting case of cylindrical pores of infinite

10 6

 6 MonteCarlo trials were performed to acquire statistical data. The thermodynamic pressure (imposed parameter) is varied between -0.5 and 10. Negative values are possible in the liquid state for low temperatures due to the existence of cohesive forces. The highest values of the thermodynamic pressure are reached at high density due to the repulsive part of the intermolecular potential. During these simulations, the fluctuating variable is the volume.

  of the same isotherm are different. In other words, one can define a liquidlike and a gaslike value for ζ. The case of the reversible isotherms is more complicated: ζ is almost constant in the low density and the high density regions of the isotherm, but varies in the intermediate region. However, as for the irreversible isotherms, one can define gaslike and liquidlike regions corresponding to the low (<0.05) and high (> 0.3) density parts of the isotherm where ζ is almost constant. Let us now discuss another point: the density of the confined fluid entering Eq. 5 depends on the exact definition of the accessible pore volume, which is not uniquely defined for nanoscopic pores. The pore volume (actually the pore radius in cylinders and pore width in slits) was then left as a free parameter, and the new corresponding confined fluid density (number of molecules divided by the new volume) is denoted ρ. An optimization

  . It is given by the solid lines inFigs 3 and

coworkers 5 ,

 5 6 who focused on hysteresis temperature determination for simple fluids in almost cylindrical pores. The data for Ar and CO 2 adsorbed in MCM-41 are reported on Fig 6. Other gases adsorbed in MCM-41 where also considered by experimentalists (nitrogen, oxygen, etc), but not reported on the figure. They were shown to follow essentially the same behavior as argon, given by a phenomenological relationship proposed by Morishige et al 6 corresponding to Eq. 1 with a = 2, (shown in Fig 6). As can be seen in Fig 6, the CO 2 data significantly deviate from this relation, possibly due to the choice of the molecular size. The data obtained for fluid adsorbed in conventional non-cylindrical porous materials are not reported. To our knowledge, little theoretical work has attempted to elucidate this behavior. Essentially two approaches are known by the author. Density Functional Theory was used by Evans and coworkers 3,4 to deduce Eq 1 with a of the order of 1. However, the model does not give more information on the exact value of a. Note that this work actually precedes the experimental work by Morishige and coworkers. Recently, Sonwane and Bhatia 7 addressed the problem of hysteresis disappearance in cylindrical pores using the Saam and Cole 8 approach and including the effect of curvature on the adsorbate chemical potential and surface tension. Their prediction, reproduced in Fig 6, is within the range of the experimental data. Examination of Fig 6 shows that the relation initially proposed by Evans and coworkers predicts a linear behavior between the reduced critical shift and the inverse pore radius, in qualitative agreement with experimental data. However, the proportionality coefficient has to be taken close to 2 for quantitative agreement (empirical model of Morishige and Sikimi 6 ). It

FIG. 1 .

 1 FIG.1. Adsoprtion / desorption isotherms at T* = 0.77, 0.90 and 1.20 in the isobaric-

FIG. 2 .FIG. 3 .FIG. 4 .FIG. 5 .FIG. 6 .TABLE 1 :TABLE 2 :

 2345612 FIG.2. Reduced confined fluid density after optimization (ρ * ) as a function of Π*+ζ *A*/V*

0.01 0.1 1 10 -0.5 0.0 reduced density Γ Γ

ACKNOWLEDGEMENTS

My colleagues A. Delville, E. Kierlik, R. Pellenq, P. Porion and G. Tarjus are gratefully acknowledged for useful discussions. The simulations were performed locally on workstations purchased thanks to grants from Région Centre (France). The Centre de Ressources Informatique de Haute Normandie (CRIHAN, St Etienne du Rouvray, France) and the Institut de Développement des Ressources en Informatique Scientifique (IDRIS-CNRS, Orsay, France) are gratefully acknowledged.


Benjamin Dufée
email: benjamin.dufee@inria.fr

Etienne Mémin

Dan Crisan

Stochastic parametrization: an alternative to inflation in Ensemble Kalman filters

We investigate the application of a stochastic dynamical model in ensemble Kalman filter methods. Ensemble Kalman filters are very popular in data assimilation because of their ability to handle the filtering of high-dimensional systems with reasonably small ensembles (especially when they are accompanied with so called localization techniques). The stochastic framework presented here relies on Location Uncertainty (LU) principles which model the effects of the model errors on the large-scale flow components. The experiments carried out on the Surface Quasi Geostrophic (SQG) model with the localized square root filter demonstrate two significant improvements compared to the deterministic framework. Firstly, as the uncertainty is a priori built into the model through the stochastic parametrization, there is no need for ad-hoc variance inflation or perturbation of the initial condition. Secondly, it yields better MSE results than the deterministic ones.

Introduction

Data assimilation uses observational data to correct the predictions of a numerical model, to get closer to the "true" state of the system. It has been applied for a long time in weather forecasting and other geophysical dynamical systems. Numerous approches from different methodological settings have been proposed to that end. They rely on optimal control methodology, functional approximation theory, stochastic processes and filtering, or geostatistics to cite only a few of them (see [START_REF] Carrassi | Data assimilation in the geosciences: An overview of methods, issues, and perspectives[END_REF] for an overview of the techniques and recent trends in data assimilation). Among the different techniques proposed, a very popular data assimilation method is the ensemble Kalman filter (EnKF), initially devised by [START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics[END_REF], and since then explored in many different variations (see the books [START_REF] Asch | Data Assimilation: Methods, Algorithms, and Applications[END_REF][START_REF] Evensen | Data assimilation: The ensemble Kalman filter[END_REF][START_REF] Reich | Probabilistic Forecasting and Bayesian Data Assimilation[END_REF] and references therein for a precise description of all these different variations).

Ensemble methods, formulated either in variational or Bayesian forms are based on constructing an ensemble of possible realizations of a given dynamical system, with no need to explicitly compute its tangent linear model as it is required in optimal control methods [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF]. This ability is very convenient in high-dimensional and complex models, as it enables in some way to free the data assimilation technique from the very likely evolution of the code with time, in terms of update, additional modules or features, requiring a new computation of the tangent linear model at each update. On the contrary, adjoint based optimization is fully attached to a given version of the numerical code of the model. For a routine use in meteorological centers, the separation between data assimilation and the numerical codes of the dynamics constitutes an undeniable advantage. As such these methods can cope, through an extension of the linear superposition principle, with non-linear models, even if in itself, the Kalman equations are theoretically set in linear and Gaussian hypothesis. Among the different versions of the proposed ensemble Kalman filter, the so called square-root filters (SRF) consist in formulating the square root of the forecast and analysis covariance matrices in terms of the ensemble members' anomaly matrix, obtained from the residual of each member with the empirical mean [START_REF] Anderson | An ensemble adjustment Kalman filter for data assimilation[END_REF][START_REF] Bishop | Adaptive sampling with the ensemble transform Kalman filter. part I: Theoretical aspects[END_REF][START_REF] Pham | Stochastic methods for sequential data assimilation in strongly nonlinear systems[END_REF][START_REF] Tippett | Ensemble square root filters[END_REF][START_REF] Whitaker | Ensemble data assimilation without perturbed observations[END_REF]. The non uniqueness of the square-root covariance has lead to several variations of the method. The ensemble transform Kalman filter (ETKF) proposed in [START_REF] Bishop | Adaptive sampling with the ensemble transform Kalman filter. part I: Theoretical aspects[END_REF] is a popular and efficient technique based on this principle. It aims at describing the posterior ensemble members as linear combinations of the prior ones in an optimal way. Indeed, square root schemes written in terms of the anomaly matrix are particular cases of ETKF.

Applying these methods to very high-dimensional geophysical systems (see [START_REF] Houtekamer | Data assimilation using an ensemble kalman filter technique[END_REF][START_REF] Houtekamer | A sequential ensemble Kalman filter for atmospheric data assimilation[END_REF][START_REF] Houtekamer | Ensemble kalman filtering[END_REF][START_REF] Houtekamer | Atmospheric data assimilation with an ensemble kalman filter: Results with real observations[END_REF] in routine weather forecast applications) requires modest size ensembles (from tens to 100-200 at best). This is a true limitation of such ensemble Kalman filters as they might lead to strong sampling errors of the empirical covariance matrices involved in the filtering. These covariances have indeed a rank of at most p-1, with p denoting the number of ensemble members. In order to tackle this limitation and to extend their rank, so called localization approaches have been introduced [START_REF] Anderson | An adaptive covariance inflation error correction algorithm for ensemble filters[END_REF][START_REF] Anderson | Localization and sampling error correction in ensemble kalman filter data assimilation[END_REF][START_REF] Hamill | Distance-dependent filtering of background error covariance estimates in an ensemble kalman filter[END_REF][START_REF] Sakov | Relation between two common localisation methods for the enkf[END_REF]. These techniques filter out spurious correlations between two far away points and are implemented either through a Schur product of the empirical covariance matrices with a finite support positive definite kernel or by localizing explicitly the area of influence of a set of observation on a given site of the spatial domain. Both methods, although formally very different, have been shown to lead in practice to equivalent results [START_REF] Sakov | Relation between two common localisation methods for the enkf[END_REF]. The latter techniques are computationally efficient as they lead to local Kalman filters, with the risk of producing a global solution that is not a solution of the partial differential equations driving the system. The Schur product conditioning of the empirical covariance is more secured from a formal point of view but cannot be used in conjunction with squared-root filtering techniques [START_REF] Sakov | Relation between two common localisation methods for the enkf[END_REF].

Ensemble Kalman filters are, by nature, stochastic filtering techniques. They rely in one way or another on randomness as a third ingredient. The randomness considered is, in the vast majority of cases, related to the initial condition of the deterministic dynamics, with possibly another additive random variable representative of the error carried by the model itself. Such perturbations correspond to the errors attached to the physical or numerical approximation performed in the constitution steps of the system's dynamics. One of the main difficulty faced by the random dynamics built from perturbations of the initial condition is that they usually exhibit an insufficient spread of the ensemble members. This flaw leads to serious problems for data assimilation as all the ensemble members may quickly drift far away from the current observation. In high dimensional spaces, it cannot correct sufficiently the ensemble toward the current observation and therefore can yield to filter divergence. The introduction of additive random variables or inflating the ensemble variance by some factor is called variance inflation and is commonly used to tackle deficient variance representation. The problem associated with such a correction is that variance is gradually added to the system across time without being explicitly dissipated. This injection of energy leads eventually to a blowup of the dynamics simulation if it is not controlled. As a result, only very low (relatively to a given dynamics) factors of inflation can be considered in practice. Another problem, more specifically related to additive inflation, is that it corresponds to a model error variable that is blindly added to the system through an increase of variance. No particular spatial structure for the error is indeed injected. The last potential pitfall associated with such a technique is that the injection of an additive random forcing may completely change the attractor or the probability distribution of the underlying system's dynamics even for a noise of small amplitude [START_REF] Chapron | Large-scale flows under location uncertainty: a consistent stochastic framework[END_REF], leading consequently to a quite different system from the one of interest. Still this is different from the multiplicative inflation mentioned before, which is the most commonly used in practice.

In this work, we explore an alternative to the randomization of the dynamics described above. Compared to the techniques mentioned previously, the error term in the dynamics is not added as an additional forcing but is, instead, part of the model. In a similar way to the Reynolds decomposition principle used in large eddies simulation setting, the stochastic formalism on which we rely, referred to as modeling under location uncertainty (LU) [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] in the following, consists in decomposing the dynamics in terms of a smooth in time component and a highly oscillating random component that is assumed uncorrelated at the characteristic time scale of the smooth component. A large-scale dynamics which takes into account the random components is then inferred from this decomposition and usual conservation laws shall hold for the stochastic model. This setting offers a stochastic parametrization, representing the impact of the small-scale flow on the global dynamics by a random small-scale velocity. This way to represent the small scales independently, and not only subordinate to the large scales, can be put in the broad context of stochastic physics. Still we only aim at addressing the subgrid scale effects, which only constitute a subset of all the "hidden physics" considered in techniques such as the stochastically perturbed parametrized tendency (SPPT) (see [START_REF] Buizza | Stochastic representation of model uncertainties in the ecmwf ensemble prediction system[END_REF][START_REF] Palmer | Stochastic weather and climate models[END_REF][START_REF] Palmer | Stochastic physics and climate modelling[END_REF][START_REF] Shutts | Convective forcing fluctuations in a cloudresolving model: Relevance to the stochastic parameterization problem[END_REF] for general insight on this). Uncertainty attached to hidden physics include among others deep convection in ocean, ocean-atmosphere interaction, cloud convection scheme, boundary layer turbulence etc. For a broad view of the different stochastic frameworks proposed in climate sciences, the reader may refer to several recent reviews on the subject [START_REF] Berner | Stochastic parameterization: Toward a new view of weather and climate models[END_REF][START_REF] Franzke | Stochastic climate theory and modeling[END_REF][START_REF] Gottwald | Stochastic climate theory[END_REF][START_REF] Palmer | Stochastic weather and climate models[END_REF] The LU framework entails deep changes in the classical formalism. Among them, ensemble members are now semi-martingales and the material derivative is replaced by its stochastic equivalent when it comes to deal with transport equations. LU has the advantage to be flow-driven and physically relevant as it truly encodes transport through the stochastic flow and relies on physical conservation laws. This stochastic parametrization has been shown to produce a higher spread than perturbation of the initial condition while keeping a good trade-off with the error representation [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF][START_REF] Bauer | Stochastic representation of mesoscale eddy effects in coarse-resolution barotropic models[END_REF][START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarsegrid computational fluid dynamics[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, Part III SQG and frontal dynamics under strong turbulence conditions[END_REF]. A similar random decomposition is also at the center of the variational stochastic framework of [START_REF] Holm | Variational principles for stochastic fluid dynamics Proceedings of the Royal society A : Mathematical[END_REF]. As our setting, this latter approach constitutes a general formal framework to derive stochastic representation of fluid flow dynamics. Both frameworks rely on a stochastic transport principle. In [START_REF] Holm | Variational principles for stochastic fluid dynamics Proceedings of the Royal society A : Mathematical[END_REF], it is dedicated to Hamiltonian dynamical systems and defined from a circulation preserving constrained variational formulation, while the one in [START_REF] Mémin | Fluid flow dynamics under location uncertainty[END_REF] is based on Newtonian principles and built upon classical physical conservation laws. This latter scheme has been used as a fundamental tool to derive stochastic representations of large-scale geophysical dynamics [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF][START_REF] Bauer | Stochastic representation of mesoscale eddy effects in coarse-resolution barotropic models[END_REF][START_REF] Brecht | Rotating shallow water flow under location uncertainty with a structure-preserving discretization[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, Part I Random transport and general models[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, Part II Quasi-geostrophy and efficient ensemble spreading[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, Part III SQG and frontal dynamics under strong turbulence conditions[END_REF] or to define large eddy simulation models of turbulent flows [START_REF] Chandramouli | 4d large scale variational data assimilation of a turbulent flow with a dynamics error model[END_REF][START_REF] Harouna | Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling[END_REF] In this study, we would like to answer the following question: Does a stochastic parameterization enable us to replace in a more secure and efficient way the classical (multiplicative) inflation techniques? A general answer is obviously difficult to provide. However the experiments carried out here on the Surface Quasi-Geostrophic (SQG) model indicate that this effort seems worth being done. The SQG model is simple as it is a 2D model, yet it captures highly non linear realistic features of ocean surface dynamics [START_REF] Blumen | Uniform potential vorticity flow: part I. theory of wave interactions and two-dimensional turbulence[END_REF][START_REF] Constantin | New numerical results for the surface quasi-geostrophic equation[END_REF][START_REF] Constantin | Formation of strong fronts in the 2-D quasigeostrophic thermal active scalar[END_REF][START_REF] Constantin | Front formation in an active scalar equation[END_REF][START_REF] Held | Surface quasigeostrophic dynamics[END_REF][START_REF] Lapeyre | Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory[END_REF]. This model combines a buoyancy transport equation with a kinematic relation between buoyancy and velocity. For this example, we present numerical evidence that LU brings many improvements to classical ensemble simulation strategy [START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarsegrid computational fluid dynamics[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, Part II Quasi-geostrophy and efficient ensemble spreading[END_REF]. We show that this parametrization is also useful in a data assimilation context, again, by comparing to variance inflation. We present numerical evidence that the Mean Square Error performances are systematically better than all the deterministic inflation parameters tested here. We will also highlight the flaws of both inflation and localization parameter tuning when comparing them to a flow-driven stochastic model and point out that the latter seems to stabilize the localization drawbacks we observe in deterministic settings.

The outline of the paper is as follows. We first briefly present the Location Uncertainty principles and its mathematical implications. 1 Section 3 explains the numerical setup and details the results mentioned just above. The ensemble methods used in our experiments are recaped in appendix. The following is a summary of the numerical experiments carried out in this article :

• Numerical validation that the deterministic inflation parameter is difficult to tune , we show divergence of the filter in finite time (Figures 3 and4).

• Performance comparison between LU and deterministic settings both in terms of MSE (Figures 5 and6) and spread quality (Figure 7).

A conclusion with some perspectives closes the paper.

2 Transport under Location Uncertainty (LU)

Location Uncertainty models the impact of the small scales on the global flow dynamics. It is a stochastic framework that relies on the following decomposition of the Lagrangian velocity of a fluid particle positioned at x t in a spatial domain

Ω ⊂ R d : dx t = v(x t , t)dt + σ(x t , t)dB t , (1)
where v is a smooth-in-time component (referred to as the large-scale component in the following) and σdB t is a highly oscillating random component, built from a (cylindrical) Wiener process B t (ie a well-defined Brownian motion taking values in a functional space) [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF]. This latter component is uncorrelated in time. Let us point out that the above relation is rigorously understood in its integral form and corresponds thus only to a practical shortcut notation. The correlation operator σ is defined through a bounded matrix kernel σ, for any

function f ∈ (L 2 (Ω)) d , by σ(x, t)f = Ω σ(x, y, t)f (y)dy. ( 2
)
From this correlation operator, the covariance matrix kernel q reads q(x, y, t)

= Ω σ(x, x , t)σ(y, x , t) T dx , (3)
and the associated covariance operator Q is given by

Q(x, t)f = Ω q(x, y, t)f (y)dy. ( 4
)
The random velocity is Gaussian and distributed as

σdB t ∼ N (0, Qdt). (5)
Moreover, at each time t, the covariance operator Q(., t) is self-adjoint, nonnegative definite and compact. Thus it admits an orthonormal eigenfunction basis {φ n (•, t)} n∈N with non-negative eigenvalues (λ n (t)) n∈N . This entails an alternative convenient spectral definition of the noise as

σ(x, t)dB t = n∈N λ n (t)φ n (x, t)dβ n t , (6)
where the β n are i.i.d standard Brownian motions. From ( 6), the noise variance tensor a is then defined by

a(x, t) = n∈N λ n (t)φ n (x, t)φ n (x, t) T . (7)
It can be noticed the variance tensor has the physical dimension of a viscosity (ie m 2 /s). The properties and structure of the noise will of course depend on the procedure used to generate the orthonormal basis functions. The one used in our experiments will be presented later. In the deterministic case, a transported tracer Θ has zero material derivative :

D t Θ = ∂ t Θ + v • ∇Θ = 0.
In the LU framework, the material derivative is replaced by the stochastic transport operator

D t Θ = d t Θ + (v * dt + σdB t ) • ∇Θ - 1 2 ∇ • (a∇Θ)dt, (8)
where

• d t Θ = Θ(x, t + δt) -Θ(x, t)
is the forward time increment of the tracer.

• The effective advection velocity is defined by

v * = v - 1 2 ∇ • a. (9)
• The term σdB t • ∇Θ is a non-Gaussian multiplicative noise corresponding to the tracer's transport by the small-scale flow.

• The last term is a diffusion term, as the variance tensor a is definite positive.

Note that the expression of the transport operator is given here for a divergencefree noise. In the case of a compressible random field, the modified advection incorporates an additional term related to the noise divergence [START_REF] Resseguier | Geophysical flows under location uncertainty, Part I Random transport and general models[END_REF]. One key property of LU (for an incompressible random term) is that under the same ideal boundary conditions as in the deterministic case, it conserves the energy of the transported random tracer Θ:

d Ω Θ 2 (x)dx = 0, (10)
and, very importantly, this energy conservation property holds pathwise (i.e for any realization of the Brownian noise) [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, Part I Random transport and general models[END_REF]. A few general remarks on the stochastic transport operator can be done at this point. Compared to the usual material derivative, it encompasses additional new terms. The viscosity is now described by the variance tensor a. It is no longer a scalar variable but a (positive definite) diffusion matrix. It is thus richer than the classical eddy viscosity models [12] derived from an analogy with the molecular friction (the so called Boussinesq's assumption). The multiplicative noise corresponds to a backscattering of energy that is exactly compensated by the loss of energy brought by the stochastic diffusion (meaning the second order differential term associated to the stochastic transport operator in equation ( 8), which is indeed diffusive, in the sense that it dissipates energy). This balance is the reason why we have the energy conservation property and can be seen as a redistribution of the tracer (in the same way as a deterministic advection equation). Finally, the modified advection corresponds to the statistical effect induced by the smallscale inhomogeneity on the large scale component structuration [START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF]. With this term, the particles are statistically encouraged to migrate from regions with high variance (ie of high kinetic turbulent energy) to regions with low variance. This modification of advection, the backscattering carried by the advection noise and the balanced diffusion matrix are the three marked ingredients of the modeling under location uncertainty. It is important to point out that this balance only holds globally. Locally, these three terms may play their own role without any balance. This scheme has been used for the modeling of large scale flow dynamics for numerous flow configurations. From a data assimilation point of view, it has been used in an optimal control setting for 3D flow reconstruction and for joint parameter estimation in a 1D shallow water model [START_REF] Chandramouli | 4d large scale variational data assimilation of a turbulent flow with a dynamics error model[END_REF][START_REF] Yang | Estimation of physical parameters under location uncertainty using an ensemble 2 -expectation-maximization algorithm[END_REF], where it was shown to provide an interesting trade-off with the so called weak constraint variational assimilation.

The LU scheme obviously depends on the noise parametrization chosen. For instance, for a homogeneous noise associated to Fourier basis functions, the variance tensor is homogeneous and constant (even diagonal for a divergence free flow). Hence, there is no modified advection. For a stationary noise, the variance field is constant in time and thus not related to the evolving large-scale components. The ability to build a flow-dependent noise enables us to improve probabilistic forecasting skills [START_REF] Brecht | Rotating shallow water flow under location uncertainty with a structure-preserving discretization[END_REF]. For the SQG dynamics several noise parametrizations have been compared and assessed through different statistical proper scores [START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarsegrid computational fluid dynamics[END_REF]. One of the main findings of this work is that a time-varying inhomogeneous parametrization, termed as SVD (Singular Value Decomposition) noise, provides the best quantification results.

SVD noise : In this paragraph, we describe the generation of the noise used in the experiments. It relies on the creation of pseudo-observations at each point of the simulation grid, and then on the diagonalization of the associated empirical matrix to extract a proper basis to support the noise. Here the domain Ω is a two-dimensional regular grid of size n x × n y . The pseudo observations are constructed from the running velocity fluctuations around a sample mean, more precisely around a velocity field composed of local spatial means computed at each grid point. At each time t and each grid point x i,j , a spatial window W i,j of size n w × n w (with n w odd), much smaller than the whole simulation grid, is built around the point, together with the model boundary conditions (periodicity, replication,...) if the current point is on the border. Then a pseudoobservation is given by a draw of the velocity in the following set :

I(x i,j , t) = {v(x k,l , t) |k, l ∈ W i,j } . ( 11
)
Proceeding to n o draws within the window, and iterating over all the grid points, a global pseudo-observation matrix V is built :

V =    v 1 1 • • • v no 1 . . . • • • . . . v 1 nxny • • • v no nxny    , ( 12
)
whose size is (2 × n x × n y ) × n o (the 2 comes from the two components of the velocity v).

Then the mean over the n o pseudo-observations V is retrieved

V = V -V (13)
and Singular Value Decomposition is applied to the fluctuation matrix V to diagonalize the corresponding second order empirical moment. This way, the matrix Φ of the left eigenvectors on which we can decompose the noise as in ( 6) is obtained.

Let us denote the simulation grid scale and L = n w the spatial scale of the sliding window used to compute the noise. The procedure described above provides a noise σ L dB t at scale L, which is artificial and it must be downscaled to the true simulation scale . [START_REF] Harouna | Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling[END_REF] proposed a rescaling of the variance tensors based on 2D turbulent cascade assumption :

a = L 2 3
a L , which relies on an estimation of the velocity fluctuations at the simulation scale . Then the effective noise for the simulation grid is

σdB t = n -1 3 w σ L dB t . (14)
[48] also show that the Uncertainty Quantification (UQ) results for the SVD noise forecast reliability are good for all the metrics presented, in particular much better than the Perturbation on Initial Conditions (PIC), often used in ensemble data assimilation techniques. Some other procedures exist to generate the noise (6), relying for example on off-line learning on high-resolution datasets, or energy-budget based method (again see [START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarsegrid computational fluid dynamics[END_REF] for further details). This one has the advantage to be flow-driven and quite simple. As the basis depends both on time and on the ensemble member at hand, some adaptations should probably be devised for realistic models. One can envisage, for instance, a combination with wavelet basis or flow-based criterion to decide across time when the noise basis should be updated.

Experiments on both deterministic and stochastic models

This section details the numerical experiments carried out in this study. The main aim is to study the benefits brought by a stochastic dynamics in an up-todate version of ensemble Kalman filter with localization and inflation (details in appendix). In particular we wish to observe whether or not the stochastic dynamics brings by itself an efficient and practical inflation procedure for ensemble Kalman filtering.

In this work we apply the LU procedure to the 2D Surface Quasi-Geostrophic (SQG) model. This model constitutes an idealized dynamics for surface oceanic currents. Yet, it involves many realistic features of real world phenomena such as front, strong multiscale eddies, driven by a 3D like turbulent energy cascade, see [START_REF] Constantin | Front formation in an active scalar equation[END_REF][START_REF] Lapeyre | Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory[END_REF] for details.

The SQG model relies on a (deterministic or stochastic) transport equation of the buoyancy field b

D t b = 0, (15)
coupled with a kinematic equation

b = N strat f 0 (-∆) 1 2 ψ (16)
and the incompressibility constraints

v = ∇ ⊥ ψ; ∇ • σdB t = 0, (17)
linking the buoyancy field to the velocity field v, where ψ is the stream function, N strat is the stratification, f 0 is the Coriolis frequency and ∇f

= (∂ x f, ∂ y f ) T , ∇ ⊥ f = (-∂ y f, ∂ x f
) respectively stand for the classical and orthogonal gradients, while ∆ denotes the Laplace operator. The stochastic dynamics simply consist in replacing the material derivative D t by the stochastic transport operator D t given by [START_REF] Berner | Stochastic parameterization: Toward a new view of weather and climate models[END_REF]. We show in Figure 1 several realizations of the stochastic dynamics.

It can be noticed that the large-scale component of the different particles remain quite close after 17 days. They mainly differ by their small-scale features and vortices (subfigures (a,b,c,d)). We also added the states of these realizations at day 72 as an example to point out that the system is chaotic and may lead to significant large scale differences (subfigures (e,f,g,h)), at least on the time range studied in this work. This system is unforced and involves in its numerical implementation a small hyperviscosity term. It is hence decaying in the long run. However, as shown in Figure 1, the system remains turbulent at the end of the time period on which we will focus in this study. Mesoscale eddies as well as submesoscale eddies, filaments and fronts can be observed for the different realizations at day 72 of Figure 1. The stochastic and deterministic simulation is run on a simulation grid, G s , of size 64 × 64, meaning that each realization x n is a 64 × 64 matrix, or equivalently a vector of size 4096 (because here only the buoyancy field is observed and simulated, the velocity being obtained through the inversion of a fractional Laplacian). The actual physical size of the domain being 1000km × 1000km, two neighbor grid points are distant of about 15km. An observation y obs on a coarser observation grid, G o , of size 16×16, is assumed to be available every day (i.e. every 600 time steps of the dynamics). It is is generated as follows:

• A trajectory (z t ) t is run from the deterministic model (PDE) at a very fine resolution grid G f , of size 512 × 512

• Then a convolution-decimation procedure [START_REF] Lindeberg | Scale-Space Theory in Computer Vision[END_REF] is applied, which is the composition of a Gaussian filter G σ and a decimation operator D subsampling one pixel out of two. The Gaussian filtering writes

G σ (z t )(x) = (g σ * z t )(x),
where g σ is a two-dimensional Gaussian function. For any observable f defined on a grid (x 11 , ..., x 2n,2n ), we define Dg on a decimated subgrid (x 11 , ..., x n,n ) by

Df (x ij ) = f (x 2i,2j ).
This convolution-decimation operator D•G σ is then performed three times in order to fit to the targeted simulation grid.

• Finally, a projection operator P is applied from G s to G o , the latter being a subset of the first, and the observation is defined, for all t, by

y obs t = P • (D • G σ ) 3 (z t ) + η t , where η t ∼ N (0, R) and R = diag(r 2 1 , ..., r 2 M ) (18)
is the diagonal observation covariance matrix. It will often be considered that R = r 2 I M , where M is the number of points on the observation grid. As a consequence, the operator P plays the role of the observation operator H in the Kalman equations.

Numerical setup : The simulations have been performed with a pseudospectral code in space [START_REF] Resseguier | New trends in ensemble forecast strategy: uncertainty quantification for coarsegrid computational fluid dynamics[END_REF]. The time-scheme is a fourth-order Runge-Kutta scheme for the deterministic PDE, and an Euler-Maruyama scheme for the SPDEs [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF]. Regardless of the resolution and stochasticity, we use a standard hyperviscosity model to dissipate the energy at the resolution cut-off. The resulting implemented dynamics is :

D t b = α∆ 2 bdt, (19)
with a hyperviscosity coefficient α = (5 × 10 29 m 8 .s -1 )M -8

x , where M x is the grid resolution (here 512 for the fine-resolution PDE used to generate the observations, and 64 for the ensemble members). The boundary conditions are double-periodic. As mentioned before, the equations are mostly handled in the Fourier space, where the following SQG relation between velocity and buoyancy Fourier transforms can be used :

v = ik ⊥ b N strat ||k|| , ( 20
)
with k is the horizontal wave-vector, k ⊥ the orthogonal horizontal wave-vector and N strat is the stratification. The test case considered in this study is the following : an ensemble of N = 100 particles is started from the very same initial condition at day 0, which consists in two cold vortices to the north and two warm vortices to the south (figure 2) : The initial field is mathematically defined on each grid point r = x y by the following formula :

b 0 (r) = F (r -r 1 ) -F (r -r 2 ) -F (r -r 3 ) -F (r -r 4 ), (21)
where the vortices initial cores are

r 1 = 250km 250km ; r 2 = 750km 250km ; r 3 = 250km 750km ; r 4 = 750km 750km
and the function F is defined by

F (r) = B 0 exp - 1 2 x 2 σ 2 x + y 2 σ 2 y ( 22
)
with B 0 = 10 -3 m.s -1 , σ x = 67 km and σ y = 133 km. We also set the Coriolis frequency to f 0 = 1.028 × 10 -4 s -1 and the stratification to N strat = 3f 0 .

In this experiment, we study the differences of efficiency of the localized Ensemble Square Root Filter (cf appendix) with inflated deterministic forecast and non inflated stochastic simulations. In both cases, the initial ensemble is generated as follows. Starting from the initial condition, the stochastic dynamics is simulated using the SVD noise for 3 days (meaning 1800 time steps for the (S)PDE), without data assimilation. This way, a random ensemble is generated, and the performances of the SVD noise indicate that this ensemble is well-spread around the truth, which will be confirmed in Figure 7. An observation is provided each day (i.e. every 600 time steps of the SPDE (or PDE)), with an observation error covariance set to r i = 10 -5 for i = 1, ..., M in equation [START_REF] Constantin | New numerical results for the surface quasi-geostrophic equation[END_REF], which corresponds to a weak (but not negligible) noise on the observation. For the remaining of the simulation (100 days), this ensemble is used for two experiments :

• Experiment 1 (deterministic dynamics with inflation): The stochasticity is shut down after day 3, the forecast is then generated by the deterministic PDE, and prior inflation is used to artificially increase the ensemble spread. Namely, before the assimilation step, given a forecast ensemble x f , the inflated ensemble x f,inf is defined by

x f,inf j = x f + α(x f j -x f ), ( 23
)
where α > 1 is the inflation parameter. Then the localized ESRF is applied to x f,inf and the same procedure is done each day at each new observation.

• Experiment 2 (stochastic dynamics without inflation): The same SPDE drives the particles for the whole simulation, and the localized ESRF is applied with an observation each day.

The localization radius was set to l obs here, where l obs 60km denotes the distance between two neighboring observational sites, as it provided the best results for both cases. The inflation procedure is very sensitive to the parameter α. It must be finely tuned to have the best results. If it is too small, then the spread is not large enough. If it is too large, then it could give rise to a divergence of the filter (cf figure 3). For our SQG configuration, it turns out that the range of validity of α is between 1 and 1.08 approximately. Starting from α = 1.09, the filter starts diverging in the long term. In the context of our model, this small range for the inflation parameter and the tuning procedure is in itself a drawback of the inflation method. Still this tuning can be very different depending on the model at hand, so we do not claim that this range is small in general. Moreover, when α is badly chosen, it brings in the additional problem of strong divergence of particles of localized ESRF, as shown in Figure 4. In this example (r loc = 3l obs and α = 1.15), starting from day 50 or so, strong gradients between neighboring grid points can be observed, and progressively lead to non physical predictions. Obviously the same kind of behaviour also occurs in our setting with a smaller value of r loc . This was theoretically expected as, in the localization procedure, the posterior ensemble members are combinations by block of linear combinations of the prior ensemble members (cf equation [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equations[END_REF] in appendix). This formula relies on the strong assumption that such a block recombination remains a solution of the underlying PDE (or SPDE), which suggests that the global equation is equivalent to a combination of local ones. Let us note that for instance for fluid dynamics equations, the presence of a pressure term, solution of an elliptic problem, is intrinsically global and theoretically prevents such a local modeling. A careless application of localization may lead to the appearance of some discontinuities or gradients due to very different decisions taken by the filter at neighboring points. This may entail, after subsequent iterations of the underlying dynamics, major errors, numerical divergence of some ensemble members, leading eventually to unphysical realizations. The long-term instability of inflation is not a good sign, especially considering we deal here with a coarse-scale diffusive simulation. For finer resolutions with less diffusion and much more pronounced non-linear features, divergence is likely to occur sooner for the same inflation value. This instability is theoretically expected, as the artificial variance injection entailed by inflation is never counterbalanced in any way in the model. At a coarse resolution, even for low inflation parameters, the typical slow growth of error at the end of the simulation seems to indicate that long-time divergence as in figures 5 and 4 should be observed extending the simulation time (Figure 6 equivalently shows that increasing the inflation parameter makes the divergence occur earlier). The divergence caused by variance inflation can be attenuated by temporal adaptive schemes with introduction of diagnostic criteria [START_REF] Lee | Preventing catastrophic filter divergence using adaptive additive inflation for baroclinic turbulence[END_REF][START_REF] Raanes | Adaptive covariance inflation in the ensemble kalman filter by gaussian scale mixtures[END_REF]. Another alternative consists in changing the hypothesis on the prior distribution accounting for the sampling errors in the ensemble, which was shown in [START_REF] Bocquet | Ensemble kalman filtering without the intrinsic need for inflation[END_REF] to make multiplicative inflation optional on Lorenz models.

An increase of variance without controlling the global energy by a balanced dissipation raises the question of the mathematical well-posedness of the numerical scheme but also on the modified physics undergone by this forcing. These two questions are far from cosmetic. The first one is related to the generality of the numerical scheme (ie it must be valid for any noise and at any resolution, etc.). The second point is about the error terms, they should not change dramatically the targeted underlying physical system [START_REF] Chapron | Large-scale flows under location uncertainty: a consistent stochastic framework[END_REF]. These two points can hardly be met by an artificial increase of variance or by non physical multiplicative noise (see [START_REF] Chapron | Large-scale flows under location uncertainty: a consistent stochastic framework[END_REF] for an example on the Lorenz63 model).

On the opposite, the LU setting brings a natural balance between the energy brought by the noise and the amount dissipated by the stochastic diffusion. In addition, as shown in Figure 5, the simulation is stable while bringing an equivalent spread as a relatively strongly inflated deterministic simulation value (with respect to the SQG dynamics studied here). It leads also to better MSE results than the deterministic setting for the different values of inflation experimented here.

Figure 5: The LU framework (in red) performs better than the deterministic cases for all the reasonable inflation parameters tested and plotted in cyan (α = 1.09), green (α = 1.08), blue (α = 1.05) and black (no inflation: α = 1.0) color.

Figure 6 shows the same comparison with a greater inflation parameter set to α = 1.20. The time window is reduced to 50 days because, as expected, increasing the inflation parameter leads to an earlier divergence, here starting from days 40-50. Still it has comparable results with LU in the first 30-35 days. This means that LU has short-term MSE performances comparable with a very strong deterministic inflation parameter, and avoids long-term divergence as well. Additionally, we examine the spread of the ensemble members around the truth and observation points. We chose three characteristic grid points, corresponding to the center, north and southwest of the southwest warm vortex of the initial condition (Figure 2). In the following figure (Figure 7), we compare the behavior of the spread of LU ensemble members (in red) and deterministically inflated ones with parameter α = 1.10 (in magenta) and α = 1.20 (in orange) around the truth (blue dots) and observation points (black dots). Note that the two chosen inflation parameters are quite strong and lead both to divergence of the filter, with the second one exhibiting a divergence much sooner than the first one.

Figure 7 shows this comparison for the first 17 days of simulation. Although both spreads of subfigures (a) and (b) seem very similar for the figures of left and right columns, we observe on the center column figure that deterministic inflation does not provide enough spread to contain the truth, while LU does (see for instance time window 8,000-10,000). It can also be noticed that the truth at time step 8400 is completely skipped by magenta and orange spreads, while LU manages to reach it. We see that increasing the inflation parameter does not counteract the flaws of the smaller inflation parameter α = 1.10 (magenta spread, Figure 7), meaning that increasing α does not entail a richer ensemble. To that extent, we can also notice that when the stochasticity is shut down at day 3, the deterministic sets of trajectories immediately become less dense than the stochastic one. This means that the neighborhood of the truth/observation is visited much more often (in time) by stochastic trajectories. This improved "density" of the stochastic ensemble in the state space is an interesting feature, not exploited yet but it could be relevant in a particle filter framework.

Deterministic inflation sometimes offers more spread than LU. We observe this situation for example at the time window 9,000-10,000 on the right column figure. But this spread seems bigger than necessary, as LU has a smaller spread, but already well-centered around the truth. This supports the idea that deterministic inflation is blind, in the sense that it is decorrelated from the dynamics and the current state of the ensemble. LU does not only provide a spread of the ensemble equivalent to a strong inflation parameter, it brings an ensemble of better quality. This ensemble is constructed from the large-scale fluctuation and, in a way, fits the physics encoded by the original dynamics.

Another indicator of this is the following spread-error consistency in Figure 8. We compared the ensemble bias absolute value with the error estimated from the ensemble spread. We did this for LU and for the deterministic setting with inflation parameters α = 1.10 and α = 1.20 at day 17, when the three models are very close in terms of global MSE (see Figure 6). The results clearly point out that the LU framework provides much more spread and a better estimation of the error. The spread induced by inflation has much less physical structure and relevance.

Let us outline that the noise experimented here is free from any parameter.

Conclusion

We have shown in this study that the introduction of a stochastic parameterization of the flow dynamics within the location uncertainty setting enables to improve localized ESRF results when compared to an ensemble simulation of the original deterministic dynamics. Our main findings are as follows :

• The LU framework provides good spread, good MSE performances, and does not require variance inflation contrary to the deterministic ensemble simulations. This constitutes a first clear advantage of such a stochastic parametrization over an artificial increase of the ensemble spread. The parameter of such inflation is known to be difficult to fix and requires a fine tuning to get good performances.

• In addition, even for low values of the inflation parameter, we have shown that such a steady increase of variance leads to filter divergence at finite time. For finer resolutions with less diffusion and much more non-linear features, this time is likely to be much shorter for the same inflation value. Augmentation of variance through inflation corresponds to an error model with no particular spatial patterns. On the contrary, the stochastic dynamics consider a more subtle transport noise directly set from the largescale on-going simulation. The fact that the injection of variance is naturally counterbalanced by a diffusion term prevents long-term divergence.

It creates also a more physical and useful spread, whereas the spread entailed by deterministic inflation is sometimes too small and sometimes bigger than necessary.

• We also point out the potential drawbacks of localization techniques when accompanied with localization in the observation space, due to the block recombination of the solution. An augmentation of spread in a blind and constant way leads the local filter to take many different decisions, yielding a strong local augmentation of the gradients. When the dynamics involves an elliptic equation, as it is the case in the surface quasi-geostrophic model, this leads to dramatic divergence in finite time. • The localized Ensemble Square Root Filter with a stochastic dynamics seems to be able to handle sparse observations in time -in the experimental context used here only one observation field per day was availablemeaning every 600 time steps of the stochastic partial differential equation. This ability could be very helpful in the context of realistic forecasting of sea surface state. This study provides encouraging results and constitutes a first step toward this goal. Time interval of 3 to 10 days together with higher resolution observation will have to be investigated.

The stochastic SQG model experimented here constitutes an interesting model with highly non-linear features and a global elliptic equation. The resolution grid considered here is relatively small (ie 64 × 64) but yet of dimension larger than usual benchmarks elaborated from the different Lorenz models. We are confident the techniques could be extended to higher dimension.

As for the intrinsic nonlinear features of the dynamics, we intend to explore other data assimilation techniques like particle filtering for the stochastic SQG forecast. The particle filtering framework is promising as it constitutes a flexible data assimilation technique in terms of prior assumptions on both the linearity (or not) of the model and the Gaussianity of the errors. Even if it is known to struggle with high-dimensional systems, we expect the additional sophisticated methods devised in [START_REF] Cotter | Data assimilation for a quasi-geostrophic model with circulation-preserving stochastic transport noise[END_REF][START_REF] Cotter | A particle filter for stochastic advection by lie transport (salt): A case study for the damped and forced incompressible 2d euler equation[END_REF] to counteract the curse of dimensionality.

Classical EnKF (with perturbed observations) Starting from a forecast ensemble denoted x f , in the classical EnKF the ensemble is updated as follows :

• The so called Kalman gain matrix is defined by

K = P f H T (HP f H T + R) -1 . (25)
• The members of the ensemble are updated by :

x a n = x f n -K(Hx f n -y obs -ξ n ), (26)
with ξ n ∼ N (0, R), which corresponds to random draws of the observation error.

Ensemble Square Root Filter (ESRF) The square root scheme, developed in [START_REF] Whitaker | Ensemble data assimilation without perturbed observations[END_REF], provides a deterministic way to generate the posterior ensemble, without sampling of the ξ n s. The method is the following:

• The innovation term Hx f -y obs is multiplied by K in order to update the mean of the ensemble :

x a = x f -K(Hx f -y obs ). (27)
• The ensemble covariance matrix is updated as well :

P a = P f -KHP f . (28)
• Each member of the ensemble is given by

x a n = x a + (P a ) 1 2 (P f ) -1 2 (x f n -x f ). (29)
Since the ensemble covariance matrices are, in practice, of very high dimension, they cannot be neither stored nor manipulated (remember that for a state of dimension 10 11 , which is usual in weather forecasting, then this gives a covariance matrix of size 10 22 , that is to say of 10 7 peta (storage unit of the state variable)). An alternative is to rewrite the filtering equations only in terms of the ensemble anomaly matrix

A f = [(x f 1 -x f ), (x f 2 -x f ), . . . , (x f N -x f )],
which is a d × N matrix, so a more convenient one to store. In the following, we provide the expression of the different covariance matrices involved in this setting. An exhaustive presentation of the full computational details may be found in [START_REF] Reinhardt | Hybrid filters and multi-scale models[END_REF].

The ensemble covariance matrix is linked to the anomaly matrix by

P f = 1 N -1 A f (A f ) T ,
and the same relation holds for the posterior (or analysis) matrices. Injecting this relation in equation [START_REF] Held | Surface quasigeostrophic dynamics[END_REF] of the update of P a , one can find, after some algebraic manipulations, that

A a (A a ) T = A f I - 1 N -1 (HA f ) T (HP f H + R) -1 HA f (A f ) T .
So, if S is such that

SS T = I - 1 N -1 (HA f ) T (HP f H + R) -1 HA f , then A a = A f S.
Using the Sherman-Morrison-Woodbury formula [START_REF] Anderson | Optimal Filtering[END_REF], we deduce that

S = I + 1 N -1 (HA f ) T R -1 HA f -1 2 . ( 30
)
Let us remark that this matrix has size N × N , and Singular Value Decomposition can be used in order to compute the inverse square root. By definition of A a , the posterior forecast is given by

x a = x a + A a = x a + A f S. (31)
It can be shown that the mean x a can be written as a linear combination of the members of the forecast ensemble, namely

x a = N i=1 ω i x f i ,
where ω i is the i-th coordinate of the column vector

ω = 1 N 1 - 1 N -1 S 2 (A f ) T H T R -1 (Hx f -y obs ).
Finally, combining this with equation ( 31), the members of the posterior ensemble can be computed, for n = 1, . . . , N ,

x a n = N i=1 ω i x f i + N i=1 (x f i -x f )S in . ( 32
)
This filter is part of what is generally called Ensemble Transform Kalman Filter (ETKF), because members of the posterior ensemble are linear combinations of the forecast ones. Here Let us mention that such superposition principle (stability of solutions by linear combinations) is in theory only valid for linear dynamics. In the context of nonlinear dynamics, this setting has to be understood as constraining the solution to live in a particular Hilbert space of Gaussian processes spanned by the ensemble members and associated to a norm defined by the inverse ensemble covariance matrix.

Localized version of the ESRF

There exist two main ways to localize the filter. The R-localization (or local analysis) works on the observation error covariance matrix R. Another way to achieve localization is to work on the ensemble covariance matrix P (B-localization or covariance localization). As equations [START_REF] Houtekamer | Data assimilation using an ensemble kalman filter technique[END_REF][START_REF] Houtekamer | A sequential ensemble Kalman filter for atmospheric data assimilation[END_REF][START_REF] Houtekamer | Atmospheric data assimilation with an ensemble kalman filter: Results with real observations[END_REF] are only based on the use of the anomaly matrices instead of the whole covariance matrices P , R-localization is the only one compatible with our equations. [START_REF] Sakov | Relation between two common localisation methods for the enkf[END_REF] give details on both localization techniques and give numerical insights showing that both methods are likely to produce similar results. The basic idea of R-localization is to update the ensemble members, not in one step and using all the grid points, but coordinate by coordinate, using for each one a local neighborhood of "relevance". By doing so, the dimension of the state space is dramatically reduced for each local update. More precisely, for each grid point k = 1, . . . , d, let us denote p k the projection on the k-th coordinate (or grid point). For each k, the aim is to compute the local coefficients d k i,n such that

p k (x a n ) = N i=1 d k i,n p k (x f i ). (33)
In order to do this, R-localization is applied, meaning that the observation error covariance matrix will be modified for each coordinate, and then we will apply the previous formulas [START_REF] Houtekamer | Data assimilation using an ensemble kalman filter technique[END_REF][START_REF] Houtekamer | A sequential ensemble Kalman filter for atmospheric data assimilation[END_REF][START_REF] Houtekamer | Atmospheric data assimilation with an ensemble kalman filter: Results with real observations[END_REF] for this new matrix.

For each k = 1, . . . , d, a diagonal localization matrix C k of size p × p is defined, where p is the number of points on which you have observations (possibly less than d if the observation grid is coarser), and the diagonal coefficients are given, for l = 1, . . . , p, by

(C k ) ll = ρ d(l, k) r loc ,
where :

• d(l, k) denotes the distance between the grid point k and the observation site l.

• r loc is the localization radius.

Figure 1 :

 1 Figure 1: Four different realizations for the stochastic dynamics at day 17 and day 72.

Figure 2 :

 2 Figure 2: Initial condition of the buoyancy field for all particles.

Figure 3 :

 3 Figure 3: Mean Square Error curves for different values of the inflation parameter α in the deterministic case (case 1): magenta (α = 1.10), cyan (α = 1.09), green (α = 1.08), blue (α = 1.05) and black (no inflation: α = 1.0). Inflation is overall very beneficial compared to the black line (no inflation), but too big values of α lead to a long-term divergence of the filter (magenta and cyan).

Figure 4 :

 4 Figure 4: Examples of diverging realizations resulting from LESRF with r loc = 3l obs 180km and α = 1.15.

Figure 6 :

 6 Figure 6: Same figure as 5 with additional inflation parameters α = 1.10 (in magenta) and α = 1.20 (in orange). The time window is reduced to 50 days as the orange curve diverges much sooner than the others.

Figure 7 :

 7 Figure 7: Comparison of spread between LU and deterministic inflation on the first 17 days of the dynamics for three points located at the center, top and bottom left of the bottom left warm vortex of the initial condition. The upper row shows in red the buoyancy values at these points for the stochastic ensemble. The two lower ones show the buoyancy values at the same points for the deterministically inflated ensembles for α = 1.10 (in magenta) and α = 1.20 (in orange). The black dots are the observations; the blue line stand for the truth.

  Deterministic inflation α = 1.10 (c) Deterministic inflation α = 1.20

Figure 8 :

 8 Figure 8: Comparison between the ensemble bias absolute value e(x) = |b(x)b truth (x)| (left maps) compared to the estimate error (1.96× the standard deviation of the ensemble) evaluated at each grid point (right maps) at day 17. The upper row shows this comparison for LU, the other two show the same comparison with the deterministic setting respectively for inflation parameters of α = 1.10 (central row) and α = 1.20 (bottom row).

  n = ω i + S i,n -

For interested reader, thorough descriptions with detailed derivation of the stochastic representation of geophysical flow dynamics can be found in[START_REF] Bauer | Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach[END_REF][START_REF] Brecht | Rotating shallow water flow under location uncertainty with a structure-preserving discretization[END_REF][START_REF] Resseguier | Geophysical flows under location uncertainty, Part II Quasi-geostrophy and efficient ensemble spreading[END_REF].

Acknowledgments

The authors acknowledge the support of the ERC project 856408-STUOD

A Appendix: Ensemble methods

In this section we recall briefly the ensemble methods on which we focus in this study. We start by introducing the notations used here. Ensembles denoted x will contain N members. Each member lives in the state space R d . The empirical mean will be denoted by x and the (unbiased) ensemble covariance matrix is defined by

We are given an observation y obs living in R M with R its observation covariance matrix.

The observation operator H maps the state space R d to the observation space R M . Let us do a recap on different ensemble methods, including the classical EnKF and the localized square root filter used in the experiments.

• ρ is the Gaspari-Cohn function, defined for any z ≥ 0 by

It basically behaves like a Gaussian, but has compact support [0, 2], so that the coefficients are set to zero as soon as d(l, k) is twice as big as the localization radius.

Then the inverse observation error covariance matrix is modified by the Schur product

By doing so, the grid points outside this area are no longer taken into account, as the observation error is set to infinity at these points (due to the cancellation of spurious correlations). Finally, R is replaced by R k in formula [START_REF] Houtekamer | Data assimilation using an ensemble kalman filter technique[END_REF], which gives a localized version S k of S, and so a localized version of equation [START_REF] Houtekamer | Atmospheric data assimilation with an ensemble kalman filter: Results with real observations[END_REF].

In order to sum everything up, in this localized version, the update formula is given by