Using Ensemble Kalman Filters to Assimilate Dynamic Ocean Topography Data into a Global Ocean Model

(1)

Using Ensemble Kalman Filters

to Assimilate Dynamic Ocean Topography Data into a Global Ocean Model

Lars Nerger

Alfred Wegener Institute for Polar and Marine Research Bremerhaven, Germany

and

Bremen Supercomputing Competence Center BremHLR Bremen, Germany

Lars.Nerger@awi.de

IGG, Uni. Bonn, 20.6.2013

(2)

Lars Nerger – Assimilating DOT data with EnKFs

Outline

  Ensemble-based Kalman filters

  Implementation aspects

  Application with global ocean model

Collaborations:

AWI: W. Hiller, J. Schröter, A. Alexandrov, P. Kirchgessner, S. Loza, T. Janjic (now DWD)

BSH: F. Janssen, S. Massmann O.A.Sys GmbH: Reiner Schnur

(3)

Application Example

Information: Model Information: Observation

Model surface temperature Satellite surface temperature

•  Generally correct, but has errors

•  all fields, fluxes, …

•  Generally correct, but has errors

•  sparse information

(only surface, data gaps, one field) Combine both sources of information

quantitatively by computer algorithm

➜  data assimilation

(4)

Data Assimilation

  Optimal estimation of system state:

• 

initial conditions (for weather/ocean forecasts, …)

• 

state trajectory (temperature, concentrations, …)

•  parameters (growth of phytoplankton, …)

•  fluxes (heat, primary production, …)

•  boundary conditions and ‘forcing’ (wind stress, …)

€

  Characteristics of system:

• 

high-dimensional numerical model -

O

⁽¹⁰⁶^-10⁹⁾

•  sparse observations

•  non-linear

(5)

Ensemble-based Kalman Filters

(6)

Ensemble-based Kalman Filter

First formulated by G. Evensen (EnKF, 1994)

Kalman filter: express probability distributions by mean and covariance matrix

EnKF: Use ensembles to represent probability distributions

observation

time 0 time 1 time 2

analysis ensemble

forecast

ensemble transformation initial

sampling state

estimate

forecast

Looks trivial!

BUT:

There are many possible

choices!

(7)

  Properties and differences are not fully understood

ETKF Which filter should one use?

Many choices - a little “ zoo ” (not complete):

EAKF EnKF(94/98)

SEIK

EnSRF SEEK

RRSQRT ROEK

MLEF EnKF(2003)

EnKF(2004)

SPKF ESSE EnKF(94/98)

SEEK

SEIK

Studied in Nerger

et al. Tellus (2005)

RHF

(8)

Data Assimilation – Model and Observations

Two components:

1. State:

Dynamical model

€

x 2 R

ⁿ

x

_i

= M

_i _1,i

[x

_i ₁

]

2. Obervations:

Observation equation (relation of observation to state x):

Observation error covariance matrix:

y 2 R

^m

y = H [x]

R

(9)

a

The Ensemble Kalman Filter (EnKF, Evensen 94)

Ensemble

Analysis step:

Update each ensemble member

Kalman filter

5 EnKF

Init

x

^a₀

⌅ R

ⁿ

, P

^a₀

⌅ R

ⁿ^⇥ⁿ

(41) { x

^a(l)₀

, l = 1, . . . , N } (42) x

^a₀

= 1

N

⇧

N

l=1

x

^a(l)₀

⇥ x

^t₀

⇥

(43)

P ˜

^a₀

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)₀

x

^a₀

⌅⇤

x

^a(l)₀

x

^a₀

⌅

T

⇥ P

^a₀

(44)

P

^a₀

= LL

^T

, L ⌅ R

ⁿ^⇥^q

(45) x

^a(i)₀

= x

^a₀

+ Lb

⁽ⁱ⁾

, b

⁽ⁱ⁾

⌅ R

^q

(46)

⇤ N (0, 1) (47)

Forecast

x

^a(l)_i

= M

_i,i ₁

[x

^a(l)_i ₁

] +

^(l)_i

(48)

Analysis

{ y

^o(l)_k

, l = 1, . . . , N } (49) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

_k^o(l)

H

_k

⌃

x

^f_k^(l)

⌥⌅

(50) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

_k^o(l)

H

_k

x

^f_k^(l)

⌅

(51) K ˜

_k

= P ˜

^f_k

H

^T_k

⇤

H

_k

P ˜

^f_k

H

^T_k

+ R

_k

⌅

1

(52) K

_k

= P

^f_k

H

^T_k

⇤

H

_k

P

^f_k

H

^T_k

+ R

_k

⌅

1

(53) H

_k

P

^f_k

H

^T_k

+ R

_k

⌅ R

^m^⇥^m

(54) P ˜

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(55)

x

^a_k

:= 1 N

⇧

N

l=1

x

^a(l)_k

(56)

P ˜

^a_k

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)_k

x

^a_k

⌅⇤

x

^a(l)_k

x

^a_k

⌅

T

(57)

5 5 EnKF

Init

x

^a₀

⌅ R

ⁿ

, P

^a₀

⌅ R

ⁿ^⇥ⁿ

(41) { x

^a(l)₀

, l = 1, . . . , N } (42) x

^a₀

= 1

N

⇧

N

l=1

x

^a(l)₀

⇥ x

^t₀

⇥

(43)

P ˜

^a₀

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)₀

x

^a₀

⌅⇤

x

^a(l)₀

x

^a₀

⌅

T

⇥ P

^a₀

(44)

P

^a₀

= LL

^T

, L ⌅ R

ⁿ^⇥^q

(45) x

^a(i)₀

= x

^a₀

+ Lb

⁽ⁱ⁾

, b

⁽ⁱ⁾

⌅ R

^q

(46)

⇤ N (0, 1) (47)

Forecast

x

^a(l)_i

= M

_i,i ₁

[x

^a(l)_i ₁

] +

^(l)_i

(48)

Analysis

{ y

_k^o(l)

, l = 1, . . . , N } (49) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

_k^o(l)

H

_k

⌃

x

^f_k^(l)

⌥⌅

(50) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

^o(l)_k

H

_k

x

^f_k^(l)

⌅

(51) x

^a(l)_k

= x

^f_k^(l)

+ K

_k

⇤

y

_k^(l)

H

_k

x

^f_k^(l)

⌅

(52) K ˜

_k

= P ˜

^f_k

H

^T_k

⇤

H

_k

P ˜

^f_k

H

^T_k

+ R

_k

⌅

1

(53) K

_k

= P

^f_k

H

^T_k

⇤

H

_k

P

^f_k

H

^T_k

+ R

_k

⌅

1

(54) H

_k

P

^f_k

H

^T_k

+ R

_k

⌅ R

^m^⇥^m

(55) P ˜

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(56)

x

^a_k

:= 1 N

⇧

N

l=1

x

^a(l)_k

(57)

P ˜

^a_k

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)_k

x

^a_k

⌅⇤

x

^a(l)_k

x

^a_k

⌅

T

(58)

5 5 EnKF

Init

x

^a₀

⌅ R

ⁿ

, P

^a₀

⌅ R

ⁿ^⇥ⁿ

(41) { x

^a(l)₀

, l = 1, . . . , N } (42) x

^a₀

= 1

N

⇧

N

l=1

x

^a(l)₀

⇥ x

^t₀

⇥

(43)

P ˜

^a₀

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)₀

x

^a₀

⌅⇤

x

^a(l)₀

x

^a₀

⌅

T

⇥ P

^a₀

(44)

P

^a₀

= LL

^T

, L ⌅ R

ⁿ^⇥^q

(45) x

^a(i)₀

= x

^a₀

+ Lb

⁽ⁱ⁾

, b

⁽ⁱ⁾

⌅ R

^q

(46)

⇤ N (0, 1) (47)

Forecast

x

^a(l)_i

= M

_i,i ₁

[x

^a(l)_i ₁

] +

^(l)_i

(48)

Analysis

{ y

^o(l)_k

, l = 1, . . . , N } (49) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

^o(l)_k

H

_k

⌃

x

^f_k^(l)

⌥⌅

(50) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

_k^o(l)

H

_k

x

^f_k^(l)

⌅

(51) x

^a(l)_k

= x

^f_k^(l)

+ K

_k

⇤

y

_k^(l)

H

_k

x

^f_k^(l)

⌅

(52) K ˜

_k

= P ˜

^f_k

H

^T_k

⇤

H

_k

P ˜

^f_k

H

^T_k

+ R

_k

⌅

1

(53) K

_k

= P

^f_k

H

^T_k

⇤

H

_k

P

^f_k

H

^T_k

+ R

_k

⌅

1

(54) H

_k

P

^f_k

H

^T_k

+ R

_k

⌅ R

^m^⇥^m

(55) P ˜

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(56)

P

^f_k

:= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(57)

x

^a_k

:= 1 N

⇧

N

l=1

x

^a(l)_k

(58)

P ˜

^a_k

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)_k

x

^a_k

⌅⇤

x

^a(l)_k

x

^a_k

⌅

T

(59)

5 5 EnKF

Init

x

^a₀

⌅ R

ⁿ

, P

^a₀

⌅ R

ⁿ^⇥ⁿ

(41) { x

^a(l)₀

, l = 1, . . . , N } (42) x

^a₀

= 1

N

⇧

N

l=1

x

^a(l)₀

⇥ x

^t₀

⇥

(43)

P ˜

^a₀

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)₀

x

^a₀

⌅⇤

x

^a(l)₀

x

^a₀

⌅

T

⇥ P

^a₀

(44)

P

^a₀

= LL

^T

, L ⌅ R

ⁿ^⇥^q

(45) x

^a(i)₀

= x

^a₀

+ Lb

⁽ⁱ⁾

, b

⁽ⁱ⁾

⌅ R

^q

(46)

⇤ N (0, 1) (47)

Forecast

x

^a(l)_i

= M

_i,i ₁

[x

^a(l)_i ₁

] +

^(l)_i

(48)

Analysis

{ y

^o(l)_k

, l = 1, . . . , N } (49) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

^o(l)_k

H

_k

⌃

x

^f_k^(l)

⌥⌅

(50) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

_k^o(l)

H

_k

x

^f_k^(l)

⌅

(51) x

^a(l)_k

= x

^f_k^(l)

+ K

_k

⇤

y

_k^(l)

H

_k

x

^f_k^(l)

⌅

(52) K ˜

_k

= P ˜

^f_k

H

^T_k

⇤

H

_k

P ˜

^f_k

H

^T_k

+ R

_k

⌅

1

(53) K

_k

= P

^f_k

H

^T_k

⇤

H

_k

P

^f_k

H

^T_k

+ R

_k

⌅

1

(54) H

_k

P

^f_k

H

^T_k

+ R

_k

⌅ R

^m^⇥^m

(55) P ˜

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(56)

P

^f_k

:= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(57)

x

^a_k

:= 1 N

⇧

N

l=1

x

^a(l)_k

(58)

P ˜

^a_k

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)_k

x

^a_k

⌅⇤

x

^a(l)_k

x

^a_k

⌅

T

(59)

5

Ensemble

covariance matrix

Ensemble mean (state estimate)

(10)

A simple test problem

  Twin experiment with nonlinear shallow water model

  Initial state estimate: temporal mean state

  Initial covariance matrix: variability around mean state

(11)

Shallow water model: filter performances

  SEEK stagnates

  same convergence behavior for EnKF and SEIK

  smaller performance for EnKF than for SEIK

  EnKF ensemble 1.5-2 times larger than SEIK ensemble for same filtering performance

Error reduction due to assimilation

Ensemble size

Nerger, Hiller, Schröter. Tellus 57A (2005) 715-735

(12)

Some results: EnKF vs. SEIK

•  EnKF94/98

•  very simple to implement

•  costly (compute analysis update in observation space)

•  observation ensemble introduces sampling errors

•  random ensemble initialization has slow convergence

•  SEIK

•  more difficult to implement

•  much faster (analysis update in ensemble space)

•  faster convergence with initialization using singular value decomposition (empirical orthogonal functions)

€

Nerger, Hiller, Schröter. Tellus 57A (2005) 715-735

What makes SEIK faster than EnKF?

(13)

Two features of the SEIK filter

1.  Avoid perturbing observations

•  Apply two step update:

1.  Update ensemble mean state

2.  Transform forecast ensemble to represent analysis P

€

2.  Ensemble transformation in ensemble space

•  Degrees of freedom of analysis: ensemble size – 1

•  EnKF uses update in observation space (usually much larger than ensemble size)

Typical for ensemble square-root Kalman filters

(14)

Efficient use of ensembles

€

Kalman gain

K ˜

_k

= ˜ P

^f_k

H

^T_k

⇣

H

_k

P ˜

^f_k

H

^T_k

+ R

_k

⌘

1

K ˜

_k

= ⇣

P ˜

^f_k

⌘

1

+ H

^T

R

¹

H

1

H

^T

R

¹

Alternative form (Sherman-Morrison-Woodbury matrix identity)

Looks worse: matrices need inversion

n ⇥ n

K ˜

_k

= X

⁰

h

(N 1)I + X

⁰^T

H

^T

R

¹

HX

⁰

i

1

X

⁰^T

H

^T

R

¹

However: with ensemble

Inversion of matrix

(Ensemble perturbation matrix )

P ˜

^f_k

= (N 1)

¹

X

⁰

X

⁰^T

N ⇥ N

X

⁰

= X X ¯

(15)

ETKF Which filter should one use?

Many choices - a little “ zoo ” (not complete):

EAKF

ETKF EnKF(94/98)

SEIK

EnSRF SEEK

RRSQRT ROEK

EnKF(2003) EnKF(2004)

ESTKF EnKF(94/98)

SEEK

SEIK

Studied in Nerger et al. Tellus (2005)

SEIK

New study:

Nerger et al., Mon. Wea. Rev.

(2012)

New filter formulation

MLEF

SPKF

ESSE

RHF

(16)

Square root of covariance matrix (ensemble size N, state dim n)

T is specific for filter algorithm:

ETKF:

T removes ensemble mean

(usually, compute directly ) Z has dimension nN

SEIK:

T removes ensemble mean and drops last column Z has dimension n(N-1)

Analysis

X

^f_k

= ⌃

x

^f_k⁽¹⁾

, . . . , x

^f_k^(N⁾

⌥

(167) X

^f_k

= ⌃

x

^f_k

, . . . , x

^f_k

⌥

(168)

Z

^f_k

= X

^f_k

X

^f_k

(169)

Z = X

^f

T (170)

P ˇ

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(171) P ˇ

^f_k

= 1

N 1 Z

^f_k

⇤

Z

^f_k

⌅

T

(172) P ˇ

^f_k

= Z

^f_k

G ⇤

Z

^f_k

⌅

T

(173) G := 1

N 1 I (174)

x

^a_k

= x

^f_k

+ Z

^f_k

w

_k

(175) w

_k

= A

_k

(H

_k

Z

^f_k

)

^T

R

_k ¹

⇤

y

^o_k

H

_k

x

^f_k

⌅

(176)

A

_k ¹

= G

¹

+ (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(177) A

_k ¹

= (N 1)I + (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(178) P ˇ

^a_k

= Z

^f_k

A

_k

(Z

^f_k

)

^T

(179) A

¹

= I + (HZ)

^T

R

¹

HZ

^f

(180)

P

^a

= ZAZ

^T

(181)

Ensemble transformation

X

^a

= X

^a

+ X

^f_k

W (182)

X

^a_k

= X

^f_k

+ Z

^f_k

W

_k

+ W

_k

⇥

(183) W

_k

= ⌃

w

_k

, . . . , w

_k

⌥

(184) P

^a_k

= 1

N 1 Z

^a_k

(Z

^a_k

)

^T

(185) Z

^a_k

= ⇥

N 1Z

^f_k

A

^1/2_k

(186)

Z

^a_k

= Z

^f_k

W

_k

(187)

W

_k

= ⇥

N 1U

_k

S

_k ^1/2

U

^T_k

(188)

U

_k

S

_k

V

_k

= A

_k ¹

(189)

15 Analysis

X ^f _k = ⌃

x ^f _k ⁽¹⁾ , . . . , x ^f _k ^(N ⁾ ⌥

(167) X ^f _k = ⌃

x ^f _k , . . . , x ^f _k ⌥

(168)

Z ^f _k = X ^f _k X ^f _k (169)

Z = X ^f T (170)

P ^f = ZZ ^T (171)

P ˇ ^f _k = 1

N 1

⇧ N

l=1

⇤ x ^f _k ^(l) x ^f _k ⌅⇤

x ^f _k ^(l) x ^f _k ⌅ T

(172) P ˇ ^f _k = 1

N 1 Z ^f _k ⇤

Z ^f _k ⌅ T

(173) P ˇ ^f _k = Z ^f _k G ⇤

Z ^f _k ⌅ T

(174) G := 1

N 1 I (175)

x ^a _k = x ^f _k + Z ^f _k w _k (176) w _k = A _k (H _k Z ^f _k ) ^T R _k ¹ ⇤

y _k ^o H _k x ^f _k ⌅

(177)

A _k ¹ = G ¹ + (H _k Z ^f _k ) ^T R _k ¹ H _k Z ^f _k (178) A _k ¹ = (N 1)I + (H _k Z ^f _k ) ^T R _k ¹ H _k Z ^f _k (179) P ˇ ^a _k = Z ^f _k A _k (Z ^f _k ) ^T (180) A = G + (HZ) ^T R ¹ HZ ⇥ ¹

(181)

P ^a = ZAZ ^T (182)

Ensemble transformation

X ^a = X ^a + X ^f _k W (183)

X ^a _k = X ^f _k + Z ^f _k W _k + W _k ⇥

(184) W _k = ⌃

w _k , . . . , w _k ⌥

(185) P ^a _k = 1

N 1 Z â _k (Z â _k ) ^T (186) Z â _k = ⇥

N 1Z ^f _k A ^1/2 _k (187)

Z ^a _k = Z ^f _k W _k (188)

W _k = ⇥

N 1U _k S _k ^1/2 U ^T _k (189)

U _k S _k V _k = A _k ¹ (190)

15 Analysis

X

^f_k

= ⌃

x

^f_k⁽¹⁾

, . . . , x

^f_k^(N⁾

⌥

(168) X

^f_k

= ⌃

x

^f_k

, . . . , x

^f_k

⌥

(169)

Z

^f_k

= X

^f_k

X

^f_k

(170)

Z = X X (171)

Z = X

^f

T (172)

P

^f

= ZZ

^T

(173)

P ˇ

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(174) P ˇ

^f_k

= 1

N 1 Z

^f_k

⇤

Z

^f_k

⌅

T

(175) P ˇ

^f_k

= Z

^f_k

G ⇤

Z

^f_k

⌅

T

(176) G := 1

N 1 I (177)

x

^a_k

= x

^f_k

+ Z

^f_k

w

_k

(178) w

_k

= A

_k

(H

_k

Z

^f_k

)

^T

R

_k ¹

⇤

y

_k^o

H

_k

x

^f_k

⌅

(179)

A

_k ¹

= G

¹

+ (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(180) A

_k ¹

= (N 1)I + (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(181) P ˇ

^a_k

= Z

^f_k

A

_k

(Z

^f_k

)

^T

(182) A = G + (HZ)

^T

R

¹

HZ ⇥

¹

(183)

P

^a

= ZAZ

^T

(184)

Ensemble transformation

X

^a

= X

^a

+ X

^f_k

W (185)

X

^a

⇥ ZW (186)

WW

^T

= A (187)

X

^a_k

= X

^f_k

+ Z

^f_k

W

_k

+ W

_k

⇥

(188) W

_k

= ⌃

w

_k

, . . . , w

_k

⌥

(189)

15 Transformation matrix in ensemble space (small matrix)

ETKF:

A has dimension N

²

G = I (identity matrix) SEIK:

A has dimension (N-1)

²

G = ( T T

^T

)

^-1

Analysis

X

^f_k

= ⌃

x

^f_k⁽¹⁾

, . . . , x

^f_k^(N⁾

⌥

(167) X

^f_k

= ⌃

x

^f_k

, . . . , x

^f_k

⌥

(168)

Z

^f_k

= X

^f_k

X

^f_k

(169)

Z = X

^f

T (170)

P

^f

= ZZ

^T

(171)

P ˇ

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(172) P ˇ

^f_k

= 1

N 1 Z

^f_k

⇤

Z

^f_k

⌅

T

(173) P ˇ

^f_k

= Z

^f_k

G ⇤

Z

^f_k

⌅

T

(174) G := 1

N 1 I (175)

x

^a_k

= x

^f_k

+ Z

^f_k

w

_k

(176) w

_k

= A

_k

(H

_k

Z

^f_k

)

^T

R

_k ¹

⇤

y

^o_k

H

_k

x

^f_k

⌅

(177)

A

_k ¹

= G

¹

+ (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(178) A

_k ¹

= (N 1)I + (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(179) P ˇ

^a_k

= Z

^f_k

A

_k

(Z

^f_k

)

^T

(180) A = G + (HZ)

^T

R

¹

HZ ⇥

1

(181)

P

^a

= ZAZ

^T

(182)

Ensemble transformation

X

^a

= X

^a

+ X

^f_k

W (183)

X

^a_k

= X

^f_k

+ Z

^f_k

W

_k

+ W

_k

⇥

(184) W

_k

= ⌃

w

_k

, . . . , w

_k

⌥

(185) P

^a_k

= 1

N 1 Z

^a_k

(Z

^a_k

)

^T

(186) Z

^a_k

= ⇥

N 1Z

^f_k

A

^1/2_k

(187)

Z

^a_k

= Z

^f_k

W

_k

(188)

W

_k

= ⇥

N 1U

_k

S

_k ^1/2

U

^T_k

(189)

U

_k

S

_k

V

_k

= A

_k ¹

(190)

15 Analysis state covariance matrix

Analysis

X

^f_k

= ⌃

x

^f_k⁽¹⁾

, . . . , x

^f_k^(N⁾

⌥

(167) X

^f_k

= ⌃

x

^f_k

, . . . , x

^f_k

⌥

(168)

Z

^f_k

= X

^f_k

X

^f_k

(169)

Z = X

^f

T (170)

P ˇ

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(171) P ˇ

^f_k

= 1

N 1 Z

^f_k

⇤

Z

^f_k

⌅

T

(172) P ˇ

^f_k

= Z

^f_k

G ⇤

Z

^f_k

⌅

T

(173) G := 1

N 1 I (174)

x

^a_k

= x

^f_k

+ Z

^f_k

w

_k

(175) w

_k

= A

_k

(H

_k

Z

^f_k

)

^T

R

_k ¹

⇤

y

_k^o

H

_k

x

^f_k

⌅

(176)

A

_k ¹

= G

¹

+ (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(177) A

_k ¹

= (N 1)I + (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(178) P ˇ

^a_k

= Z

^f_k

A

_k

(Z

^f_k

)

^T

(179)

A

¹

= I + (HZ)

^T

R

¹

HZ (180)

P

^a

= ZAZ

^T

(181)

Ensemble transformation

X

^a

= X

^a

+ X

^f_k

W (182)

X

^a_k

= X

^f_k

+ Z

^f_k

W

_k

+ W

_k

⇥

(183) W

_k

= ⌃

w

_k

, . . . , w

_k

⌥

(184) P

^a_k

= 1

N 1 Z

^a_k

(Z

^a_k

)

^T

(185) Z

^a_k

= ⇥

N 1Z

^f_k

A

^1/2_k

(186)

Z

^a_k

= Z

^f_k

W

_k

(187)

W

_k

= ⇥

N 1U

_k

S

_k ^1/2

U

^T_k

(188)

U

_k

S

_k

V

_k

= A

_k ¹

(189)

15 Computations in ensemble-spanned space

Ensemble transformation based on square root of A

Very efficient:

Transformation matrix computed in space of dim. N or N-1

Analysis

X

^f_k

= ⌃

x

^f_k⁽¹⁾

, . . . , x

^f_k^(N⁾

⌥

(167) X

^f_k

= ⌃

x

^f_k

, . . . , x

^f_k

⌥

(168)

Z

^f_k

= X

^f_k

X

^f_k

(169)

Z = X

^f

T (170)

P

^f

= ZZ

^T

(171)

P ˇ

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

^T

(172) P ˇ

^f_k

= 1

N 1 Z

^f_k

⇤

Z

^f_k

⌅

T

(173) P ˇ

^f_k

= Z

^f_k

G ⇤

Z

^f_k

⌅

T

(174) G := 1

N 1 I (175)

x

^a_k

= x

^f_k

+ Z

^f_k

w

_k

(176) w

_k

= A

_k

(H

_k

Z

^f_k

)

^T

R

_k ¹

⇤

y

_k^o

H

_k

x

^f_k

⌅

(177)

A

_k ¹

= G

¹

+ (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(178) A

_k ¹

= (N 1)I + (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(179) P ˇ

^a_k

= Z

^f_k

A

_k

(Z

^f_k

)

^T

(180) A = G + (HZ)

^T

R

¹

HZ ⇥

¹

(181)

P

^a

= ZAZ

^T

(182)

Ensemble transformation

X

^a

= X

^a

+ X

^f_k

W (183)

X

^a

⇥ X

^f

L (184)

LL

^T

= A (185)

X

^a_k

= X

^f_k

+ Z

^f_k

W

_k

+ W

_k

⇥

(186) W

_k

= ⌃

w

_k

, . . . , w

_k

⌥

(187)

15 Analysis

X

^f_k

= ⌃

x

^f_k⁽¹⁾

, . . . , x

^f_k^(N⁾

⌥

(167) X

^f_k

= ⌃

x

^f_k

, . . . , x

^f_k

⌥

(168)

Z

^f_k

= X

^f_k

X

^f_k

(169)

Z = X

^f

T (170)

P

^f

= ZZ

^T

(171)

P ˇ

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(172) P ˇ

^f_k

= 1

N 1 Z

^f_k

⇤

Z

^f_k

⌅

T

(173) P ˇ

^f_k

= Z

^f_k

G ⇤

Z

^f_k

⌅

T

(174) G := 1

N 1 I (175)

x

^a_k

= x

^f_k

+ Z

^f_k

w

_k

(176) w

_k

= A

_k

(H

_k

Z

^f_k

)

^T

R

_k ¹

⇤

y

^o_k

H

_k

x

^f_k

⌅

(177)

A

_k ¹

= G

¹

+ (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(178) A

_k ¹

= (N 1)I + (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(179) P ˇ

^a_k

= Z

^f_k

A

_k

(Z

^f_k

)

^T

(180) A = G + (HZ)

^T

R

¹

HZ ⇥

¹

(181)

P

^a

= ZAZ

^T

(182)

Ensemble transformation

X

^a

= X

^a

+ X

^f_k

W (183)

X

^a

⇥ ZL (184)

LL

^T

= A (185)

X

^a_k

= X

^f_k

+ Z

^f_k

W

_k

+ W

_k

⇥

(186) W

_k

= ⌃

w

_k

, . . . , w

_k

⌥

(187)

15

(17)

The T matrix

SEIK and ETKF use different projections T

➜  results in slightly different ensemble transformations

➜  SEIK is slightly faster than ETKF

ETKF uses minimal ensemble transformation – desirable feature!

Analysis

X

^f_k

= ⌃

x

^f_k⁽¹⁾

, . . . , x

^f_k^(N⁾

⌥

(167) X

^f_k

= ⌃

x

^f_k

, . . . , x

^f_k

⌥

(168)

Z

^f_k

= X

^f_k

X

^f_k

(169)

Z = X

^f

T (170)

P ˇ

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(171) P ˇ

^f_k

= 1

N 1 Z

^f_k

⇤

Z

^f_k

⌅

T

(172) P ˇ

^f_k

= Z

^f_k

G ⇤

Z

^f_k

⌅

T

(173) G := 1

N 1 I (174)

x

^a_k

= x

^f_k

+ Z

^f_k

w

_k

(175) w

_k

= A

_k

(H

_k

Z

^f_k

)

^T

R

_k ¹

⇤

y

_k^o

H

_k

x

^f_k

⌅

(176)

A

_k ¹

= G

¹

+ (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(177) A

_k ¹

= (N 1)I + (H

_k

Z

^f_k

)

^T

R

_k ¹

H

_k

Z

^f_k

(178) P ˇ

^a_k

= Z

^f_k

A

_k

(Z

^f_k

)

^T

(179) A

¹

= I + (HZ)

^T

R

¹

HZ

^f

(180)

P

^a

= ZAZ

^T

(181)

Ensemble transformation

X

^a

= X

^a

+ X

^f_k

W (182)

X

^a_k

= X

^f_k

+ Z

^f_k

W

_k

+ W

_k

⇥

(183) W

_k

= ⌃

w

_k

, . . . , w

_k

⌥

(184) P

^a_k

= 1

N 1 Z

^a_k

(Z

^a_k

)

^T

(185) Z

^a_k

= ⇥

N 1Z

^f_k

A

^1/2_k

(186)

Z

^a_k

= Z

^f_k

W

_k

(187)

W

_k

= ⇥

N 1U

_k

S

_k ^1/2

U

^T_k

(188)

U

_k

S

_k

V

_k

= A

_k ¹

(189)

15

(18)

Error Subspace Transform Kalman Filter (ESTKF)

Combine advantages of SEIK and ETKF

Redefine T:

1.  Remove ensemble mean from all columns 2.  Subtract fraction of last column from all others 3.  Drop last column

L. Nerger et al., Monthly Weather Review 140 (2012) 2335-2345

Features of the ESTKF:

•  Same ensemble transformation as ETKF

•  Slightly cheaper computations

•  Direct access to ensemble-spanned error space

(19)

Requirements for applying ensemble Kalman filters

“Pure” ensemble-based Kalman filters have usually bad performance

•  e.g. due to small ensemble size

Improvements through

•  Covariance inflation

•  Localization

•  Model error simulation

S: Analysis region

D: Corresponding data region Localization

(20)

Implementation Aspects

(21)

Computational and Practical Issues

Data assimilation with ensemble-based Kalman filters is costly!

Memory: Huge amount of memory required (model fields and ensemble matrix)

Computing: Huge requirement of computing time (ensemble integrations)

Parallelism: Natural parallelism of ensemble integration exists (needs to be implemented)

„Fixes “ : Filter algorithms do not work in their pure form („fixes “ and tuning are needed)

because Kalman filter optimal only in linear case

(22)

Implementing Ensemble Filters & Smoothers

Ensemble forecast

•  can require model error simulation

•  naturally parallel

Analysis step of filter algorithms operates on abstract state vectors (no specific model fields)

Analysis step requires information on observations

•  which field?

•  location of observations

•  observation error covariance matrix

•  relation of state vector to observation

➜ Analysis step can be implemented independently of model!

(23)

Offline mode – separate programs Model

Aaaaaaaa Aaaaaaaa aaaaaaaa a

Start

Stop

read ensemble files analysis step

Aaaaaaaa Aaaaaaaa aaaaaaaaa

Start

Stop Do i=1, nsteps

Initialize Model

generate mesh Initialize fields

Time stepper

consider BC Consider forcing

Post-processing

For each ensemble state

•  Initialize from restart files

•  Integrate

•  Write restart files

•  Read restart files (ensemble)

•  Compute analysis step

•  Write new restart files

Assimilation program

write model restart files

⬅ generic

(24)

Extending a Model for Data Assimilation

Aaaaaaaa Aaaaaaaa aaaaaaaaa

Start

Stop Do i=1, nsteps

Initialize Model

Time stepper

Post-processing

Aaaaaaaa Aaaaaaaa aaaaaaaaa

Start

Stop Do i=1, nsteps

Initialize Model

Time stepper

Post-processing

Model

Ensemble forecast

Analysis step

Initialization Extension for data assimilation

Aaaaaaaa Aaaaaaaa aaaaaaaaa

Start

Stop

Initialize Model

Time stepper

Post-processing init_parallel_asml

Do i=1, nsteps get_state_asml

init_asml

put_state_asml Filter-Analysis

Online assimilation program:

➜ Avoid expensive writing and reading of ensemble files

(25)

Features of online program

•  minimal changes to model code when combining model with filter algorithm (adding 4 routines)

•  model not required to be a subroutine

•  no change to model numerics

•  control of assimilation program coming from model

•  filter method encapsulated in subroutine

•  simple switching between different filters and data sets

•  complete parallelism in model, filter, and ensemble integrations

Implementation structure can be implemented in a generic

framework (for online and offline modes)

(26)

PDAF: A tool for data assimilation

PDAF - Parallel Data Assimilation Framework

  an environment for ensemble assimilation

  provide support for ensemble forecasts

  provide fully-implemented filter algorithms

  for testing algorithms and real applications

  useable with virtually any numerical model

  makes good use of supercomputers

Open source:

Code and documentation available at http://pdaf.awi.de

L. Nerger, W. Hiller, Computers & Geosciences 55 (2013) 110-118

Using Ensemble Kalman Filters to Assimilate Dynamic Ocean Topography Data into a Global Ocean Model