L Mathelin
email: mathelin@limsi.fr

C Desceliers

M Y Hussaini

Stochastic data assimilation of the random shallow water model loads with uncertain experimental measurements ⋆

Keywords: Parametric estimation, stochastic inverse method, Polynomial Chaos, data assimilation, maximum likelihood

This paper is concerned with the estimation of a parametric probabilistic model of the random displacement source field at the origin of seaquakes in a given region. The observation of the physical effects induced by statistically independent realizations of the seaquake random process is inherent with uncertainty in the measurements and a stochastic inverse method is proposed to identify each realization of the source field. A statistical reduction is performed to drastically lower the dimension of the space in which the random field is sought and one is left with a random vector to identify. An approximation of the vector components is determined using a Polynomial Chaos decomposition, solution of an optimality system to identify an optimal representation. A second order gradient-based optimization technique is used to efficiently estimate this statistical representation of the unknown source while accounting for the non-linear constraints in the model parameters. This methodology allows the uncertainty associated with the estimates to be quantified and avoids the need for repeatedly solving the forward model.

Introduction

Scientific context

It is no surprise that statistical estimation of stochastic processes has received considerable attention from researchers in the past and still continues to be a very active field as real-world problems are stochastic in nature. It is almost always necessary to quantify a stochastic process in the sense of giving it a tractable representation for further use in subsequent analysis or modeling. Such a representation allows one to reduce the stochastic process to a finite set of parameters describing it as closely as possible in a certain statistical sense. Indeed, techniques for propagating and quantifying uncertainty in numerical simulations are now widespread and they often require the description of the corresponding boundary and/or initial conditions. Their description may come from observational data of experiments, which leads to the problem of the identification of a stochastic process properties. In the general case, this stochastic process must be described in a specific space for its use as input data in a subsequent model. Often, one cannot directly observe the random process but one can only observe the effects it induces. The problem is then to identify the underlying stochastic process (in other words, infer its statistical properties) from a set of observational data.

In this work, we are interested in deriving an approximate probabilistic description of the sea bottom displacement field caused by earthquakes and modeled as a random process S(x, t; θ) ∈ Ω x × T × Θ, where Ω x is the spatial domain, T is the temporal domain and Θ is the space of elementary events θ. The stochastic framework is the probability space (Θ, σ Θ , P Θ ), where σ Θ is a σ-algebra and P Θ is a probability measure. The sea bottom motion field is different from event to event but the local geophysical properties cause the field to remain close to an "average" field. For instance, the amplitude and spatial extent of the motion field may vary but the global patterns remain similar because they are related to active tectonic plate locations, such as the contact line between two plates. A statistically well characterized approximation of the displacement random field can be subsequently used for simulating the uncertain future evolution of the tectonic plates or as an input parameter for inferring magma motions underneath the outer layer of the Earth. However, in the identification step of the sea bottom displacement stochastic field, the model we rely on to predict the output field (which is observed) from the input source field involves poorly known parameters and is then itself inherent with uncertainty. Further, the observation sensor on the ocean monitoring satellite is of finite accuracy and also introduces uncertainty in the measurements. Thus, we also aim to quantify the impact of the random parameters related to the physical model and sensor uncertainty on the estimate of the source field S(x, t; θ).

Among the usual methods, one of the most popular approaches is to rely on the maximum likelihood principle combined with an a priori functional form (e.g., Polynomial Chaos) for the source field probabilistic description (e.g., see [START_REF] Eggermont | Maximum Penalized Likelihood Estimation[END_REF], [START_REF] Ghanem | Polynomial chaos in stochastic finite elements[END_REF], [START_REF] Ghanem | Stochastic finite elements. A spectral approach[END_REF]). The probabilistic model of the source-field and all uncertain variables involved in the problem are then statistically sampled, the model is solved for the output field and the constitutive parameters of the source field are adjusted so that the likelihood of the observed data is maximized. This procedure leads to the optimal, in the likelihood sense, estimate of the source field in the chosen functional form. However, it does not allow one to provide confidence intervals to the identified form, nor can it give insights into the parameters most affecting the identification accuracy. Further, as the optimization step requires a large number of cost function evaluations, this approach is thought to become intractable in the general case.

Another popular approach for identifying the most probable parameter values of a system is the Bayesian approach. It relies on conceptual grounds similar to the maximum likelihood approach, but incorporates prior information in the parameters to estimate, usually leading to a better conditioned problem. Further, it results in probability distribution functions rather than in point-estimates. It has been successfully applied to various fields (see for instance [START_REF] Wang | A Bayesian inference approach to the inverse heat conduction problem[END_REF], [START_REF] Marzouk | Stochastic spectral methods for efficient Bayesian solution of inverse problems[END_REF], [START_REF] Marzouk | Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems[END_REF], [START_REF] Zabaras | A scalable framework for the solution of stochastic inverse problems using a sparse grid collocation approach[END_REF], [START_REF] Soize | Reduced chaos decomposition with random coefficients of vector-valued random variables and random fields[END_REF] or [START_REF] Koutsourelakis | A multi-resolution, non-parametric, Bayesian framework for identification of spatially-varying model parameters[END_REF]). However, it is not easy to incorporate several sources of uncertainty into the Bayesian framework, as necessary here. Further, the resulting identified field is often highly sensitive to the prior distribution, conflicting with the motivation of this work where we aim at providing confidence intervals for the identified source field.

In fact, there is a specific point to bear in mind: in the problem at hand, the uncertainty in the model and the sensor is epistemic as opposed to aleatory for the random process S(x, t; θ). This different nature of uncertainty allows one to separate their impact and one can make use of this to quantify the uncertainty associated with the identified stochastic process. A maximum likelihood-related approach does not account for this distinction and would result in a point-wise estimate. Further, the estimate would be erroneous because it is biased by the uncertainty in the model and the sensor.

Proposed method

In this paper, we propose to solve the stochastic data assimilation problem with the objective to identify a model of the stochastic input linked to the stochastic model of the physical problem at hand. This means that the stochastic model of the input will depend on the probabilistic models of both the uncertain constitutive parameters of the physical problem and the output (sensor uncertain measurements). Indeed, the poorly known parameters are parameterized with a set of known random variables ξ so that the probabilistic model of the identified source field is determined up to this ξ-parameterization. Similarly, the identified source field is determined up to the stochastic parameterization of the measurement model.

The identification of the source stochastic properties is achieved through a combination of inverse problems and statistical reduction. For the m-th realization θ m of a seismic event, the solution of an inverse problem allows one to estimate the deterministic source field S m ≡ S(x, t; θ m ) giving rise to the m-th realization of the observed data. A statistical reduction procedure is applied to the collection of identified source fields S m at a given time t m , m = 1, . . . , m obs , and allows one to approximate it in a Karhunen-Loève-like expansion form. One could think of letting the ξ-parameterization affect the empirical correlation matrix, which would then result in eigenvectors and eigenvalues directly parameterized in terms of ξ. However, this approach is hardly feasible since there is no guarantee of smoothness of the resulting eigenvectors and eigenvalues with respect to ξ. Instead, the parameterization is incorporated directly in the correlation matrix determination: the ξ-parameterization is described in a given functional form, say polynomial, so that it reduces to a finite set of variables (say, the polynomial coefficients). It results that the m-th identified source field is modeled as uncertain and corresponds to the uncertain input which, combined with uncertainty in the constitutive parameters of the physical model and uncertainty in the measurements, will result in the m-th measured data field. Upon describing these latter uncertainties with Polynomial Chaos (PC), the identified stochastic field is identified with a PC expansion and results in a vector-valued field of the PC coefficients.

The correlation matrix is then derived from the set of m obs vector-valued fields and the resulting eigenvectors lie in a Ω x × Ω ξ space. A functional form, e.g., Polynomial Chaos-based, is then prescribed for an approximation of the corresponding random variables, and the procedure is then conceptually similar to that of [START_REF] Desceliers | Maximum likelihood estimation of stochastic chaos representations from experimental data[END_REF] and [START_REF] Desceliers | Identification of chaos representations of elastic properties of random media using experimental vibration tests[END_REF] with a specific distinction being that no uncertainty was considered in the measurements in these papers. As a consequence, their approach was to maximize the likelihood of a deterministic field in contrast with the present work where uncertainty in the measurements and in some exogenous parameters leads to rely on a random field.

In the end, the identified probabilistic model of the source comprises both the epistemic and the aleatoric uncertainty. One is then able to estimate the uncertainty with which the stochastic process at the origin of the ocean surface motion is determined, hence defining the confidence one may have in the identified model based on the confidence one has in the physical model and the observational data. Once an optimal identification is found, the impact of the poorly known, ξ-parameterized, sensor and model parameters onto the identified source field probabilistic model can be quantified from the eigenvectors' dependence on ξ.

The paper is organized as follows. Section 2 presents the general seaquakes source problem considered in this paper for statistical identification. The different sources of uncertainty involved in the problem and their models are detailed in section 3. In section 4, the stochastic inverse problem is derived. It allows one to derive the source field corresponding to a given observed realization of the ocean surface field. This step is carried out for each realization m (seismic event). From a collection of m obs retrieved source fields (1 ≤ m ≤ m obs ), an optimization procedure is used for the final identification of the statistical source field properties and is presented in section 5. The numerical implementation of the solution method is given in section 6. Some results are finally shown in section 7 both for the stochastic inverse technique and the full stochastic identification problem. Final conclusions are drawn in section 8.

Shallow water model

Formulation

Due to tectonic forces, the ocean bottom is sometimes subjected to sudden motions when and where the material can not sustain the enormous pressure and slips and deforms to release the strain. These seaquakes induce ocean surface perturbations that can be monitored and subsequently used to estimate some parameters of the corresponding ocean bottom displacement field, here modeled as a source term S. The ocean surface behavior is assumed to be governed by the Shallow Water Equations, hereafter denoted SWE, and the bottom motion field is represented by three source terms denoted by S u , S v and S. The shallow water flow is then described by the following set of equations:

D u D t = f v -g G ∂w ∂x -b D u + S u , (1)
D v D t = -f u -g G ∂w ∂y -b D v + S v , (2)
∂w ∂t = - ∂ ((H + w) u) ∂x - ∂ ((H + w) v) ∂y + S, (3)
where D/Dt denotes the substantive derivative and the fluid density is assumed to be constant as well as the free surface pressure. Here, f is the term corresponding to the Coriolis force, b D is the viscous drag coefficient, u = {u, v} is the fluid velocity vector, w is the deviation of the ocean surface from its position at rest, g G is the gravity constant and H is the ocean depth field. Depending on the latitude Φ lat , one has f = 2 ω f sin(Φ lat ), where ω f is the angular rotational rate. Without loss of generality, it is chosen that the source field only acts along the w direction (S u = 0, S v = 0) and that the drag and the Coriolis forces can be neglected. Boundary conditions are prescribed along the edge of the domain

Ω x ⊂ R 2 as u • n = 0, ( 4
)
where n is the local normal vector to the boundary.

Solution of the shallow water model

The problem is discretized using a set of non-overlapping spectral elements and results in an N x -degree-of-freedom spatial representation for the ocean surface perturbation and the source field can be written as

S Nx (x, t) = Nx k=1 S k (t) e k (x), (5)
where e 1 (x), . . . , e Nx (x) are the usual finite element interpolation functions related to a finite element mesh of domain Ω x . More details about the numerical treatment of the SWE problem are provided in section 7.1. The above equation is rewritten as

S Nx (x, t) = e(x) T S(t), (6)
where e(x) ≡ (e 1 (x) . . . e Nx (x)) T and S(t) ≡ S 1 (t) . . . S Nx (t)

T .

(7)

The SWE solution w Nx (x, t) depends on x, t, S and H so that a functional g can be formally introduced, such as:

w Nx (x, t) = g(x, t; S, H). ( 8
)
3 Uncertainty sources

Ocean depth field

From the SWE, it is seen that the ocean depth field H directly impacts the surface perturbations, w. However, the local ocean depth field is poorly known and limited knowledge is available to describe it. It is modeled as an uncertain field depending on m H stochastic germs ξ 1 , . . . , ξ m H that are modeled as m H independent uniform random variables lying on [-1, 1]. Without loss of generality, one here takes m H = 1:

H(x; ξ 1 , . . . , ξ m H ) ≡ H(x) + c 0 ξ 1 H 1 (x), (9)
with H(x) the mean depth field and c 0 ≥ 0 drives the variance of H(x). More complex descriptions may be used such as a Karhunen-Loève decomposition if the covariance kernel is known.

Stochastic Shallow Water Model

To improve on the SWE model, intrinsic uncertainties related to the seaquake source have to be taken into account. As seen above, the goal is to identify the bottom displacement field statistics given only observations of the ocean surface. The representation of the source field is then modeled by a random field S indexed by Ω x × T written as

S Nx (x, t) = e(x) T S(t), (10)
where e(x) is the vector used in the representation introduced in section 2.2 and S is now an R Nx -valued stochastic process indexed by T .

Consequently, from Eq. ( 8), the representation of the ocean surface deviation from its position at rest is now a random field W Nx indexed by Ω x × T written as W Nx (x, t) = g(x, t; S, H).

(11)

Satellite measurements

The surface displacement of the ocean is modeled as a random field W indexed by Ω x × T and which is correlated with the ocean depth and the source of the seaquake. A tractable and practical way for getting information about the ocean bottom displacement is through observations of the ocean surface response to a seaquake event. The available observations comprise measurements of the height of a large area of the ocean surface in time as obtained from dedicated satellites such as Topex-Poseidon or the recently launched Jason-2. Let w 1 obs (x, t), . . . , w m obs obs (x, t) be measurements of ocean height at a given time t and at position x from m obs statistically independent seismic events. These measurements constitute the experimental database and are used to retrieve information about the source of the seaquakes. The set of observations is modeled as m obs statistically independent realizations of a random field W obs (S, H) indexed by Ω x × T . As with any measure, the surface observations are uncertain and it is assumed that the uncertainties related to the measurement process are known and that the random field W sat (S, H) can be written as

W obs (S, H) = W (S, H) + c sat ξ, ( 12
)
where ξ is a uniform random variable with values on [-1, 1] and the constant c sat represents the amount of uncertainty associated with the measures. Let the random field W sat (S, H, ξ) indexed by Ω x × T be the random observation of surface ocean displacement field constructed with the stochastic shallow water model. We then have

W sat (S, H, ξ) = g(x, t; S, H) + c sat ξ, (13)
4 Optimal random source field for each observation of the experimental database

Algebraic representation

We introduce the q-th order Polynomial Chaos expansion (hereafter denoted PC, [START_REF] Ghanem | Stochastic finite elements. A spectral approach[END_REF]) of S m (t) for all t ∈ T writing as

S m (t) ≈ q α,|α|=0 a m α (t) L α (ξ), (14)
where

ξ = ( ξ 1 , . . . , ξ m H , ξ) is an R m H +1 -valued vector of uniform random vari- ables on [-1, 1], α = (α 1 . . . α m H +1 ) with |α| ≡ m H +1 k=1 |α k | are the N m H +1 - valued multi-indexes and a m α (t) are the R Nx -valued coefficients of the PC ex- pansion of S m (t).
The polynomials L α are the multivariate Legendre polynomials such that

L α (ξ) = m H +1 j=1 L α j ({ξ} j ) , (15)
where L α j is a Legendre polynomial of degree α j . Then Eqs. ( 10) and ( 14) yields

S Nx (x, t) ≈ e(x) T q α,|α|=0
a m α (t) L α (ξ). ( 16)

Cost function for identification

For each experimental observation w m obs , m = 1, . . . , m obs , we introduce an optimal stochastic process S m,opt indexed by T , written as

S m,opt (t) = q α,|α|=0 a m,opt α (t) L α (ξ), ( 17
) such that W sat (S m,opt , H, ξ) -w m obs 2
S is minimal, where, for any second-order random field F indexed by Ω x × T , the norm • S is defined as

F 2 S = T Ωx E F (t, x)| 2 dx dt
and where E(•) is the mathematical expectation operator. Consequently, the vector a m,opt (t) = (a m,opt

α (t) , α ; |α| ≤ q) ∈ R N a chaos is solution of an optimization problem a m,opt = arg min a m ∈A W sat   q α,|α|=0 a m α L α (ξ), H, ξ   -w m obs 2 S , ( 18
)
where A is the set of all mappings a m from T to R N a chaos such that a m (t) = (a m α (t) , α ; |α| ≤ q). The length of a m,opt is N a chaos = (P + 1) N x where (P + 1) = (q + m H + 1)!/(q! (m H + 1)!) is the number of terms in the PC expansion, Eq. ( 14). Formulating Eq. ( 18) in words, one looks for the set of coefficients a m,opt such that the corresponding stochastic source field, combined with the uncertainty in the ocean depth field and the measurement, leads to the same (deterministic) ocean surface displacement field for all statistical realizations. Then, for each experimental observations, we construct the stochastic field S m,opt for m = 1, . . . , m obs , such that

S m,opt (x, t) = e(x) T q α,|α|=0 a m,opt α (t) L α (ξ). ( 19
)
The unconstrained optimization problem, Eq. ( 18), is here solved using a second-order quasi-Newton memory-limited technique, [START_REF] Gilbert | Some numerical experiments with variable-storage quasi-Newton algorithms[END_REF], relying on the gradient of the functional computed from the solution of the adjoint Shallow Water Equations. Iteratively solving the primal and dual SWE allows to converge to a m,opt within just a few tens of iterations.

No regularization term was found necessary thanks to the weakly non-linear character of the governing equations.

5 Global optimal random source field for the experimental database

Algebraic representation

For each observation w m obs , m = 1, . . . , m obs , the polynomial expansion of a random field S m,opt is constructed by solving Eq. ( 18). While the random fields S, S Nx and S m,opt are indexed on the time and space domain, we are only interested in describing their spatial structure at the time t m when their mean mechanical energy is maximum. For a given stochastic source S m,opt , t m is then

t m ≡ arg max t∈T E{S m,opt (x, t)} 2 Ωx . (20)
In the rest of the paper, we will hence drop the dependence in time and all time dependent quantities will be evaluated at t = t m .

In this section, a global optimal representation of the random fields S m,opt is constructed by considering each coefficient a m,opt as the m-th independent realization of a random vector A such that A ≡ (A α , α ; |α| ≤ q), where A α is an R Nx -valued random vector. We then define S opt as

S opt (x) = e(x) T q α,|α|=0 A α L α (ξ). (21)
In this section, an algebraic representation of random vector A α is constructed using Polynomial Chaos expansion. Since such a construction can be computationally expensive, a statistical reduction of the random vector A is first introduced. Different approaches may be thought of to derive a suitable representation of a random vector. Among them, one can cite Principal Component Analysis (also termed Karhunen-Loève), as will be considered in this work; factor analysis, useful when noise is present in the identification process; projection pursuit, which basically tries to find "meaningful" directions in the data [START_REF] Webb | Statistical pattern recognition[END_REF]); and Independent Component Analysis (ICA), where a basis leading to the least dependent components is searched for [START_REF] Hyvärinen | Survey on Independent Component Analysis[END_REF]). In the present case, the Karhunen-Loève approach is the method of choice since one is basically interested in reducing the size of the dataset to the most significant modes in the L 2 -sense rather than deriving a basis spanned by a potentially large number of independent random variables.

Let [C] be the (N a chaos × N a chaos ) positive-definite symmetric real matrix estimated for large value of m obs > N a chaos by

[C] ≃ 1 m obs m obs m=1 a m,opt -A a m,opt -A T , (22)
with

A ≡ 1 m obs m obs m=1 a m,opt . (23)
Let λ 1 ≥ . . . ≥ λ N a chaos > 0 and ψ 1 , . . . , ψ N a chaos respectively be the eigenvalues and the normalized eigenvectors of matrix [C] such that ψ 1 Ωx = . . . = ψ N a chaos Ωx = 1. Then, a N red -term statistically reduced representation of the random vector A is expressed in terms of the random vector

B T = (B 1 . . . B N red ) such that A ≈ A + [ψ] [Λ] B, (24)
where

[ψ] is an N a chaos × N red orthonormal matrix with [ψ] jk = {ψ j } k , [Λ] is a N red × N red matrix with [Λ] jk = √ λ k [I]
jk and where [I] is the identity matrix of R N red . We then have

B = [Λ] -1 [ψ] T (A -A) . ( 25
)
From the zero-mean orthonormal properties of [ψ], it follows that:

E (B) = 0, E B B T = [I]. ( 26
)
Since a m,opt is assumed to be the m-th realization of A, the realization B(θ m ) of B is such that, for all m = 1, . . . , m obs ,

B(θ m ) = [Λ] -1 [ψ] T a m,opt -A . ( 27
)
In order to derive a tractable statistical representation of B, it is decomposed into a PC expansion:

B ≈ B ≡ p β,|β|=0 b β H β (ζ), (28)
where, for a given ν ∈ N, ζ is an R ν -valued normal random variable of probability density p ζ , β T ≡ (β 1 . . . β ν ) with |β| = ν k=1 |β k | ≤ p are the N ν -valued multi-indexes and H β are the unit-norm multivariate Hermite polynomials, such that

H β (ζ) = ν j=1 H β j ({ζ} j ) , (29)
with H β j the Hermite polynomial of degree β j . Then, from Eqs. ( 26), it can be deduced that the constraints b β (t) must obey:

b 0 = 0, p β,|β|=0 b β b β T = [I]. ( 30
∈ R N b chaos , N b chaos ≡ N red (p + ν)!/ (p! ν!)
, such that the most likely (highest probability) approximated representation of B, Eq. ( 28), is that derived from the observations through Eq. ( 27). In other words,

b ≡ arg max b ′ ∈ B m obs m=1 p B (B m ; b ′ ) = arg max b ′ ∈ B L(b ′ ), (31)
which constitutes an optimization problem constrained by the fact that b lies on the manifold B defined by Eq. ( 30).

Finally, we can summarize the main equations at hand,

S opt (x) = e(x) T q α,|α|=0 A α L α (ξ), [Eq. (21)] , (32)
with A = (A α , |α| ≤ q) and 24) and (28)] ,

A ≈ A + [ψ] [Λ] p β,|β|=0 b β H β (ζ), [Eqs. (
where b ≡ (b β , β ; |β| ≤ p) is solution of the optimization problem defined by Eq. ( 31).

6 Statistical identification

Regularized optimization problem

The probability density function p B (B(θ m ); b ′ ) reads

p B (B(θ m ); b ′ ) = R ν δ 0   B(θ m ) - p β,|β|=0 b ′ β H β (s)   p ζ (s) ds, (34)
where δ 0 is the N red -dimensional Dirac distribution,

δ 0 (s) : R N red ∋ s → δ 0 (s) = N red i=1 δ(s i ) ∈ R. ( 35
)
As seen from Eqs. ( 31) and ( 34), the determination of the optimum b ′ implies the computation of ν-D integrals. The integrand is inexpensive to evaluate since it only involves basic operations on the PC series, and a Monte-Carlo method can be used to compute the integrals. For the subsequent use of gradient-based methods, a derivable integrand is desirable. Consider a surrogate delta operator defined as

δ 0 ( s) ≈ a a 2 + d( s) 2 n δ , d( s) ∈ R + , s ∈ R N red , (36)
with 0 < a ≪ 1, n δ ∈ N and where d( s) is the distance in R N red between B and B(θ m ). Note that this is in a similar spirit as the approximation of probability distribution functions using a kernel density estimation approach.

However, as the dimension N red of the random vector B grows, the support of the delta function δ 0 or its surrogate decreases exponentially and thus requires an increasing number of Monte-Carlo samples to accurately evaluate the integrals. This makes this approach difficult even for moderate length random vectors, say for N red 5. Effectively sampling such integrals has attracted a great deal of attention in the past two decades and usually relies on an adaptive defensive importance sampling based on some approximation of the integrand, see Owen & Zhou (1998) among others. Another approach is chosen here and a modified expression of the cost function is defined in Eq. ( 31) is considered:

log L ≡ - m obs m=1 log L m , L m ≡ N red j=1 p B j (B j (θ m ); b ′ ). ( 37
)
This alternative formulation prevents very low probabilities from having to be estimated as is the case when the support of the integrand decreases, thus softening the requirement on the number of Monte-Carlo samples necessary to estimate the integrals with a given accuracy. Of course, the optima from the two formulations of the cost function are the same when the random components B j , 1 ≤ j ≤ N red , are independent. This alternative definition of the likelihood will be used throughout the rest of the paper.

This kind of optimization problem is often solved with a simplex approach based on point estimations of the cost function for different sets of parameters b ′ . Here we take advantage of the regularity of the cost function expression to derive an efficient gradient-based second order optimization technique.

A potential issue is that the solution of the optimization problem lies in a high-dimensional R N b chaos space, bringing specific difficulties. In particular, the cost function response surface may exhibit local extrema and the optimization problem achieves only local convexity. The solution method relies on a combination of techniques. A good initial point candidate is given by matching the first statistical moments, [START_REF] Eggermont | Maximum Penalized Likelihood Estimation[END_REF]. However, this can only be done up to a limited order since it involves increasingly higher order inner products E H β i H β j H β k . . . of the PC basis, which are costly to evaluate. Further, the accuracy of the experimental high order statistical moments estimation is poor since the number m obs of observations is limited. The first four moments are considered and an initial point candidate is determined. Because the optimization space is of dimension N b chaos , which is much higher than four, we also rely on a collection of randomly chosen vectors on the constraint manifold B. From the best, in the sense of L, initial point found, a simulated annealing procedure is carried-out until it is reasonable to think that the basin of attraction of a good minimum (rather than a maximum, note the minus sign in Eq. ( 37)) has been reached. The optimization problem is then assumed locally strictly convex and a gradient-based approach is applied.

From Eqs. ( 34) and (37), the gradient of the cost function L can be computed as

∇ {b ′ β } j log L = m obs m=1 N red j=1 2 n δ a p B j (B j (θ m ); b ′ ) R ν d 2 n δ -1 j (s) H β (s) a 2 + d 2 n δ j (s) 2 p ζ (s) ds, (38)
with |β| ≤ p and

d 2 j (s) = B j -B j (θ m ) 2 . ( 39
)
To solve the optimization problem involving the non-linear constraint defining the manifold B, a SQP technique achieving quadratic convergence is employed. It relies on a memory-limited second-order quasi-Newton algorithm [START_REF] Gilbert | Some numerical experiments with variable-storage quasi-Newton algorithms[END_REF]) and estimates the Hessian matrix without computing it. The optimal point in R N a chaos is found through a multi-grid-like procedure as the PC order is initially set to a low value to find a coarse estimate of the stochastic source field. Once this estimate is found, a finer discretization is used from that initial point to determine the accurate source field. This approach reduces the risk of ending-up in a local, as opposed to a global, minimum of the Lagrangian and lightens the overall computational burden.

Algorithm

The global algorithm for solving the stochastic identification problem is summarized as follows:

(1) for each observation of the ocean surface field w m obs (x, t) after the mth seismic event, solve the corresponding optimal stochastic problem, Eq. ( 18), to determine a m,opt = a m,opt (t m ), ( 2) from the collection of m obs vectors a m,opt , 1 ≤ m ≤ m obs , carry-out a statistical reduction and determine B(θ m ), Eq. ( 27), (3) compute the cost function from Eq. ( 37) for different randomly chosen initial vector-valued functions b ′ onto the constraint manifold B and retain the best (lowest L) as the initial point, (4) from this selected initial point, determine the optimal set b, solution of Eq. ( 31), through a SQP-based solution method, (5) use Eqs. ( 32)-( 33) for constructing the optimal random field S opt .

Remark The approach followed in this paper involves two main steps. In a first step, a random field is constructed, through an inverse procedure, in order to represent the observed surface displacement for each 1, . . . , m obs observations (section 4, Eq. ( 18), step (1) of the algorithm). This first inverse problem requires the use of a mechanical solver (solving the governing equations of the physical problem at hand). The second step consists in constructing a unique random field for which the set of realizations constructed in step (1) would maximize the likelihood (section 5, Eq. ( 31), steps (2)-( 4)). In these last steps, the estimation of the likelihood does not require the use of the mechanical solver. The only difficulty lies on the high number of parameters required to represent the solution.

Alternatively, one could avoid this first step by directly searching for the parameters that would maximize the likelihood of the m obs observations. As for the proposed methodology, one of the difficulties would lie on the high number of parameters required to represent the solution. But, in such method, additional difficulties arise due to the necessary intensive use of the mechanical solver involved in the evaluation of the likelihood at each iteration of the optimization procedure. This approach is hence not thought to be tractable in the general case, unless the direct problem is very significantly cheaper than the inverse procedure and/or if the inverse step is overwhelmingly ill-posed.

7 Convergence analysis and results

Shallow Water numerical solver and solution method

To identify the seaquake source field, the inverse Shallow Water problem is solved. We rely on a 10 × 10-element mesh for the physical space, discretized with p x = p y = 6-th order Legendre cardinal polynomials for the fluid velocity. The unknowns are interpolated with these Legendre cardinal functions collocated at the Gauss-Lobatto points. The surface height, w, is discretized with a lower order polynomial to prevent spurious pressure modes from occurring (Q N -Q N -2 scheme), see [START_REF] Iskandarani | A staggered spectral element model with application to the oceanic shallow water equations[END_REF]. An additive Schwarz preconditioner is used to accelerate the solution process; see [START_REF] Douglas | An additive Schwarz preconditioner for the spectral element ocean model formulation of the shallow water equations[END_REF].

Solving the inverse SWE requires the solution of the adjoint equations which brings a potential issue since they involve the primal fields u, v and w. These fields may be large and their storage over the time horizon of the optimization procedure may introduce some difficulties. Different techniques exist to cope with this issue such as the check-pointing strategy. The interested reader may refer to [START_REF] Biros | Parallel Lagrange-Newton-Krylov-Schur methods for PDE-constrained optimization. Part I: The Krylov-Schur solver[END_REF] for further details. However, such an issue is not encountered in the present work due to its moderate size and the primal fields can be stored without specific problems.

As mentioned in section 3.1, m H = 1 was retained for describing the uncertain ocean depth field. The random vector ξ thus reduces to a random variable ξ and the random vector ξ then reduces to ξ ξ T . Other retained parameters include c sat = 10 -3 (Eq. ( 13)), q = 2 (Eq. ( 16)) and m obs = 1514.

Inverse problem verification for S m,opt

For the purpose of verification, let us first consider the case c 0 = 0 in Eq. ( 9) and c sat = 0 in Eq. ( 13). Determining S m,opt thus reduces to identifying the sea bottom displacement field giving rise to the observed surface perturbation: this is an inverse problem (IP) which solution is given by Eq. ( 18). This verification step is necessary since the IP is an essential ingredient of the whole statistic identification problem considered in this paper.

7.2.1

Identification of S m,opt for c sat = 0 and c 0 = 0

In this case, no uncertainty affects the physical model or the observation sensor and the inverse problem is deterministic. The surface is initially at rest and a time-varying source is applied at point x 0 = {x 0 , y 0 }:

S(x, t) = e -1 σ 2 S x-x 0 2 2 × t 2 1 + t 4 1 -e -20 t , (40)
with σ 2 S = 0.005 and x 0 = {0.4, 0.9}. All variables are non-dimensional and we are interested in the source field at the time of its maximum amplitude, corresponding to t = t m ≃ 4. Using the deterministic SWE model, the corresponding synthetic observation field w obs may be derived and mimics an observation w m obs that would have risen from an actual seismic event.

The inverse problem is initialized with an arbitrary source field at all time, as shown in Fig. 1. In general, one often has some expertise about the most likely location of the seaquake but the goal here is to investigate the efficiency of the inverse method and its robustness. The solution field of the IP is plotted in Fig. 2 (right) together with the true solution (left) and the agreement is seen to be very satisfying. No regularization terms were found necessary to achieve good identification results, even not a simple Tikhonov approach. This mainly comes from the fact that the problem at hand is only weakly non-linear and exhibits a one-to-one correspondence between the surface fluctuation w and the source S for a whole basin of considered parameters set. The problem is then bijective in a subdomain of these parameters' space. Further, the whole solution domain is observed and the null space of the inverse SWE operator is thus essentially not in the primal space. The inverse problem is then well-posed and essentially does not require regularization.

7.2.2

Identification of S m,opt for c sat = 0 and c 0 = 0 A more complex, and somehow more realistic, case is now considered for verifying the determination of the source field through a stochastic inverse problem (SIP): the ocean depth H is considered uncertain and the problem is then to determine a stochastic source field which, through the (uncertain) SWE equations, and given the uncertainty in the measure, would lead to surface data as close as possible to the observed (deterministic) data.

Relying on a similar synthetic observed field as for the previous test but centered around x 0 = {0.3, 0.9}, the mean of the identified source field (PC mode 0) at t = 4 is plotted in Fig. 3 together with the first, Fig. 4 (left), and second, Fig. 4 (right), stochastic modes of the identified stochastic source field decomposition. The initial field was the same as that plotted in Fig. 1 (deterministic). Again, the identified source is centered around the true solution. In this case, non-zero PC modes are also present as a consequence of the uncertainty of the problem.

The identified source field is stochastic instead of being deterministic, since, to lead to a given (deterministic quantity) surface observation at all times within the observation horizon T despite an uncertain depth field and observation device, the source term in the SWE equations must itself be uncertain. It represents the information one can derive on the source field from a given observation w m obs provided the uncertainty in the sensor is well characterized as well as the uncertain model describing the effects of the source field onto the ocean surface field.

The validity of the implementation is confirmed in Figs. 5, the resulting surface motion at a particular time (t = 300) is computed from the SWE equations with the actual, i.e., deterministic and given by Eq. ( 40), source field S(x, t) (left plots) and the identified source field issued from the SIP (right plots). Both surfaces are issued from the same stochastic code and only the source terms differ. In the left plots, the source term is deterministic and corresponds to the true, synthetic, source. On the right plots, the source term is that identified by the SIP. The resulting sea surface motion is seen to be similar in the statistical average sense (Fig. 5). However, the actual (deterministic) source term leads to an uncertain surface field (non-zero higher stochastic modes) while the identified (uncertain) source leads to a quasi deterministic surface (stochastic modes higher than 0 are essentially null), as shown in Figs. 6 and7. This clearly shows that a deterministic source field cannot induce a deterministic surface displacement field through the stochastic SWE model. This example and the essentially null uncertainty (quantified by the magnitude of the stochastic modes of non-zero index) in the predicted surface is a strong argument supporting the validity of the stochastic inverse problem strategy.

The subsequent global statistical identification of the seaquake random process can then reliably rely on the source estimation step. In particular, it shows that the uncertain source field, with first modes plotted in Figs. 3 and4, is indeed a versatile description of the source field which leads to a deterministic surface field given the uncertainty in the sensor and the physical model at hand.

Convergence analysis

Convergence analysis of the random vector B

In this section, the convergence of the stochastic identification process in terms of the Polynomial Chaos order p of the PC expansion of B, see Eq. ( 28), and stochastic germ dimension ν of ζ is investigated. For sake of clarity, we focus on the first random vector component B 1 . The fourth statistical moment for different PC orders p and a 2-D stochastic germ (ν = 2) is given in Table 1. The moment converges to an asymptotic value as the approximation relies on a finer discretization. The corresponding probability density function of the vector component is plotted in Fig. 8 for p = 5 and already exhibits a reasonable agreement with the pdf estimated from the data.

Identified source statistics

The mean-square convergence of the statistically reduced decomposition (Eq. 27) resulting from the identification procedure can be appreciated from Fig. 9 in terms of the eigenvalues spectrum. The spectrum decays quickly as higher index eigenvalues are considered. Since the variance of the approximation error σ 2 ε can be expressed as

σ 2 ε = N a chaos i=N red +1 λ i , (41)
with N red being the number of modes retained, the approximation of the identified sources rapidly converges in the mean-square sense as more modes are accounted for in the reduction. The approximation series can then be truncated to a limited number of terms. From the eigenvalue spectrum, it is seen that retaining only the first four modes (N red = 4) already accounts for most of the variance:

4 i=1 λ i N a chaos i=1 λ i ≃ 94 %. ( 42
)
Once the coefficients b β are found, we focus on the statistics of the identified source field, which is the main quantity of interest in this work. To this end, the auto-covariance function of the identified source random field is plotted in Fig. 10 together with that directly estimated from the collection of observations. The auto-covariance function writes

C S (x, x ′ ) ≡ E S opt (x) -E S opt (x) S opt (x ′ ) -E S opt (x ′ ) . (43)
It is a 4-D object and it is thus here plotted as a function of x only for the sake of clarity. As expected from Fig. 9, as more modes are accounted for, the agreement between the experimental and the identified random field quickly improves as a few modes are sufficient to reach a decent approximation.

The retained approximation parameters then include N red = 4, ν = 4, p = 5 and the identified stochastic source field can now be fully characterized. It is finally represented in Fig. 11 in terms of its 20% and 80% quantiles (nested plots). It is seen to be essentially located in a particular region of space, with a slight elongation along the x-axis, and exhibits a large uncertainty.

Since the identification procedure is subject to uncertainty (due to sensor and direct model uncertainty), the accuracy one can identify the underlying stochastic source field with is intrinsically limited. The method employed in this work allows one to distinguish between the intrinsic uncertainty of the physical system at hand (seaquakes) and that introduced by the identification procedure relying on limited knowledge of the physical system. As an illustration of the confidence one can have in the identified source description, the variance of S opt solely arising from the uncertainty introduced by the finite accuracy of the sensor and the poor knowledge of the ocean depth field, i.e., in the case E ζ 2 = 0, is plotted in Fig. 12. The uncertainty brought to the identified source field is seen to be reasonably low in the present case. Its spatial distribution resembles that of the source field, while exhibiting a more pronounced elongation along the x-axis. This plot shows how accurate and quantitative the estimation can be over the whole physical domain.

Concluding remarks

An efficient statistical estimation procedure has been proposed for use with realistic, uncertainty affected, data obtained from the indirect observations of statistically independent realizations of a random process that needs to be identified. The measurements are done on some quantity affected by the random process at hand and are subject to uncertainty arising from the experimental observation device and a poor knowledge of the physical model. For each statistical realization of the random process, a stochastic inverse problem is solved to evaluate the source field giving rise to that observation. A statistical description of the resulting set of source fields is then sought. A statistical reduction is performed to lower the dimensionality of the solution space and the random process to be identified is approximated under the form of a truncated Principal Component series. The resulting random vector is approximated with a Polynomial Chaos decomposition and an optimality principle, reminiscent of a maximum likelihood approach, is invoked to derive a good estimation of the PC development of the random vector components. To find the maximum of the manifold-constrained regularized cost function, a Sequential Quadratic Programming algorithm was used, providing a second order convergence rate and allowing for non-linear, equality and/or inequality, constraints on the parameters to identify.

The methodology derived in this work was applied to the estimation of the stochastic properties of a seaquakes source. Input data were available from uncertainty subjected satellite observations of the ocean surface after each seismic event and a statistical description of the ocean displacement field was determined. Results show that the statistical behavior of the identified random field (auto-covariance, probability density function) reproduces that given by empirical estimators from the experimental data.

Compared to the usual maximum likelihood-based method, and at the price of an additional effort in the development, the present inverse problem-based approach avoids the need for multiple resolutions of the forward model, a step that is CPU challenging for many problems of practical interest. Further, in addition to the estimation of the source field statistical properties, the precise uncertainty associated with this estimation is available and allows one to define, say, confidence bounds. Future developments of this methodology include improvements in the search of a global optimal estimation and an alternative, numerically more efficient, optimality principle.

  ) 5.2 Cost function for identification Let b = (b β , β ; |β| ≤ p) be the vector which consists of all the coefficients of the Polynomial Chaos expansion of B. Let b → p B j (b; b) be the probability density function of the j-th random component B j = { B} j . Then, the coefficients b β of the polynomial chaos expansion in Eq. (28) are identified as the vector b

Fig. 1 .Fig. 2 .

 12 Fig. 1. Chosen initial source field S(x, •).

Fig. 3 .Fig. 4 .

 34 Fig. 3. Plot of the mean identified source field (PC mode 0): S m,opt α

Fig. 5 .Fig. 6 .

 56 Fig. 5. Mean surface field at intermediate time (t = 300) for the actual (deterministic) source field (left) and the IP-computed source field (right). Both surfaces are issued from the same stochastic code and only the source terms differ.

Fig. 7 .

 7 Fig. 7. Stochastic mode 2 surface field at intermediate time (t = 300) for the actual source field (left) and the SIP-computed source field (right).

Fig. 8 .

 8 Fig. 8. Probability density function of the first random vector component B for a p = 5-PC order. ν = 2. The pdf estimated from data is plotted for comparison (histogram).

Fig. 9 .

 9 Fig.9. Spectrum of the identified sources covariance kernel. The eigenvalues decay quickly and most of the statistical variance is accounted for by the first few modes.

  Fig.10. Auto-covariance C S (x 1 , x 2 , y 1 , y 2 ) 1-D field for different numbers N red of modes in the reduced representation. From top left to bottom right, N red = 1, . . . , 4; ν = 2; p = 3. x 2 ≃ 0.38, y 1 = y 2 ≃ 0.62. The auto-covariance estimated from the experimental data is plotted with a solid black line.

Fig. 11 .

 11 Fig. 11. Confidence interval for S opt (x). The two nested surfaces plotted correspond to the 20% and 80% quantiles.

Fig. 12 .

 12 Fig. 12. Plot of the variance of S opt (x) in the case E ζ 2 = 0.

  Fourth statistical moment of the first random vector component m 1 4 ≡ E (B 1 ) 4 for different PC orders p. ν = 2, N red = 4.A similar study is carried-out for different germ dimensions ν at a given PC order (p = 5). Results are gathered in Table2in terms of the fourth statistical moment of the first component B 1 . Again, the statistical moment is seen to converge when the stochastic discretization is enriched enough. A germ dimension ν 4 leads to a reasonably converged estimation. Similar findings (not shown) hold for the other components of B. Fourth statistical moment of the first random vector component m 1 4 ≡ E (B 1 ) 4 for different PC dimension ν. p = 5, N red = 4.

	p	2	3	4	5	6	7	8	9
	m 1 4	13.09 13.23 13.19 13.26 13.25 13.31 13.32 13.32
	Table 1
		ν	1	2	3	4	5	6
		m 1 4	13.42 13.26 13.34 13.44 13.45 13.45
	Table 2

Acknowledgement

The first two authors would like to thank Christian Soize for fruitful discussions. LM also gratefully acknowledges the financial support of the Office of the Provost during his stay at Florida State University, the US NSF and the French National Agency for Research under projects ASRMEI JC08#375619 and TYCHE (ANR-2010-BLAN-0904).


F Arnal

R Dhaou

J Fasson

J Bernard
email: julien.bernard@toulouse.viveris.com

D Barvaux
email: didier.barvaux@toulouse.viveris.com

E Dubois

P Gélard

Handover management for hybrid satellite/terrestrial networks

Keywords: MIH, Satellite, 3GPP, LTE, Handover, Mobility, Test-bed

Initially envisaged to support handover between different wireless 802.x network technologies, the IEEE 802.21 standard also appears as the good candidate for handover management in future integrated satellite / terrestrial systems. This paper presents an analysis of how this standard could be implemented in the frame of a realistic scenario and taking into account the current trends in wireless network and mobility architectures. Our solution is then evaluated by means of emulation over a DVB-RCS representative testbed, and based on an experimental MIH implementation. We finally show that seamless handover can nearly be achieved with very short service outages.

Introduction

These last years show a tremendous evolution of GEO systems towards the integration of terminal mobility. This association with wireless terrestrial networks is a relatively new idea. Several proprietary hybrid systems integrating GEO satellites and ground components arise. The design of these integrated/hybrid systems takes into account physical, MAC, and network layers issues [START_REF] Kota | Satellite component of NGN: Integrated and hybrid networks[END_REF]. The combination of these two technologies can achieve for broadband services full coverage and high capacity. However the convergence at network level is still an open issue.

The main networking trial is certainly seamless handover. Managing handover in a couple of heterogeneous technologies may rely either on specific optimised or on standard and generic solution such as Media Independent Handover (MIH). The first solution may be optimised for each case of technologies. The last solution evolves more naturally and easily in order to integrate new access technologies. However there is a clear lack of experience in the deployment of this kind of solution. Therefore, our position here is to propose realistic hybrid network architecture and to see how MIH could be integrated.

The paper is organized as follow: Section 2 presents existing works related to mobility and MIH backgrounds. Section 3 identifies a possible use case for satellite/terrestrial integration and describes a relevant network architecture and design choices for the MIH usage. At last, Section 4 presents the evaluation of this solution through experimentations lead on an emulation platform.

Inter-system mobility architectures and MIH backgrounds

The major standard organisations support mobility features: MIH provides abstract services to the higher layers using a unified interface (L2.5 functionalities). MIH is not designed to take handover decision but rather provides the communication tools to assist the handover process. MIH defines three different services, i.e. Media Independent Event Service (MIES), Media Independent Command Service (MICS) and Media Independent Information Service (MIIS). MIES provides events triggered by changes in the link characteristic and status. MICS provides the upper layers necessary commands to manage and control the link behaviour to accomplish handover functions. MIIS provides information about the neighbouring networks and their capabilities. Extensive details on the IEEE 802.21 standard can be found in [START_REF]IEEE Standard for Local and metropolitan area networks-Part 21: Media Independent Handover Services[END_REF]. One major question is certainly the reason to use MIH in satellite context. Indeed, the needs of vertical handover schemes [START_REF] Márquez-Barja | An overview of vertical handover techniques: Algorithms, protocols and tools[END_REF] have been shown in several architectures based on the 802.21 standard [START_REF] Salhani | Terrestrial Wireless Networks and Satellite Systems Convergence[END_REF][START_REF] Qureshi | Adding Support for Satellite Interfaces to 802.21 Media Independent Handover[END_REF][START_REF] Hu | Mobility Extension for Broadband Satellite Multimedia[END_REF].

Existing studies on MIH, implementations and testbeds

To our knowledge, there is no reference implementation of the IEEE 802.21 standard. Nevertheless, FP6/FP7 IST projects have contributed to a prototype on GNU/Linux [START_REF] Piri | Towards a GNU/Linux IEEE 802.21 Implementation[END_REF]. Another prototype proposes an implementation of a mobile VoIP using SCTP with automatic address reconfiguration [START_REF] Chen | SCTP-based handoff based on MIH triggers information in campus networks[END_REF]. An experimental and analytical study was conducted [START_REF] Chiba | Performance Analysis of Next Generation Mobility Protocols for IMS/MMD Networks[END_REF] to determine the adequacy of protocols for mobility networks. This study presents the handover procedures (based on MIPv6 and PMIPv6) in WIFI/CDMA networks. PMIPv6 based on reactive handover method has been also observed in [START_REF] Choi | Seamless handover scheme for proxy mobile IPv6 using smart buffering[END_REF]. In [START_REF] Melia | Analysis of effect of mobile terminal speed on WLAN/3G vertical handovers[END_REF], a 802.21 client is simulated in a heterogeneous environment and the paper analyzed the effect of speed on terminal handovers between WiFi and 3G, but without validation of these simulations. Last, MIH has been proposed for future aeronautical communications [START_REF] Ayaz | Architecture of an IP-based Aeronautical Network to Integrate Satellite and Terrestrial Data Links[END_REF].

The use of MIH was also discussed in several prototypes. An implementation of 802.21 with a limited number of features was presented in [START_REF] Dutta | Seamless Handover across Heterogeneous Networks -An IEEE802.21 Centric Approach[END_REF]. In [START_REF] Buburuzan | Integration of Broadcast Technologies with Heterogeneous Networks -An IEEE 802.21 Centric Approach[END_REF] the integration of broadcast technology in heterogeneous networks using MIH is discussed. [START_REF] Tauil | Realization of IEEE 802.21 services and pre-authentication framework[END_REF] presents a solution for handover between heterogeneous networks that integrates MIH and Mobile Proxy Agent by defining an API performing the MIH Service Access Point. Eventually there is an extended research activity on seamless handover. Models of MIH are integrated to simulation tools such as Qualnet [START_REF] Latkoski | SDL+QualNet: A Novel Simulation Environment for Wireless Heterogeneous Networks[END_REF] and ns-2 [START_REF] Rouil | MIH transport using cross-layer optimized stream control transmission protocol[END_REF]. Open source implementations such as openMIH [START_REF] Lopez | OpenMIH, an Open-Source MIH Implementation and Its Application to Proactive pre-Authentication[END_REF] and Odtone [START_REF] Corujo | Using an open-source IEEE 802.21 implementation for network-based localized mobility management[END_REF] have been also set up.

Scenario description

Mobile Services A hybrid system aims at offering access to users' services wherever they are by enlarging the system coverage. However, such systems only make sense if they offer a mobile service. For a single user to a whole group of users in the same vehicle, applications are miscellaneous, from voice applications to Internet access.

Access networks

Between the different types of satellites, GEO systems offer the best opportunity to cope with the terrestrial part in terms of coverage and cost. Apart private solutions, the main standards proposing mobility for GEO are GMR, the mobile version of DVB Return Channel for satellites (DVB-RCS+M) and DVB for Satellite Handheld (DVB-SH) that focus on mobile broadcast services. GMR technology has been designed for extension of the cellular world to satellite and thus it focuses on voice services. Some integrated systems, like Terrestar [START_REF] Vojcic | Network of Mobile Networks; Hybrid Terrestrial-Satellite Radio[END_REF] are nowadays using GMR technology to interconnect it to a mobile operator core network. Indeed GMR standard has evolved, following the "all-IP" wave. Therefore a LTE solution adapted to satellite appears as the right direction for GMR standard. In Ku/Ka bands, the use of DVB-RCS+M is highly relevant for vehicular services (i.e. mobile collective terminal) and this will be retained for our reference scenario, including for the test-bed implementation.

Concerning the terrestrial networks, we could consider any kind of terrestrial networks that may provide mobility. If LTE and Wimax are two promising solution, the limits of test-bed have induced the use of a simpler technology, WIFI.

Network architecture

A major element in network architecture is the different roles undertaken. We assume a mobile service provider (MSP) ensures the IP connectivity of its mobile clients (final users and/or hybrid terminal) and may provide some other kinds of services. This role is generally linked to the core network of final user. MSP is in charge of interaction with the manager of the different access networks that its clients are using. If needed, it may rent bandwidth at user's demand. The following Fig. 1 illustrates a distribution of roles through a DVB-RCS+M / LTE architecture. A group of users are connecting through a LAN in a train to their services. The hybrid terminal can connect to an E-nodeB or to the satellite Gateway. Finally, a Mobile IP / NEMO architecture is chosen for network mobility support.

The terrestrial and satellite networks are access networks between the user and his MSP. The MSP may propose the first visible router to user and it manages the AAA functionalities. However, some of these functionalities can be delegated to access network, using MSP Point of Presence (PoP). Moreover, MSP role can directly be undergone by one of the two access networks, especially for the 3GPP one. This option is more convenient for LTE proposing a whole set of services like IP Multimedia Subsystem (IMS) as shown on the figure. For satellite systems or Wimax the core network is less developed than cellular one. So this integration demands more change to their network architectures.

Policing functionalities can be distributing on both sides and centralized by a global AAA entity in the MSP. In the case of a LTE network this function is managed by the Policy and Charging Resource Function (PCRF) whereas for others access networks the control at IP level can be made in implementing the IETF COPS standard. Policy Enforcement Point (PEP) and Policy Decision Point (PDP) are then introduced in the architecture, the former being deployed in the access network while the later is directly managed by the Mobile Service Provider (MSP). The interaction between these functionalities and the handover signaling is however out of the scope of this paper.

Handover decision

MIH supports both mobile-initiated handover as well as network-initiated handover, both alternatives can be preferred according to the exact deployment. In our scenario, the terminal chooses its mode of connection and the timing to perform the handover. A simple decision algorithm is implemented as an independent process (as a MIH user) in the terminal and based on predefined Signal-To-Noise Ratio thresholds. Thanks to this modular architecture, the handover decision and/or policies is open and could be easily extended with much sophisticated algorithms.

MIH functionalities

The MIH Function (MIHF) aims at referencing any event and at informing its users. In our proposal, MIHF is present in mobile and network elements, as it is the most general solution. For Wimax, MIHF is integrated on the access points whereas for LTE it may be the PDN-GW. However if the MSP takes part in the handover decision, it must have an MIHF so as to communicate with the handover decision function. MIHF is communicating with other MIHF in the different networks and they collect information, inform or receive order from the MIH users on the same entities. Another interesting MIH user is the MIIS server which can be used to select the available networks. However this option has not been chosen here.

MIH primitive selection

Many primitives are defined in the MIH standard although all are not necessarily required to be implemented. The selection of the right primitives may allow a seamless handover. The selection has been lead by the need to obtain a seamless handover and the convenience of the deployment for the tests. The following primitives have been used for handover indication:

• MIH_Link_Configure_Thresholds is used to fix the thresholds for signaling the need of a handover before the link is going done.

•

MIH_Link_Parameters_Report proposes to inform the different entities when the threshold is reached. It allows the preparation of the handover.

•

MIH_Link_Going_Down is used by the L2 interface to announce the near link outage.

Implementation

Description of the emulation test-bed, implementation and scenario

The core platform used is the Platine emulation tool [START_REF] Baudoin | Overview of Platine emulation testbed and its utilization to support DVB-RCS/S2 evolutions[END_REF], representative of a DVB-S2/RCS access network and running under a Linux environment. NEMO was integrated to Platine from the UMIP v0.4 software upgraded with the NEPL patch [22]. Standard Linux Route Advertisements (RADVD) are used for address autoconfiguration. Finally the ODTONE v0.2 beta stack was used for MIH implementation, modified by minor adaptations for our needs. In order to model the LTE/3GPP connection we have used a WIFI link which was more convenient for performing tests in our lab environment. In addition, only the handover from WIFI to satellite could have been tested, as few extra developments were further needed for Satellite to WIFI handover in our platform.

The (emulated) mobile network is connected through a hybrid Satellite Terminal that can be connected to the satellite gateway through the emulated satellite link, or to the WIFI access point. The NEMO instance and the MIH local function are therefore deployed on the Satellite Terminal, while the satellite and WIFI gateways host the network-side MIH Functions. . Each threshold is mapped to a given MIH event (LINK_GOING_DOWN and LINK_PARAMETER_REPORT). The first event, when received, is used to start the Platine logging process. If successful, Route Advertisements can be received for the satellite connection through the virtual Platine network interface. Subsequently, when the second threshold is crossed, it triggers the NEMO handover phase. We capture this event to send an inter-process signal (SIGUSR2) to the NEMO stack in order it performs the Binding Update / Acknowledgement procedure. For future improvements, this interworking will be developed in a proper module that will also implement separately the Decision algorithm for flexible evolutions.

Application flows are composed of a CBR audio and a VBR video stream (RTP/UDP), each generated by VLC at a host located in the mobile network, and sent to a Correspondent Node (CN). This receiver is located at a neighboring network interconnected to both WIFI and satellite networks.

Results

A Wireshark client running at the CN measures the incoming rates. As shown in Fig. 3, a short outage (0.6 seconds) is experienced just after the handover at second 90. At the receiver, the audio flow starts to be corrupted, but the video quality remains intact.

Fig. 3. Incoming bit rates variations measured at Correspondent Node

Our investigation shows that no packet loss happens during the handover, which is a first point of satisfaction. With comparison to a basic scenario test conducted in the same conditions but without MIH, that cannot be presented here due the lack of space, a packet loss sequence lasting nearly for 4.5 seconds was observed.

Immediately after the handover, we understood in observing an increase of the satellite queue utilization that the introduced delay is not related to any mobility latency but to the satellite access layer allocations loop (Rate-Based Dynamic Capacity allocations were used here). In addition, the sudden propagation delay change (tenths of milliseconds with WIFI, 300ms with satellite) also adds a contribution to the observed delay. However since the video client largely buffers the receiving flow, this latency effect remains transparent, whereas audio samples are by far less buffered.

For the remaining period, audio remains corrupted and not video. This is due to the absence of any QoS support on the overloaded satellite connection. Buffering at the satellite link layer continues is never completely absorbed. Except this effect, we conclude this MIH integration is successful and achieves good performances, though only based on an experimental stack.

Conclusion

This paper has demonstrated that MIH, that appear complex at first glance, is actually well applicable to any hybrid IP satellite network, such as the future LightSquared system in the US based on LTE. Coupled with Mobile IP network mobility, it is able to show very good performances provided that the loss of signal can be detected before the link completely goes down. In our future works, derived scenarios with the support of QoS will be presented in an IMS environment, in order to illustrate how resources could be provisioned before the mobile switches to another network. We also intend to complete our integration work by supporting Satellite to WiFi handover, and maybe to implement handover decisions based on real SNR and/or location measurements.

Fig. 1 .

 1 Fig. 1. An example of network architecture for the DVB-RCS+M/LTE reference scenario

Fig. 2 .

 2 Fig. 2. MIH sequence diagram between the MIH user, MIH functionality and test-bed Signal-to-Noise Ratio (SNR) variations are emulated by means of a pre-computed scenario stored in a file indicating the virtual SNR measurements. They are polled each 5 seconds by the mobile MIH Function (LINK_SAP interface). Two SNR thresholds are configured at the MIH stack (at the MIH startup, with

  • IETF proposes network, transport and session level mobility solutions, therefore usable for inter-system handover. The MIPSHOP and MEXT working groups focus on mobility issues in heterogeneous networks. Several mobile IP based solutions and extensions were proposed and a big trend has been raised on the development of "Smart Routers" products that are based on such solutions. IEEE has formed the 802.21 group, where MIH was formerly proposed for handover between 802.x access technologies. Later, the standard has been extended for 3GPP networks support.

Transport level mobility solutions are based either on TCP (MP-TCP, Freeze-TCP and Mobile Socket Service) or on SCTP. SIP-based mobility was also proposed. • 3GPP/3GPP2 architecture not only supports inter-technology mobility but also vertical handover (e.g. with WiFi and Wimax). Two operation modes are defined: the non-optimized and the optimized handovers. For instance, VoIP services, in a 2G/3G network, are based on the first mode whereas no standardized solutions exist for the second one. The network level mobility is based on MIPv4, DSMIPv6 and PMIP (Proxy MIP). • From another hand,


Didier Auroux
email: auroux@unice.fr

Jacques Blum

The "Back and Forth Nudging" algorithm for oceanographic data assimilation

Introduction

In this paper, we generalize the so-called "nudging" algorithm in order to identify the initial condition of a dynamical system from experimental observations. The standard nudging algorithm consists in adding to the state equations of a dynamical system a feedback term, which is proportional to the difference between the observation and its equivalent quantity computed by the resolution of the state equations.

The Back and Forth Nudging algorithm consists in solving first the forward nudging equation and then a backward equation, where the feedback term which is added to the state equations has an opposite sign to the one introduced in the forward equation. The initial state of this backward resolution is the final state obtained by the standard nudging method. After resolution of this backward equation, one obtains an estimate of the initial state of the system. We iteratively repeat these forward and backward resolutions (with the feedback terms) until convergence of the algorithm. This algorithm is presented in section 2. Then, we present in section 3 a way to choose the forward and backward nudging matrices in order to have well-posed problems, and to obtain the convergence of the algorithm. Then we report in section 4 the results of numerical experiments on several oceanographic models, including comparison and hybridization with the 4D-VAR algorithm. Finally, some conclusions are given in section 5.

Back and Forth Nudging (BFN) algorithm

In order to simplify the notations, we assume that the model equations have been discretized in space by a finite difference, finite element, or spectral discretization method. The time continuous model satisfies dynamical equations of the form:

, 0 ), ( T t X F dt dX < < = (1)
with an initial condition . In this equation, represents all the linear or nonlinear operators of the model equation, including the spatial differential operators. We will denote by

0 ) 0 ( x X = F
H the observation operator, allowing us to compare the observations with the corresponding , deduced from the state vector . We do not particularly assume obs X

))

( ( t X H ) (t X
H to be a linear operator.

If we apply nudging to the model (1), we obtain

( ) , 0 , ) ( ) ( T t X H X K X F dt dX obs < < - + = (2)
with the same initial condition, and where K is the nudging (or gain) matrix. Note that it may also be a nudging scalar coefficient in some simple cases. The model then appears as a weak constraint, and the nudging term forces the state variables to fit as well as possible to the observations. In the linear case (where is a matrix, and F H is a linear operator), the forward nudging method is nothing else than the Luenberger observer [START_REF] Luenberger | Observers for multivariable systems[END_REF], also called asymptotic observer, where the matrix K can be chosen so that the error goes to zero when time goes to infinity. Unfortunately, in most geophysical applications, the assimilation period is not long enough to have the nudging method giving good results. Then, one can apply the nudging technique to the model, but backwards in time, in order to come back to the initial state.

The back and forth nudging algorithm, see e.g. [START_REF] Auroux | A nudging-based data assimilation method for oceanographic problems: the Back and Forth Nudging (BFN) algorithm[END_REF], consists of first solving the forward nudging equation and then the backward nudging equation (i.e. the model equation backwards in time, with a nudging feedback to the observations). The initial condition of the backward integration is the final state obtained after integration of the forward nudging equation. At the end of this process, one obtains an estimate of the initial state of the system. We repeat these forward and backward integrations (with the feedback terms) until convergence of the algorithm:

( ) ( ) ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = > > - - = = < < - + = ≥ - ), ( ) ( , 0 , ) ( ' ) ( ~), 0 ( ) 0 ( , 0 , ) ( ) ( 1 1 T X T X t T X H X K X F dt X d X X T t X H X K X F dt dX k k k k obs k k k k k obs k k (3) with the notation 0 0 ) 0 ( ~x X = .
Note that this algorithm can be compared with the 4D-VAR algorithm (Le [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF], which also consists in a sequence of forward and backward resolutions. But in our algorithm, it is useless to linearize (or differentiate) the system and the backward system is not the adjoint equation but the model equation, with an extra feedback term that stabilizes the ill-posed backward resolution.

Choice of the nudging matrices and interpretation

Variational interpretation of the nudging

The standard nudging method has been widely studied in the past decades (see e.g. [START_REF] Bao | An adjoint examination of a nudging method for data assimilation[END_REF]). Thus, there are several ways to choose the nudging matrix K in the forward part of the algorithm. One can for example consider the optimal nudging matrix , as discussed in [START_REF] Vidard | Determination of optimal nudging coefficients[END_REF]. In such an approach, a variational data assimilation scheme is used in a parameter estimation mode to determine the optimal nudging coefficients. This choice theoretically provides the best results for the forward part of the BFN scheme, but the computation of the optimal gain matrix is costly.

opt K When , the forward nudging problem (2) simply becomes the direct model (1). On the other hand, setting

0 = K +∞ = K
forces the state variables to be equal to the observations at discrete times, as is done in [START_REF] Talagrand | A study of the dynamics of four-dimensional data assimilation[END_REF]. These two choices have the common drawback of considering only one of the two sources of information (model and data).

Let us assume that we know the statistics of errors on observations, and denote by R the covariance matrix of observation errors. This matrix is involved in all standard data assimilation, either variational or sequential. Usually, it is impossible to know the exact statistics of errors, and thus only an approximation of R is available, assumed to be symmetric positive definite.

We assume here the direct model to be linear (or linearized). We consider a temporal discretization of the forward nudging problem (2), using for example an implicit scheme. If we denote by n X the solution at time and

n t 1 + n X
the solution at time , and , then equation ( 2) is equivalent to the following optimization problem:

1 + n t n n t t t - = Δ +1 ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - - Δ + Δ - - - = - + HX X HX X R t X FX t X X X X X obs obs n n X n ), ( 2 , 2 , 2 1 min arg 1 1 (4) if
is symmetric and if the nudging matrix is set to be F

1 - = R H K T .
The first two terms correspond exactly to the energy of the discretized direct model, and the last term is the observation part of the variational cost function. This variational principle shows that at each time step, the nudging state is a compromise between minimizing the energy of the system and the distance to the observations. As a consequence, there is no need to consider an additional term ensuring an initial condition close to the background state like in variational algorithms, neither for stabilizing or regularizing the problem, nor from a physical point of view. One can simply initialize the BFN scheme with the background state, without any information on its statistics of errors.

Pole assignment method and backward nudging matrix

The goal of the backward nudging term is both to have a backward data assimilation system and to stabilize the integration of the backward system, as this system is usually ill posed. The choice of the backward nudging matrix is then imposed by this stability condition.

If we consider a linearized situation, in which the system and observation operators ( and F H , respectively) are linear, and if we make the change of time variable , then the backward equation can be rewritten as

t T t - = ' ( ) . ' 0 , ' ' T t HX X K FX dt dX obs < < - - = - (5)
Then, the matrix to be stabilized is

H K F ' - -
, i.e. the eigenvalues of this matrix should have negative real parts. The pole assignment result (see e.g. [START_REF] Trélat | Contrôle optimal: théorie et applications[END_REF]) claims that if the system is observable, then there exists at least one matrix such that

' K H K F ' - -
is a Hurwitz matrix, i.e. all its eigenvalues are in the negative halfplane. However, such a matrix may be hard to compute, as it usually requires the resolution of a Riccati equation.

Numerical experiments

Numerical choice of the nudging matrices

All numerical experiments have been performed with an easy-to-implement nudging matrix:

where is a positive scalar gain, and ,

1 - = R kH K T k
R is the covariance Paper.273.3 matrix of observation errors. This choice does not require a costly numerical integration of a parameter estimation problem for the determination of the optimal coefficients. Choosing

L H K T =
, where is a square matrix in the observation space, has another interesting property: if the observations are not located at a model grid point, or are a function of the model state vector, i.e. if the observation operator L H involves interpolation/extrapolation or some change of variables, then the nudging matrix K will contain the adjoint operations, i.e. some interpolation/extrapolation back to the model grid points, or the inverse change of variable.

As in the forward part of the algorithm, for simplicity reasons we make the following choice for the backward nudging matrix:

. The only parameters of the BFN algorithm are then the coefficients and . In the forward mode, is usually chosen such that the nudging term remains small in comparison with the other terms in the model equation. The coefficient is usually chosen to be the smallest coefficient that makes the numerical backward integration stable.

1 ' - = R H k K T k ' k 0 > k ' k

Experimental approach

The experimental approach consists of performing twin experiments with simulated data. First, a reference experiment is run and the corresponding data are extracted. From now on this reference trajectory will be called the exact solution. Experimental data are supposed to be obtained every grid-points of the model, and every time steps. The simulated data are then optionally noised with a Gaussian white noise distribution, and provided as observations to the assimilation scheme. The first guess of the assimilation experiments is chosen to be either a constant field or the reference model state some time before the beginning of the assimilation period. Finally, the results of the assimilation process are compared with the exact solution.

Shallow water model

The shallow water model (or Saint-Venant's equations) is a basic model, representing quite well the temporal evolution of geophysical flows. This model is usually considered for simple numerical experiments in oceanography, meteorology or hydrology. The shallow water equations are a set of three equations, describing the evolution of a two-dimensional horizontal flow. These equations are derived from a vertical integration of the three-dimensional fields, assuming the hydrostatic approximation, i.e. neglecting the vertical acceleration. There are several ways to write the shallow water equations, considering either the geopotential or height or pressure variables. We consider here the following configuration:

⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ = ∂ + ∂ + ∂ Δ + - = ∂ + + + ∂ Δ + - = ∂ + + - ∂ , 0 ) ( ) ( , ) ( , ) ( 0 0 hv hu h v rv h B u f v u ru h B v f u y x t y t x t ν ρ τ ς ν ρ τ ς (6) Paper.273.4 where u v y x ∂ - ∂ = ς is the relative vorticity, ) ( 2 1 2 2 * v u h g B + + =
is the Bernoulli potential, is the reduced gravity, is the Coriolis parameter (in the * g f β -plane approximation), 0 ρ is the water density, r is the friction coefficient, and ν is the viscosity (or dissipation) coefficient. The unknowns are and the horizontal components of the velocity, and the geopotential height. Finally, u v h τ is the forcing term of the model (e.g. the wind stress) [START_REF] Durbiano | Vecteurs caractéristiques de modèles océaniques pour la réduction d'ordre en assimilation de données[END_REF].

In the following experiments (see [START_REF] Auroux | The back and forth nudging algorithm applied to a shallow water model, comparison and hybridization with the 4D-VAR[END_REF] for all numerical results), the observations are available every 5 gridpoints and 24 time steps. The assimilation period lasts 720 time steps, corresponding to 15 days, and the forecast period ends 2880 time steps (or 2 months) after the initial time. As shown in Table 1, the BFN converges in 5 iterations, and we compared it with the 4D-VAR, either at convergence (16 iterations), or stopped after a similar time computation (5 iterations). Several conclusions can be drawn from this table. First the BFN is more efficient than the 4D-VAR in the same time (5 iterations of each). Then, if we compare the results when both algorithms have achieved convergence, the initial condition identified by the 4D-VAR is closer to the true initial state, but the BFN algorithm provides a better solution between time T and 4T, which is a keypoint for the quality of predictions.

Table 1 -Relative error of the forecast solutions corresponding to the BFN (5 iterations, converged), 4D-VAR (16 iterations, converged) and 4D-VAR (5 iterations) identified initial conditions, for the three variables, at various times: initial time, end of the assimilation period, and end of the prediction period.

Figure 1 shows a comparison between the BFN and 4D-VAR algorithms, and also a hybrid scheme. As the BFN provides at each iteration a new estimation of the initial condition, it is very easy to hybridize it with another data assimilation method, e.g. the 4D-VAR. We compare in Figure 1 the evolution of the identified trajectories by the BFN and 4D-VAR algorithms (stopped after 5 iterations only), and also the BFN+4D-VAR hybrid scheme in the same computation time. We compare then 5 iterations of BFN, 5 iterations of 4D-VAR, and 2 iterations of BFN followed by 3 iterations of 4D-VAR. We can see that during the assimilation period, the best solution is provided by the hybrid algorithm. Note that it would have been improved a little bit by the 4D-VAR at convergence, but we stop all algorithms after 5 iterations here. But during the forecast period, it is still provided by the BFN algorithm, except at the beginning. The results are comparable for u and v, except that the BFN gives a better solution a little earlier. In all the experiments performed, with several levels of noise and different background states, the hybrid method was the best during the assimilation period, and sometimes it was also the best at the beginning of the prediction period, but in most cases, the BFN solution is the best at the end. We can deduce from these results that a very cheap and simple way to largely improve the 4D-VAR algorithm is to perform very few (2 or 3) iterations of BFN before starting the 4D-VAR minimization process.

Quasi-geostrophic ocean model

We have also considered a layered quasi-geostrophic ocean model. This model arises from the primitive equations (conservation laws of mass, momentum, temperature and salinity), assuming first that the rotational effect (Coriolis force) is much stronger than the inertial effects. The Rossby number, ratio between the characteristic time of the Earth's rotation and the inertial time, must then be small compared to 1. Second, the thermodynamic effects are completely neglected in this model. Quasigeostrophy assumes that the horizontal dimension of the ocean is small compared to the size of the Earth, with a ratio of the order of the Rossby number. We finally assume that the depth of the basin is small compared to its width. In the case of the Atlantic Ocean, not all these assumptions are valid, notably the horizontal extension of the ocean. But it has been shown that the quasi-geostrophic approximation is fairly robust in practice, and that this approximate model reproduces quite well the ocean circulations at mid-latitudes, such as the jet stream (e.g. Gulf Stream in the case of the North Atlantic Ocean) and ocean boundary currents. The model system is then composed of coupled equations resulting from the conservation law of the potential vorticity. The equations can be written as: Finally, is the forcing term of the model, the wind stress applied to the ocean surface. We refer to [START_REF] Blayo | Assimilation of TOPEX/POSEIDON altimeter data into a circulation model of the North Atlantic[END_REF] for more details about this model and its equations.

n [. , 0 ] ), ( ) ( ) ) ( ( 1 1 6 4 1 T k F A k A Dt f D k n n k k × Ω = Ψ ∇ + ΔΨ + + Ψ δ δ θ (

F

The numerical experiments have been made on a three-layered square ocean. The basin has horizontal dimensions of 4000 km x 4000 km and its depth is 5 km. The layers' depths are 300m for the surface layer, 700m for the intermediate layer, and 4000m for the bottom layer. The ocean is discretized by a Cartesian mesh of 200×200×3 grid points. The time step is 1.5 h. The initial conditions are chosen equal to zero for a six-year ocean spin-up phase, the final state of which becoming the initial state of the data assimilation period. Then the assimilation period starts (time t=0) with this initial condition, and lasts 5 days (time t=T), i.e. 80 time steps. Data are available every 5 gridpoints of the model, with a time sampling of 7.5 h (i.e. every 5 time steps). We refer to [START_REF] Auroux | A nudging-based data assimilation method for oceanographic problems: the Back and Forth Nudging (BFN) algorithm[END_REF] for all these numerical results. Figure 2 shows, for each of the three layers, the RMS relative difference between the BFN trajectory (computed using the last BFN state at time t=0 as initial condition) and the reference trajectory as a solid line, and between the 4D-VAR trajectory (computed using the last initial state produced by the minimization process) and the reference trajectory as a dotted line, versus time. The first 5 days correspond to the assimilation period, and the next 15 days correspond to the forecast period.

The first point is that the BFN and 4D-VAR reconstruction errors have similar global behaviours: first decreasing at the beginning, during the assimilation period, and then increasing over the forecast period (see e.g. [START_REF] Pires | On extending the limits of variational assimilation in nonlinear chaotic systems[END_REF] for a detailed study of the 4D-VAR estimation error). Even if the reconstruction error on the initial condition is much higher with the BFN algorithm than with 4D-VAR, the BFN error decreases much more steeply and for a longer time than the 4D-VAR one, and increases less quickly at the end of the forecast period. The quality of the initial condition reconstruction is better using the 4D-VAR algorithm, but the BFN algorithm provides a comparable final estimation, and even a better forecast.

Another interesting point is that, even if only surface observations are assimilated, the identification of the intermediate and bottom layers is quite good. The 4D-VAR algorithm was already known to propagate surface information to all layers [START_REF] Luong | A variational method for the resolution of a data assimilation problem in oceanography[END_REF], but it is also the case for the BFN algorithm, even if this algorithm is less efficient on the bottom layer.

Conclusions

We have proved on several idealized situations that, provided that the feedback term is large enough as well as the assimilation period, we have convergence of the algorithm to the real initial state.

Numerical results have been performed on Burgers, shallow-water and layered quasi-geostrophic models. Numerical tests are currently investigated on a full primitive ocean model (NEMO). Twice less iterations than the 4D-VAR are necessary to obtain the same level of convergence. This algorithm is hence very promising in order to obtain a correct trajectory, with a smaller number of iterations than in a variational method, with a very easy implementation.

Fig. 1 -

 1 Fig. 1 -Relative difference between the true solution and the forecast trajectory corresponding to the BFN, 4D-VAR and BFN-preprocessed 4D-VAR identified initial conditions, after 5 iterations, versus time, for the height variable in the case of noisy observations.

  7) Ω is the circulation basin, is the stream function at layer , k Ψ k k θ is the sum of the dynamical and thermal vorticities at layer , is the Coriolis force, and the dissipative terms correspond to the lateral friction and the bottom friction dissipation.

Fig. 2 -

 2 Fig. 2 -Evolution in time of the RMS relative difference between the reference trajectory and the identified trajectories for the BFN (solid line) and 4D-VAR (dotted line) algorithms, versus time, for each layer: from surface (1) to bottom (3).