An iterative ensemble Kalman smoother Marc Bocquet 1 , 2 Pavel Sakov - PowerPoint PPT Presentation

An iterative ensemble Kalman smoother Marc Bocquet 1 , 2 Pavel Sakov 3 1 Universit´ e Paris-Est, CEREA, joint lab ´ Ecole des Ponts ParisTech and EdF R&D, France 2 INRIA, Paris-Rocquencourt Research center, France 3 Bureau of Meteorology, Australia (bocquet@cerea.enpc.fr) 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 1 / 26

Reminders: the EnKF and the EnKF-N Reminder: failure of the raw ensemble Kalman filter (EnKF) ◮ EnKF relies for its analysis on the first and second-order empirical moments: N N x = 1 1 ∑ ∑ ( x k − x )( x k − x ) T . x k , P = (1) N N − 1 n =1 k =1 Lorenz ’95 N=20 ∆ t=0.05 ◮ With the exception of Gaussian and linear systems, the EnKF fails to pro- vide a proper estimation on most sys- EnKF without inflation RMSE analysis tems. 1 EnKF with inflation λ =1.02 EnKF-N ◮ To properly work, it needs clever but ad hoc fixes: localisation and in- 0.1 flation. 0 1000 2000 Analysis cycle ◮ In a perfect model context, the finite-size EnKF (EnKF-N) avoids tuning inflation. 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 2 / 26

Reminders: the EnKF and the EnKF-N Reminder: principle of the EnKF-N ◮ The prior of EnKF and the prior of EnKF-N: � � − 1 2 ( x − x ) T P − 1 ( x − x ) p ( x | x , P ) ∝ exp � � − N � � 2 , � ( x − x )( x − x ) T + ε N ( N − 1) P p ( x | x 1 , x 2 ,..., x N ) ∝ (2) � with ε N = 1 (mean-trusting variant), or ε N = 1+ 1 N (original variant). ◮ Ensemble space decomposition (ETKF version of the filters): x = x + Aw . ◮ The variational principle of the analysis (in ensemble space): J ( w ) = 1 2 ( y − H ( x + Aw )) T R − 1 ( y − H ( x + Aw ))+ N − 1 w T w 2 � � J ( w )= 1 2 ( y − H ( x + Aw )) T R − 1 ( y − H ( x + Aw ))+ N ε N + w T w 2 ln . (3) 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 3 / 26

Reminders: the EnKF and the EnKF-N Reminder: the EnKF-N algorithm Requires: The forecast ensemble { x n } n =1 ,..., N , the observations y , and error 1 covariance matrix R Compute the mean x and the anomalies A from { x k } k =1 ,..., N . 2 Compute Y = HA , δ = y − Hx 3 Find the minimum: 4 � � �� ( δ − Yw ) T R − 1 ( δ − Yw )+ N ln ε N + w T w w a = min w Compute x a = x + Aw a . 5 � � − 1 Y T R − 1 Y + N ( ε N + w T a w a ) I N − 2 w a w T a Compute Ω a = 6 ( ε N + w T a w a ) 2 Compute W a = { ( N − 1) Ω a } 1 / 2 U 7 k = x a + AW a Compute x a 8 k 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 4 / 26

Context Iterative Kalman filters: context ◮ The iterative extended Kalman filter [Wishner et al., 1969; Jazwinski, 1970] IEKF ◮ The iterative extended Kalman smoother [Bell, 1994] IEKS Much too costly + needs the TLM and the adjoint − → ensemble methods ◮ The iterative ensemble Kalman filter [Sakov et al., 2012; Bocquet and Sakov, 2012] IEnKF ◮ The iterative ensemble Kalman smoother [This talk. . . ] IEnKS It’s TLM and adjoint free! Don’t want to be bothered by inflation tuning? ◮ The finite-size iterative ensemble Kalman filter [Bocquet and Sakov, 2012] IEnKF-N ◮ The finite-size iterative ensemble Kalman smoother [This talk. . . ] IEnKS-N 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 5 / 26

Context Iterative ensemble Kalman smoother (IEnKS): context ◮ An extension of the iterative ensemble Kalman filter (IEnKF), a fairly recent idea: [Gu & Oliver, 2007] : Initial idea. [Kalnay & Yang, 2010-2012] : A closely related idea. [Sakov, Oliver & Bertino, 2012] : The “pi` ece de r´ esistance” [Bocquet & Sakov, 2012] : Bundle scheme + ensemble transform form. ◮ Related but not to be confused with the iterative ensemble Smoother (IEnS) in the oil reservoir modelling smoothers, where cycling is not an issue. ◮ Assumptions of the present study: Perfect model. In the rank-sufficient regime. Localisation is more challenging (but possible) in this context (Pavel’s talk). Looking for the best performance. Numerical cost secondary. 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 6 / 26

The interative ensemble Kalman smoother Theory Iterative ensemble Kalman smoother: the cycling ◮ L : length of the data assimilation window; S : shift of the data assimilation window in between two updates. y y L−3 L−2 t L−1 t t L−3 L−2 y L y L−1 S ∆ t t L t t t t L−1 0 1 L y y L+1 L+2 S ∆ t t L+1 t t L+1 L+2 L t ∆ ◮ This may or may not lead to overlapping windows. Here, we study the case S = 1, which is close to quasi-static conditions [Pires et al., 1996] . ◮ Let us first focus on the single data assimilation (SDA) scheme. 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 7 / 26

The interative ensemble Kalman smoother Theory SDA IEnKS: a variational standpoint ◮ Analysis IEnKS cost function in state space p ( x 0 | y L ) ∝ exp( − J ( x 0 )): J ( x 0 ) =1 2 ( y L − H L ◦ M L ← 0 ( x 0 ))) T R − 1 L ( y L − H L ◦ M L ← 0 ( x 0 ))) + 1 2 ( x 0 − x 0 ) P − 1 0 ( x 0 − x 0 ) . (4) ◮ Reduced scheme in ensemble space, x 0 = x 0 + A 0 w , where A 0 is the ensemble anomaly matrix: � J ( w ) = J ( x 0 + A 0 w ) . (5) ◮ IEnKS cost function in ensemble space: J ( w ) =1 2 ( y L − H L ◦ M L ← 0 ( x 0 + A 0 w )) T R − 1 � L ( y L − H L ◦ M L ← 0 ( x 0 + A 0 w )) + 1 2( N − 1) w T w . (6) 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 8 / 26

The interative ensemble Kalman smoother Theory SDA IEnKS: minimisation scheme ◮ As a variational reduced method, one can use Gauss-Newton [Sakov et al., 2012] , Levenberg-Marquardt [Bocquet and Sakov, 2012; Chen and Oliver, 2013] , etc, minimisation schemes (not limited to quasi-Newton). ◮ Gauss-Newton scheme: w ( p +1) = w ( p ) − � H − 1 ( p ) ∇ � J ( p ) ( w ( p ) ) , x ( p ) = x (0) + A 0 w ( p ) , 0 0 � � J ( p ) = − Y T y L − H L ◦ M L ← 0 ( x ( p ) +( N − 1) w ( p ) , ∇ � ( p ) R − 1 0 ) L H ( p ) = ( N − 1) I N + Y T � ( p ) R − 1 L Y ( p ) , Y ( p ) = [ H L ◦ M L ← 0 A 0 ] ′ ( p ) . (7) ◮ One alternative to compute the sensitivities: the bundle scheme. It simply mimics the action of the tangent linear by finite difference: �� � � I N − 11 T Y ( p ) ≈ 1 x ( p ) 1 T + ε A 0 ε H L ◦ M 1 ← 0 . (8) N 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 9 / 26

The interative ensemble Kalman smoother Theory IEnKS: ensemble update and the forecast step ◮ Generate an updated ensemble from the previous analysis: √ 0 1 T + H − 1 / 2 E ⋆ 0 = x ⋆ N − 1 A 0 � where U1 = 1 . (9) U ⋆ ◮ Just propagate the updated ensemble from t 0 to t S : E S = M S ← 0 ( E 0 ) . (10) In the quasi-static case: S = 1. 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 10 / 26

The interative ensemble Kalman smoother Theory IEnKS: introducing the MDA scheme ◮ Suppose we could assimilate the observation vectors several times. . . β β β y y y L−1 L 0 L−3 L−2 −2 t L−1 t t L−3 L−2 β y L β y y L−1 β 0 L−1 L S ∆ t 0 t L t t t t L−1 0 1 L y β y β y L−1 L β 0 L+1 L+2 S ∆ t 2 t L+1 t t L+1 L+2 L t ∆ ◮ This leads to overlapping windows. Here, we study the quasi-static case S = 1. ◮ This is called multiple data assimilation (MDA) scheme. 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 11 / 26

The interative ensemble Kalman smoother Theory IEnKS: the MDA approach ◮ Two flavours of Multiple Data Assimilation: The splitting of observations: Following the partition 1 = ∑ L k =1 β k , the observation vector y with prior error covariance matrix is split into L partial observation y β k , with prior error covariance matrix β − 1 k R . It is a consistent approach in the Gaussian/linear limit, and one hopes it is still so in nonlinear conditions. The multiple assimilation of each observation with its original weights. It is correct but the filtering/smoothing pdf (essentially) becomes a power of the searched pdf! ◮ An extra step in the analysis. MDA IEnKS does not approximate per se the filtering pdf, but a more complex pdf. To approach the correct filtering/smoothing pdf, one needs an extra step, that we called the balancing step which re-weights the observations within the data assimilation window, and perform a final analysis. 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 12 / 26

The interative ensemble Kalman smoother Numerical applications The Lorenz ’95 model ◮ The toy-model [Lorenz and Emmanuel 1998] : It represents a mid-latitude zonal circle of the global atmosphere. M = 40 variables { x m } m =1 ,..., M . For m = 1 ,..., M : d x m = ( x m +1 − x m − 2 ) x m − 1 − x m + F , d t where F = 8, and the boundary is cyclic. Chaotic dynamics, topological dimension of 13, a doubling time of about 0 . 42 time units, and a Kaplan-Yorke dimension of about 27 . 1. ◮ Setup of the experiment: Time-lag between update: ∆ t = 0 . 05 (about 6 hours for a global model), fully observed, R = I . 12.5 35 10.0 30 7.5 25 5.0 20 2.5 15 0.0 10 -2.5 5 -5.0 -7.5 0 0 100 200 300 400 500 8 th EnKF workshop, Bergen, Norway, 27-29 May 2013 M. Bocquet 13 / 26