Transcript Examples of Four-Dimensional Data Assimilation in Oceanography Ibrahim Hoteit University of Maryland October 3, 2007

Examples of Four-Dimensional Data
Assimilation in Oceanography
Ibrahim Hoteit
University of Maryland
October 3, 2007
Outline

4D Data Assimilation



4D-VAR and Kalman Filtering
Application to Oceanography
Examples in Oceanography

4D-VAR Assimilation
Tropical Pacific, San Diego, …

Filtering Methods
Mediterranean Sea, Coupled models, Nonlinear filtering ...

Discussion and New Applications
1
Data Assimilation

Goal: Estimate the state of a dynamical system

Information:
 Imperfect dynamical Model:
xk 1  M k 1,k ( xk ) k 1
x state vector, 

N(0,Q) model error
M k 1,k transition operator form k to k  1

Sparse observations:

yko  Hk ( xk )   k

y observation vector,  N(0, R) observation error
H k observational operator

A priori Knowledge: xb and its uncertainties B

o
2
Data Assimilation

Data assimilation: Use all available information to determine the
best possible estimate of the system state



Observations show the real trajectory to the model
Model dynamically interpolates the observations
a
3D assimilation: Determine an estimate x of the state at a
given time given an observation y o by minimizing
T
J ( x)   y  Hx  R 1  y o  Hx    xb  x  B1  xb  x 
o

4D assimilation: Determine x1a , , xNa given y1o ,
 4D-VAR and Kalman Filtering
, yNo
3
4D-VAR Approach

Optimal Control: Look for the model trajectory that best fits the
observations by adjusting a set of “control variables”  minimize
J (c)   k 1  y  H k xk  R 1  yko  H k xk   cb  c  B1 cb  c 
N
o
k
T
with the model as constraint: xk 1  M k 1,k ( xk ) k 1

c is the control vector and may include any model parameter
(IC, OB, bulk coefficients, etc) … and model errors
 Use
a gradient descent algorithm to minimize J
 Most efficient way to compute the gradients is to run the
adjoint model backward in time
4
Kalman Filtering Approach

Bayesian estimation: Determine pdf of xk given y1:ok   y1o , , yko 
Px
k

y1o:k
Minimum Variance (MV) estimate (minimum error on average)
xka  E  xk / y1:ok 

Maximum a posteriori (MAP) estimate (most likely)
xka  Arg max Px
k

y1:ok
Kalman filter (KF): provides the MV (and MAP) estimate for linearGaussian Systems
5
The Kalman Filter (KF) Algorithm
 Initialization Step: x0 and P0
Analysis Step (observation)
Forecast Step (model)
 Forecast state
xkf  M k ,k 1xka1
 Forecast Error covariance
Pkf  Mk ,k 1 Pka1 Mkk1  Qk
 Kalman Gain
Gk  Pkf HkT [ Hk Pkf HkT  R k ]1
 Analysis state
xka  xkf  Gk [ yko  Hk xkf ]
 Analysis Error covariance
Pka  Pkf  Gk Hkf Pkf
6
Application to Oceanography

4D-VAR and the Kalman filter lead to the same estimate at the end of
the assimilation window when the system is linear, Gaussian and perfect

Nonlinear system:
 4D-VAR cost function is non-convex  multiple minima
 Linearize the system  suboptimal Extended KF (EKF)

System dimension ~ 108:
 4D-VAR control vector is huge
 KF error covariance matrices are prohibitive

Errors statistics:
 Poorly known R, Q, B
 Non-Gaussian: KF is still the MV among linear estimators
7
4D Variational Assimilation

ECCO 1o Global Assimilation System

Eddy-Permitting 4D-VAR Assimilation

ECCO Assimilation Efforts at SIO

Tropical Pacific, San Diego, …
In collaboration with the ECCO group, especially
Armin Köhl*, Detlef Stammer*, Patrick Heimbach**
*Universitat Hamburg/Germany, **MIT/USA
8
ECCO 1o Global Assimilation System

Model:



Data:




MITGCM (TAF-compiler enabled)
NCEP forcing and Levitus initial conditions
Altimetry (daily): SLA TOPEX, ERS
SST (monthly): TMI and Reynolds
Profiles (monthly) : XBTs, TAO, Drifters, SSS, ...
Climatology (Levitus S/T) and Geoid (Grace mission)

Assimilation scheme: 4D-VAR with control of the initial
conditions and the atmospheric forcing (with diagonal weights!!!)

ECCO reanalysis: 1o global ocean state and atmospheric forcing from
1992 to 2004, …and from 1952  2001 (Stammer et al. …)
9
ECCO Solution Fit
Equatorial Under Current (EUC)
Johnson
ECCO
10
ECCO Tropical Pacific Configuration

Regional: 26oS  26oN, 1/3o, 50 layers, ECCO O.B.

Data:
 TOPEX, TMI SST, TAO, XBT, CTD, ARGO, Drifters; all
at roughly daily frequency
MITGCM

Tropical
Climatology: Levitus-T
and S,Pacific
Reynolds SST and GRACE

Control: Initial conditions , 2-daily forcing, and weekly O.B.

Smoothing: Smooth ctrl fields using Laplacian in the horizontal
and first derivatives in the vertical and in time

First guess: Levitus (I.C.), NCEP (forcing), ECCO (O.B.)
OB = (U,V,S,T)
11
Eddy-Permitting 4D-VAR Assimilation

The variables of the adjoint model exponentially increase in time
 Typical behavior for the adjoint of a nonlinear chaotic model
 Indicate unpredictable events and multiple local minima
 Correct gradients but wrong sensitivities
 Invalidate the use of a gradient-based optimization algorithm

Assimilate over short periods (2 months) where the adjoint is stable
Replace the original unstable adjoint with the adjoint of a tangent linear
model which has been modified to be stable (Köhl et al., Tellus-2002)


Exponentially increasing gradients were filtered out using
larger viscosity and diffusivity terms in the adjoint model
12
HFL gradients after 45 days with increasing viscosities
Visc = 1e11 & Diff = 4e2
20*Visc & 20*Diff
10*Visc & 10*Diff
30*Visc & 30*Diff
13
Initial temperature gradients after 1 year (2000)
10*Visc & 10*Diff
20*Visc & 20*Diff
14
Data Cost Function Terms
1/3;39
1/6;39
1;39
1;23
15
Control Cost Function Terms
1/3;39
1/6;39
1;39
1;23
16
Fit to Data
1/3;39
1/6;39
1;39
1;23
17
Assimilation Solution (weekly field end of August)
18
What Next …

Fit is quite good and assimilation solution is reasonable

Extend assimilation period over several years

Add new controls to enhance the controllability of the
system and reduce errors in the controls

Improve control constraints …

Some references
 Hoteit et al. (QJRMS-2006)
 Hoteit et al. (JAOT-2007)
 Hoteit et al. (???-2007)
19
Other MITGCM Assimilation Efforts at SIO

1/10o CalCOFI 4D-VAR assimilation system
Predicting the loop current in the Gulf of Mexico …

San Diego high frequency CODAR assimilation

●
Assimilate hourly HF radar data and other data
Adjoint effectiveness at small scale
 Information content of surface velocity data

●
MITGCM with 1km resolution and 40 layers
●
Control: I.C., hourly forcing and O.B.
●
First guess: one profile T, S and TAU (no U, V, S/H-FLUX)
●
Preliminary results: 1 week, no tides
20
Model Domain and
Radars Coverage
Time evolution of the
normalized radar cost
1/3;39
1/6;39
1;39
1;23
21
Assimilation Solution: SSH / (U,V) & Wind Adj.
1/3;39
1/6;39
1;39
1;23
22
What Next …

Assimilation over longer periods

Include tidal forcing

Coupling with atmospheric model

Nesting into the CalCOFI model
23
Filtering Methods

Low-Rank Extended/Ensemble Kalman Filtering

SEEK/SEIK Filters

Application to Mediterranean Sea

Kalman Filtering for Coupled Models

Particle Kalman Filtering
In collaboration with
D.-T. Pham*, G. Triantafyllou**, G. Korres**
*CNRS/France, **HCMR/Greece
24
Low-rank Extended/Ensemble Kalman Filtering

Reduced-order Kalman filters: Project x on a low-dim subspace
xk  Lxk  Pk  L Pk L
T
 x (n); x (r )

 L (n, r ); P ( r , r )
 Kalman correction along the directions of L

Reduced error subspace Kalman filters: Pk has low-rank
Pk  Lk Pk Lk T

Ensemble Kalman filters: Monte Carlo approach to
Pk 
1
N
i
i
i
i T
(
x

x
)(
x

x
 i1 k k k k )
N
 Correction along the directions of [ x1k  xki
x1k  xkN ]
25
Singular Evolutive Kalman (SEEK) Filters
 Low-rank (r) error subspace Kalman filters: x0 ,P0  L0U0 LT0
Analysis
Forecast
x  M k (x
f
k
a
k 1
Lk  M k Lk 1

)
1
k
U U
1
k 1
  Hk Lk  R Hk Lk
1
i
T
x  x  LkU k  Hk Lk  R k1[ yko  H k ( xkf )]
a
k
f
k
T
A “collection” of SEEK filters:
 SEEK: Extended variant M  Μ
 SFEK: Fixed variant M  Id  L  L0
 SEIK: Ensemble variant with (r+1) members only! (~ETKF)
 Inflation and Localization
26
The Work Package WP12 in MFSTEP

EU project between several European institutes

Assimilate physical & biological observations into coupled
ecosystem models of the Mediterranean Sea:

Develop coupled physical-biological model for regional
and coastal areas of the Mediterranean Sea

Implement Kalman filtering techniques with the physical
and biological model

… Investigate the capacity of surface observations (SSH,
CHL) to improve the behavior of the coupled system
27
The Coupled POM – BFM Model
One way coupled: Ecology
does not affect the physics
x p  n, u, v, u, v, T, S, q2 , q2l
O2o,O3o, N1p, N3n, N5s, N4n,



e
x  P1,2,3,4 (c,n,p,s,i),Z4,5,6 (c,n,p),R 6(c,n,p,s), 


 R1(c,n,p),B1(c,n,p),Q1,6,7 (c,n,p,s)

28
A Model Snapshot
1/10o Eastern Mediterranean configuration 25 layers
Elevation and Mean Velocity
Mean CHL integrated 1-120m
29
Assimilation into POM

Model = 1/10o Mediterranean configuration with 25 layers

Observations = Altimetry, SST, Profiles T & S profiles,
Argo data, and XBTs on a weekly basis

SEIK Filter with rank 50 (51 members)

Initialization = EOFs computed from 3-days outputs of a
3-year model integration

Inflation factor = 0.5

Localization = 400 Km
30
Assimilation into POM
Mean Free-run RMS Error
SSH RMS Misfits
Mean Forecast RMS Error
Free-Run
Forecast
Obs Error = 3cm
Mean Analysis RMS Error
Analysis
31
Salinity RMS Error Time Series
FerryBox data at Rhone River
SSH 07/12/2005
07/12/05: SATELLITE SSH
FREE RUN
FORECAST
ANALYSIS
6/11/2015
ECOOP KICK-OFF
32
Assimilation into BFM

Model = 1/10o Eastern Mediterranean with 25 layers
with perfect physics

Observations = SeaWiFS CHL every 8 days in 1999

SFEK Filter = SEEK with invariant correction subspace

Correction subspace = 25 EOFs computed from 2-days
outputs of a one year model integration

Inflation factor = 0.3

Localization = 200 Km
33
Assimilation into BFM
CHL RMS Misfits
Forecast
Free-Run
Analysis
34
CHL Cross-Section at 34oN
Ph Cross-Section at 28oE
35
Kalman Filtering for Coupled Models
Physical System

p
p
p
p

 xk  M k ( xk 1 )   k
 p
p
p
p
y

H
(
x
)



k
k
k
 k


Ecological System
e
e
e
p
e

x

M
(
x
,
x
)


 k
k
k 1
k
k
 e
e
e
e

 yk  H k ( xk )   k
MAP: Direct maximization of the joint conditional density
 xp 
 e   Max Px p xe y p ye
x 
 standard Kalman filter estimation problem
 Joint approach: strong coupling and same filter (rank) !!!
Dual Approach
Decompose the joint density into marginal densities

Px p xe

y p ye
 Pxe
P
x p y p ye x p y p ye
Compute MAP estimators from each marginal density
x p  Arg max Px p
y p ye
&
x e  Arg max Pxe
x p y p ye
 Separate optimization leads to two Kalman filters …
 Different degrees of simplification and ranks for each filter
 significant cost reduction
 Same
e
x from the joint or the dual approach
p
e
 The physical filter assimilate y and y
RRMS for state vectors

Twin-Experiments 1/10o Eastern Mediterranean (25 layers)

Joint: SEEK rank-50

 RRMS
.
X REF  X FILTER

Dual: SEEK rank-50 for physics
SFEK rank-20 for biology
X REF  X
Physics
Biology
Ref
Dual
Joint
Ref
Dual
Joint
38
What Next …

Joint/Dual Kalman filtering with real data

State/Parameter Kalman estimation

Better account for model errors

Some references

Hoteit et al. (JMS-2003), Triantafyllou et al. (JMS-2004),
Hoteit et al. (NPG-2005), Hoteit et al. (AG-2005),
(Hoteit et al., 2006), Korres et al. (OS-2007)
39
Nonlinear Filtering - Motivations

The EnKF is “semi-optimal”; it is analysis step is linear

The optimal solution can be obtained from the optimal nonlinear
filter which provides the state pdf given previous data

Particle filter (PF) approximates the state pdf by mixture of Dirac
N
functions i 1 wi x , but suffers from the collapse (degeneracy)
of its particles (analysis step only update the weights wi )
i

Surprisingly, recent results suggest that the EnKF is more stable
than the PF for small ensembles because the Kalman correction
attenuates the collapse of the ensemble
40
The Particle Kalman Filter (PKF)

The PKF uses a Kernel estimator to approximate the pdfs of the
nonlinear filter by a mixture of Gaussian densities
Px / y
1:k 1

 i1 wki k 1 G  x  xki k 1; ik k 1 
N
The state pdfs can be always approximated by mixture of
Gaussian densities of the same form:
 Analysis Step:
Kalman-type: EKF analysis to update x i and Pi
Particle-type: weight update (but using  instead of R )
i
 Forecast Step: EKF forecast step to propagate x and P i

Resampling Step: …
41
Particle Kalman Filtering in Oceanography

It is an ensemble of extended Kalman filters with weights!!
 Particle Kalman Filtering
 requires simplification of the particles error covariance matrices


The EnKF can be derived as a simplified PKF
Hoteit et al. (MWR-2007) successfully tested one low-rank
PKF with twin experiments

What Next …
 Derive
and test several simplified variants of the PKF
 Assess the relevance of a nonlinear analysis step:
 comparison with the EnKF
 Assimilation of real data …
42
Discussion and New Applications

Advanced 4D data assimilation methods can be now applied to
complex oceanic and atmospheric problems

More work is still needed for the estimation of the error
covariance matrices, the assimilation into coupled models, and
the implementation of the optimal nonlinear filter

New Applications:
 ENSO prediction using neural models and Kalman filters
 Hurricane reconstruction using 4D-VAR ocean assimilation!
 Ensemble sensitivities and 4D-VAR
 Optimization of Gliders trajectories in the Gulf of Mexico
 …
THANK YOU
43
44
4D-VAR or (Ensemble) Kalman Filter?
4D-VAR
Easier to understand
EnKF
Easier to incorporate a complex
background covariance matrix
More portable
(easier to implement?)
Support different degrees of
simplifications
Low-rank estimates of the error
cov. matrices (better forecast!)
Dynamically consistent solution
Still room for improvement …
No low-rank deficiency
 4D-VAR or EnkF? …
45
Sensitivity to first guess (25 Iterations)
46
Comparison with TAO-Array RMS
RMS Zonal Velocity (m/s)
1/3;39
1/6;39
RMS Meridional Velocity (m/s)
1;39
1;23
47
San Diego HF Radar Currents Assimilation

Assimilate hourly HF radar data and other data

Goals:
 Adjoint effectiveness at small scale
 Information content of surface velocity data
 Dispersion of larvae, nutrients, and pollutants

MITGCM with 1km resolution with 40 layers

Control: I.C., hourly forcing, and O.B.

First guess: one profile, no U and V, and no forcing!

Preliminary results: 1 week, no tides
48
Cost Function terms
1/3;39
1/6;39
1;39
1;23
49
Assimilation Solution: U & V
1/3;39
1/6;39
1;39
1;23
50
Assimilation Solution
1/3;39
1/6;39
1;39
1;23
51
Gulf of Mexico

Loop current prediction
 Observations HF radar, Gliders, ADCP, …

Adjoint effectiveness …

1/10o with 50 layers

Ctrl: I.C. (S,T), daily forcing, and weekly O.B.

Proof of concept assimilating SSH, Levitus and
Reynolds …
52
What Next …

Ensemble forecasting

Ensemble Kalman filtering

Optimization of observation systems
53
Twin-Experiments Setup
Model = 1/10o Eastern Mediterranean with 25 layers

1996
2000
Spin up
4 years
05/03/02
EOFs
2 years
05/03/02
REF – OBS
3 months

Pseudo-obs: SSH and CHL surface data every 3 days

Initialization: start from mean state of the 2 years run

Free-run: run without assimilation starting from mean state

Evaluation: RMS misfit relative to the misfit from mean state
RRMS 
X REF  X FILTER
X REF  X
54
Low-rank Extended/Ensemble Kalman Filtering
EnKF
a ,i a low-dim
f ,i
f ,i
 Reduced-order Kalmanf ,ifilters:
subspace
● Project xxon
●

x

G
[
y

Hx
xk
k
k
k ]
●
●
x T(r ) f T
1
x (Pnf);H
T
xka,i1 ● ● xModel
G
[H
P
H

R
]
●
 Sxk ● ●Pk 
S●Pk S

k
●
●
●
●
a
n, r ); xPa ,(i r), r )
● ●
 S(cov(
●
P
●● ●
● ●●
k
● Data
●
●
●
●
●
of S
● ● ●Analysis along the directions
● ●
SPF
●
f
f ,i
a
f
f ,i

x

x
●  Reduced
k
k
error
subspace
Kalman
filters:
has
low-rank
x

x

G
[
y

Hx
●
k Pk
k
k ]
 f
f ,i
f
T
f
T
1
P

cov(
x
)
T
G

P
H
[H
P
H

R
]
 Pk  S P Sk
k
k k k
Resampling
Pa  P f  G H P f
 Ensemble Kalman filters: Monte Carlo approach to  Pk
Pk 

1
N
i
i
i
i T
(
x

x
)(
x

x
 i1 k k k k )
N
Analysis along the directions of [ x1k  xki
x1k  xkN ]
55
Low-Rank Deficiency

Issues
 Error covariance matrices are underestimated
 Few degrees of freedom to fit the data

Amplification by an inflation factor
P P
f

f
Localization
of the covariance matrix (using Schur product)
56
Joint Approach
 Direct maximization of the joint conditional density
standard Kalman filter estimation problem acting on
 xkp   M kp ( xkp1 )  kp 
xk   e    e e    e 
 x   M ( x )   
 k   k k 1   k 
and assimilating
 ykp   H kp ( xkp )    kp 
yk   e    e e    e 
 y   H (x )    
 k   k k   k 
 Issues: Strong coupling and same filter (rank)
57
Dual Approach – Some Facts
e
 Only the second marginal density depends on x , this means
same x e from the joint or the dual approach
p
e
 x does not depend on x , more in line with the one-way
coupling of the system
p
e
 The physical filter assimilates both y and y : assimilation
of y e guaranties consistency between the two subsystems
e
 The ecological filter assimilates only y , but it is forced with
the solution of the physical filter
 The linearization of the observational operator in the physical
filter is a complex operation because of the dependency of the
ecology on x p, it was neglected in this preliminary application
Twin-Experiments Setup
Model = 1/10o Eastern Mediterranean with 25 layers

1996
2000
Spin up
4 years
05/03/02
EOFs
2 years
05/03/02
REF – OBS
3 months

Pseudo-obs: SSH and CHL surface data every 3 days

Initialization: start from mean state of the 2 years run

Free-run: run without assimilation starting from mean state

Evaluation: RMS misfit relative to the misfit from mean state
RRMS 
X REF  X FILTER
X REF  X
59
The Optimal Nonlinear Filter

As the Kalman filter, it operates as a succession of forecast and
analysis steps to update the state pdf:

Forecast Step: Integrate the analysis pdf with the model


pk k 1 ( x | y1:k 1 )   n  x  M k (u ); Qk pk 1 (u | y1:k 1 )du

Analysis Step: Correct the predictive pdf with the new data
pk ( x | y1:k ) 

1
bk
pk k 1 ( x | y1:k 1 )   yk  H k ( x); Rk 
Particle Filter approximates the state pdf by mixture of Dirac
functions  N wi x , but suffers from degeneracy.
i 1
i
60
New Directions/Applications

New Applications
 ENSO prediction using surrogate models and Kalman filters
 Hurricane reconstruction using 4DVAR ocean assimilation!
 Ensemble sensitivities and 4DVAR

Other Interests

Optimal Observations

Estimate Model and Observational Errors

Estimate Background Covariance Matrices in 4DVAR

Study the behavior of the different 4DVAR methods with
highly nonlinear models
61