Back to the future : the Back and Forth Nudging algorithm - - PowerPoint PPT Presentation

Didier Auroux Mathematics Institute of Toulouse, France Near future : University of Nice auroux@math.univ-toulouse.fr auroux@unice.fr Work in collaboration with : J. Blum (U. Nice), M. Nodet (U. Grenoble)

Back to the future : the “Back and Forth Nudging” algorithm

Conference on Applied Inverse Problems, Vienna, Austria. July 20-24, 2009

Motivations

Environmental and geophysical studies : forecast the natural evolution retrieve at best the current state (or initial condition) of the environment. Geophysical fluids (atmosphere, oceans, . . .) : turbulent systems = ⇒ high sensitivity to the initial condition = ⇒ need for a precise identification (much more than observations) Environmental problems (ground pollution, air pollution, hurricanes, . . .) : problems of huge dimension, generally poorly modelized or observed Data assimilation consists in combining in an optimal way the observations of a system and the knowledge of the physical laws which govern it. Main goal : identify the initial condition, or estimate some unknown parame- ters, and obtain reliable forecasts of the system evolution.

AIP 2009, Vienna, July 23 2009 1/28

Data assimilation

t Observations Model combination model + observations ⇓ identification of the initial condition

• f a geophysical system

Fundamental for a chaotic system (atmosphere, ocean, . . .) Issue : These systems are generally irreversible Goal : Combine models and data Typical inverse problem : retrieve the system state from sparse and noisy

• bservations

AIP 2009, Vienna, July 23 2009 2/28

4D-VAR

Observations t Xb X0 Model governed by a system of ODE :      dX dt = F(X), 0 < t < T X(0) = X0 Xobs(t) : observations of the system, H : observation operator, Xb : background state (a priori estimation of the initial condition), B and R : covariance matrices of background and observation errors. J(X0) = 1 2(X0 − Xb)T B−1(X0 − Xb) + 1 2 T [Xobs(t) − H(X(t))]T R−1 [Xobs(t) − H(X(t))] dt

AIP 2009, Vienna, July 23 2009 3/28

Optimality system

Optimization under constraints : L (X0, X, P) = J(X0) + T

• P, dX

dt − F(X)

• dt

Direct model :    dX dt = F(X) X(0) = X0 Adjoint model :      −dP dt = ∂F ∂X T P + HT R−1 [Xobs(t) − H(X(t))] P(T) = 0 Gradient of the cost-function : ∂J ∂X0 = B−1(X0 − Xb) − P(0) Optimal solution : X0 = Xb + BP(0)

[Le Dimet-Talagrand (86)] AIP 2009, Vienna, July 23 2009 4/28

Forward nudging

     dX dt = F(X)+K(Xobs − H(X)), 0 < t < T, X(0) = X0, where K is the nudging (or gain) matrix. In the linear case (where F is a matrix), the forward nudging is called Luenberger or asymptotic observer. – Meteorology : Hoke-Anthes (1976) – Oceanography (QG model) : Verron-Holland (1989) – Atmosphere (meso-scale) : Stauffer-Seaman (1990) – Optimal determination of the nudging coeffcients : Zou-Navon-Le Dimet (1992), Stauffer-Bao (1993), Vidard-Le Dimet-Piacentini (2003)

AIP 2009, Vienna, July 23 2009 5/28

Forward nudging : linear case

Luenberger observer, or asymptotic observer (Luenberger, 1966)        dX dt = FX+K(Xobs − HX), d ˆ X dt = F ˆ X, Xobs = H ˆ X. d dt(X − ˆ X) = (F−KH)(X − ˆ X) If F − KH is a Hurwitz matrix, i.e. its spectrum is strictly included in the half-plane {λ ∈ C; Re(λ) < 0}, then X → ˆ X when t → +∞.

AIP 2009, Vienna, July 23 2009 6/28

Backward nudging

How to recover the initial state from the final solution ? Backward model :      d ˜ X dt = F( ˜ X), T > t > 0, ˜ X(T) = ˜ XT . If we apply nudging to this backward model :      d ˜ X dt = F( ˜ X)−K(Xobs − H ˜ X), T > t > 0, ˜ X(T) = ˜ XT .

AIP 2009, Vienna, July 23 2009 7/28

BFN : Back and Forth Nudging algorithm

Iterative algorithm (forward and backward resolutions) : ˜ X0(0) = Xb (first guess)      dXk dt = F(Xk)+K(Xobs − H(Xk)) Xk(0) = ˜ Xk−1(0)      d ˜ Xk dt = F( ˜ Xk)−K′(Xobs − H( ˜ Xk)) ˜ Xk(T) = Xk(T) If Xk and ˜ Xk converge towards the same limit X, and if K = K′, then X satisfies the state equation and fits to the observations.

AIP 2009, Vienna, July 23 2009 8/28

Choice of the direct nudging matrix K

Implicit discretization of the direct model equation with nudging : Xn+1 − Xn ∆t = FXn+1 + K(Xobs − HXn+1). Variational interpretation : direct nudging is a compromise between the mini- mization of the energy of the system and the quadratic distance to the obser- vations : min

X

1 2X − Xn, X − Xn − ∆t 2 FX, X + ∆t 2 R−1(Xobs − HX), Xobs − HX

• ,

by chosing K = HT R−1 where R is the covariance matrix of the errors of observation.

[Auroux-Blum (08)] AIP 2009, Vienna, July 23 2009 9/28

Choice of the backward nudging matrix K′

The feedback term has a double role :

• stabilization of the backward resolution of the model (irreversible system)
• feedback to the observations

If the system is observable, i.e. rank[H, HF, . . . , HF N−1] = N, then there exists a matrix K′ such that −F − K′H is a Hurwitz matrix (pole assignment method). Simpler solution : one can define K′ = k′HT , where k′ is e.g. the smallest value making the backward numerical integration stable.

AIP 2009, Vienna, July 23 2009 10/28

Example of convergence results

Viscous linear transport equation :    ∂tu − ν∂xxu + a(x)∂xu = −K(u − uobs), u(x, t = 0) = u0(x) ∂t˜ u − ν∂xx˜ u + a(x)∂x˜ u = K′(˜ u − uobs), ˜ u(x, t = T) = uT (x) We set w(t) = u(t) − uobs(t) and ˜ w(t) = ˜ u(t) − uobs(t) the errors.

• If K and K′ are constant, then ∀t ∈ [0, T] :

w(t) = e(−K−K′)(T −t)w(t)

(still true if the observation period does not cover [0, T])

• If the domain is not fully observed, then the problem is ill-posed.

Error after k iterations : wk(0) = e−[(K+K′)kT ]w0(0) exponential decrease of the error, thanks to :

• K + K′ : infinite feedback to the observations (not physical)
• T : asymptotic observer (Luenberger)
• k : infinite number of iterations (BFN)

AIP 2009, Vienna, July 23 2009 11/28

Shallow water model

∂tu − (f + ζ)v + ∂xB = τx ρ0h − ru + ν∆u ∂tv + (f + ζ)u + ∂yB = τy ρ0h − rv + ν∆v ∂th + ∂x(hu) + ∂y(hv) =

• ζ = ∂xv − ∂yu is the relative vorticity ;
• B = g∗h + 1

2(u2 + v2) is the Bernoulli potential ;

• g∗ = 0.02 m.s−2 is the reduced gravity ;
• f = f0 + βy is the Coriolis parameter (in the β-plane approximation), with f0 =

7.10−5 s−1 and β = 2.10−11 m−1.s−1 ;

• τ = (τx, τy) is the forcing term of the model (e.g. the wind stress), with a maximum

amplitude of τ0 = 0.05 s−2 ;

• ρ0 = 103 kg.m−3 is the water density ;
• r = 9.10−8 s−1 is the friction coefficient.
• ν = 5 m2.s−1 is the viscosity (or dissipation) coefficient.

AIP 2009, Vienna, July 23 2009 12/28

Shallow water model

2D shallow water model, state = height h and horizontal velocity (u, v) Numerical parameters : (run example) Domain : L = 2000 km × 2000 km ; Rigid boundary and no-slip BC ; Time step = 1800 s ; Assimilation period : 15 days ; Forecast period : 15 + 45 days Observations : of h only (∼ satellite obs), every 5 gridpoints in each space direction, every 24 hours. Background : true state one month before the beginning of the assimilation period + white gaussian noise (∼ 10%) Comparison BFN - 4DVAR : sea height h ; velocity :u and v.

AIP 2009, Vienna, July 23 2009 13/28

Convergence - perfect obs.

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 100 200 300 400 500 600 700 RMS error Time steps h 0.05 0.1 0.15 0.2 0.25 100 200 300 400 500 600 700 RMS error Time steps u 0.05 0.1 0.15 0.2 0.25 0.3 100 200 300 400 500 600 700 RMS error Time steps v

Relative difference between the BFN iterates (5 first iterations) and the true solution versus the time steps, for h, u and v.

AIP 2009, Vienna, July 23 2009 14/28

Comparaison - perfect obs.

Relative error h u v Background state 37.6% 21.5% 30.3% BFN (5 iterations, converged) 0.44% 1.78% 2.41% 4D-VAR (5 iterations) 0.64% 3.14% 4.47% 4D-VAR (18 iterations, converged) 0.61% 2.43% 3.46% Relative error of the background state and various identified initial conditions for the three variables.

AIP 2009, Vienna, July 23 2009 15/28

Noisy observations

0.05 0.1 0.15 0.2 500 1000 1500 2000 2500 3000 RMS error Time steps h u v 0.05 0.1 0.15 0.2 500 1000 1500 2000 2500 3000 RMS error Time steps h u v

Relative difference bet- ween the true solution and the forecast trajec- tory corresponding to the BFN (top) and 4D-VAR (bottom) identified initial conditions at convergence, versus time, for the three variables and in the case

• f noisy observations.

AIP 2009, Vienna, July 23 2009 16/28

BFN-preprocessed 4D-VAR

0.01 0.02 0.03 0.04 0.05 0.06 500 1000 1500 2000 2500 3000 RMS error Time steps BFN 4D-VAR BFN + 4D-VAR

Relative difference between the true solution and the forecast trajectory corres- ponding 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.

AIP 2009, Vienna, July 23 2009 17/28

Quasi-geostrophic ocean model

The oceanic model used in this study is based on the quasi-geostrophic ap- proximation obtained by writing the conservation of the potential vorticity. The vertical structure of the ocean is divided into N layers. Each one has a constant density ρk with a depth Hk (k = 1, . . . , N). We get a coupled system

• f N equations :

Dk(θk(Ψ) + f) Dt + δk,N C1∆ΨN − C3∆3Ψk = Fk in Ω×]0, T[, ∀k = 1, . . . , N; where :

• Ω ⊂ I

R2 is the oceanic basin, and [0, T] the time interval for the study ;

• Ψk is the stream function in the layer k ;
• C1∆ΨN is the dissipation on the bottom of the ocean ;
• C3∆3Ψk is the parametrization of internal and subgrid dissipation ;

AIP 2009, Vienna, July 23 2009 18/28

Example

• θk(Ψ) is the potential vorticity in the layer k, given by :

   θ1(Ψ) . . . θN(Ψ)    = [∆ − [W]]    Ψ1 . . . ΨN    where [W] is a N × N tridiagonal matrix ;

• f is the Coriolis force ;
• Dk.

Dt is the Lagrangian derivative in layer k, given by : Dk. Dt = ∂. ∂t − ∂Ψk ∂y ∂. ∂x + ∂Ψk ∂x ∂. ∂y = ∂. ∂t + J(Ψk, .), where J(., .) is the Jacobian operator J(ϕ, ξ) = ∂ϕ ∂x ∂ξ ∂y − ∂ϕ ∂y ∂ξ ∂x;

• Fk is a forcing term. In this model only the tension due to the wind, denoted

τ, is taken into account : F1 = Rotτ and Fk = 0, ∀k ≥ 2.

AIP 2009, Vienna, July 23 2009 19/28

Example

We consider altimetric measurement of the surface of the ocean given by satel- lite observations (Topex-Poseidon, Jason). The observed data is the change in the surface of the ocean. According to the quasi geostrophic approximation it is proportional to the stream function in the surface layer : hobs = f0 g Ψobs

1

Therefore we will assimilate surface data in order to retrieve the fluid circula- tion especially in the deep ocean layers. The control vector is the initial state on the N layers : u =

• Ψk(t = 0)
• k=1,...,N ∈ Uad

The state vector is

• Ψk(t)
• k=1,...,N

AIP 2009, Vienna, July 23 2009 20/28

QG - BFN convergence

0.02 0.04 0.06 0.08 0.1 0.12 0.14 5 10 15 20 RMS relative difference BFN iterations 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 5 10 15 20 RMS relative error BFN iterations

RMS relative difference between two consecutive BFN iterates (a) and bet- ween the BFN iterates and the exact solution (b) versus the number of BFN iterations.

AIP 2009, Vienna, July 23 2009 21/28

QG - BFN vs 4D-VAR

10000 1e+06 1e+08 1e+10 1e+12 1e+14 5 10 15 20 BFN iterations J ||gradJ1|| ||gradJ2|| ||gradJ3|| 10000 1e+06 1e+08 1e+10 1e+12 1e+14 5 10 15 20 4D-VAR iterations J ||gradJ1|| ||gradJ2|| ||gradJ3||

Evolution of the 4D-VAR cost function and of the 3 gradients of the cost func- tion in the 3 ocean layers for the BFN iterates versus the number of BFN ite- rations (top) and for the 4D-VAR ite- rates versus the number of 4D-VAR ite- rations (b).

AIP 2009, Vienna, July 23 2009 22/28

QG - BFN convergence (noisy obs.)

0.02 0.04 0.06 0.08 0.1 0.12 0.14 5 10 15 20 RMS relative difference BFN iterations 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 5 10 15 20 RMS relative error BFN iterations

RMS relative difference between two consecutive BFN iterates (a) and bet- ween the BFN iterates and the exact solution (b) versus the number of BFN iterations in the case of noised observa- tions.

AIP 2009, Vienna, July 23 2009 23/28

QG - BFN vs 4D-VAR (noisy obs.)

10000 1e+06 1e+08 1e+10 1e+12 1e+14 5 10 15 20 BFN iterations J gradJ1 gradJ2 gradJ3 10000 1e+06 1e+08 1e+10 1e+12 1e+14 5 10 15 20 4D-VAR iterations J gradJ1 gradJ2 gradJ3

Evolution of the 4D-VAR cost function and of the 3 gradients of the cost func- tion in the 3 ocean layers for the BFN iterates versus the number of BFN ite- rations (top) and for the 4D-VAR ite- rates versus the number of 4D-VAR ite- rations (b), in the case of noisy obser- vations.

AIP 2009, Vienna, July 23 2009 24/28

QG - BFN vs 4D-VAR (noisy obs.)

0.01 0.02 0.03 0.04 0.05 0.06 5 10 15 20 RMS relative error BFN layer 1 4D-VAR layer 1 0.01 0.02 0.03 0.04 0.05 0.06 5 10 15 20 RMS relative error Time (days) BFN layer 3 4D-VAR layer 3 0.01 0.02 0.03 0.04 0.05 0.06 5 10 15 20 RMS relative error BFN layer 2 4D-VAR layer 2

Evolution in time of the RMS relative difference between the reference trajec- tory and the identified trajectories for the BFN (solid line) and 4D-VAR (dot- ted line) algorithms, versus time, for each layer : from surface to bottom.

AIP 2009, Vienna, July 23 2009 25/28

QG - BFN convergence (model error)

0.02 0.04 0.06 0.08 0.1 0.12 0.14 5 10 15 20 RMS relative difference BFN iterations 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 5 10 15 20 RMS relative error BFN iterations

RMS relative difference between two consecutive BFN iterates (a) and bet- ween the BFN iterates and the exact solution (b) versus the number of BFN iterations in the case of imperfect mo- del.

AIP 2009, Vienna, July 23 2009 26/28

QG - BFN vs 4D-VAR (model error)

10000 1e+06 1e+08 1e+10 1e+12 1e+14 5 10 15 20 BFN iterations J gradJ1 gradJ2 gradJ3 10000 1e+06 1e+08 1e+10 1e+12 1e+14 5 10 15 20 4D-VAR iterations J gradJ1 gradJ2 gradJ3

Evolution of the 4D-VAR cost function and of the 3 gradients of the cost func- tion in the 3 ocean layers for the BFN iterates versus the number of BFN ite- rations (a) and for the 4D-VAR iterates versus the number of 4D-VAR itera- tions (b) in the case of imperfect model.

AIP 2009, Vienna, July 23 2009 27/28

Conclusions

• Easy implementation (no linearization, no adjoint state, no minimization

process)

• Very efficient in the first iterations (faster convergence)
• Lower computational and memory costs than other DA methods
• Stabilization of the backward model
• Excellent preconditioner for 4D-VAR (or Kalman filters)

Under investigation :

• Identification of the entire solution of a full primitive model from the

knowledge of the SSH (or SST) only

• Efficient resolution of Riccati-like equations for a better backward nudging

matrix

AIP 2009, Vienna, July 23 2009 28/28

Observability condition

Let χ(x) be the time during which the characteristic curve with foot x lies in the support of K. Then the system is observable if and only if min

x χ(x) > 0.

Partial observations in space : half of the domain is observed.

0.2 0.4 0.6 0.8 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x T=0.5 T=0.3 T=0.15 T=0.05 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x T=1 T=0.5 T=0.3 T=0.15 T=0.05

Decrease rate of the error after one iteration of BFN as a function of the space variable x, for various final times T.

Linear case (left) : theoretical observability condition = T > 0.5 Nonlinear case (right) : numerical observability condition = T > 1

AIP 2009, Vienna, July 23 2009 29/28