Data Assimilation- A11. 4DVAR Shunji Kotsuki
Environmental Prediction Science Laboratory
Center for Environmental Remote Sensing (CEReS), Chiba University(shunji.kotsuki@chiba-u.jp)
with many thanks to Mrs. F. Kawasaki and T. SaitoDA Lectures A (Basic Course)‣ (1) Introduction and NWP‣ (2) Deterministic Chaos and Lorenz-96 model‣ (3) A toy model and Bayesian estimation‣ (4) Kalman Filter (KF)
‣ (5) 3D Variational Method (3DVAR)‣ (6) Ensemble Kalman Filter (PO method)‣ (7) Serial Ens. Square Root Filter (Serial EnSRF)‣ (8) Local Ens. Transform Kalman Filter (LETKF)‣ (9) Innovation Statistics‣ (10) Adaptive Inflations‣ (11) 4D Variational Method (4DVAR)Today’s goals‣ Lecture‣ what is the 4D-Var?
‣ what is the cost function of 4DVAR?
‣ what is adjoint and back propagation?

‣ Training Course‣ to implement 4DVAR

Review: 3DVAR3DVAR Equations
Kalman Gain

Analysis Error Covariance

Analysis Update Equation

4DVAR

Conceptual Image of 4DVAR
4D-VarObservations

Forecastw DAFirst Guess

Forecast
w/o DA
0Cost Function

time

Data Assimilation Window
RegularizationMisfit to 4d Obs||

incremental form
|

|

=|

: nonlinear model from t=0 to t=i,

: nonlinear obs. operator at t=i|

4DVAR Equation
Tangent Linear Model

Cost Function (scalar)|

Jacobian ()numerical model||obs operator|

backpropagations
Hessian (×)

Problem to be solved|

Nonlinear, TLM & ADJ models
UpdatedtrajectoryFirst GuessBackgroundtrajectory
t=0
Nonlinear Models:

t=i|

,

Tangent Linear Models (TLM):

Adjoint Model (ADJ):| ,| ,Update neededtimeFlow-dependent B in 4DVAR
Analysis of 3DVARassimilating obs at time kstatic background (is used0Analysis of 4DVARk||||
||

||

|

|

|||||

|

flow-dependent background (|||is usedtime

Iterative Solver(BFGS method)

How to solve 4DVAR?
Newton method(wiki)

Problem to be solved

Newton Method
after jth iteration (), to obtainthat satisfies:

Talyer expansion gives

And obtain

whenHere, Hessian matrix (0×to update) is given by Quasi-Newton method:(1) without having
explicitly(2) to approximate
by a positive definite matrix

for inversion

iteratively

Quasi-Newton Method
Secant Conditions
Suppose we haveLet

be a positive definite matrix thatapproximates the Hessian matrix (i.e.,, and would like to updateTaylor)for the next iteration
Approximation

Secant Equation
Secant condition

whereThere are many matrices that satisfy the Secant Condition (e.g., DFP, BFGS, SR1).
BFGS method

Usually(widely used in 4DVAR for NWP)

Algorithm of BFGS method1

to set initial condition (2to compute the update direction by3

to compute the update parameterand)

by Armijo Condition4to update5

to update6

to repeat steps 2-5 until

𝐐 𝐬 𝐐 𝐬(𝐬 ) 𝐐 𝐬𝐮 (𝐮 )(𝐬 ) 𝐮e.g.

iteration(there may be other efficient condition such as Wolfe’s condition)

Armijo Condition
To obtain an appropriate update parameter

Armijo Condition

１．set parameters

Algorithm,

２．set initial condition:
３．end the algorithm if
４．If not, update

satisfies the condition
and go back Step 3

Outer and Inner Loops

ECMWF 4DVAR (2015)
IFS Documentation (2015)

JMA 4DVAR (2022)

NL: Non linear
TL: Tangent linear model
AD: Adjoint model
https://www.jma.go.jp/jma/jma-eng/jma-center/nwp/outline2022-nwp/pdf/outline2022_02.pdf

Solver of JMA 4DVAR (2022)‣ The limited memory Broyden-FletcherGoldfarb-Shanno (L-BFGS) algorithm (Liu andNocedal 1989)‣ with Veersé’s preconditioner (Veersé et al. 2000)

to be updated

JMS’s NWP System (2022)

JMA (2022)

Outer loop: High-res. NL

Outer loop: High- & Low-res. NL(O1) high-res. NL model(O2) departure:|(O3) get TLM and ADJ models based on
the trajectory of high-res. NL model(O1) high-res. NL model(O2) departure:|(O3) low-res. NL model(O4) get TLM and ADJ models based on
the trajectory of low-res. NL model

Inner loop: Low-res. TLM & ADJ

Inner loop: Low-res. TLM & ADJ|(I1) Cost and its gradient||(I2) BGGS to update(I2) BGGS to update(O4) update analysis(O4) update analysis① iteration|① iteration(I1) Cost and its gradient② iteration

ECMWF (2015)② iteration

4DVARs in ECMWF & JMA

Use Low-res. NL that is consistent with TLM & ADJ

Training CourseDA Study w/ 40-variable Lorenz-96
Lorenz-96 model (Lorenz 1996)

Advection term

For j=1,…,N, Xj=Xj+N
Dissipation term Forcing termText Books① Training Description

pswd: ceres‣ https://kotsuki-lab.com/internal-pages/

Advanced Task 4

Advanced Task 4

4DVAR Implementation
StateBFGS--: background trajectory
One DA Window －: updated trajectory

Background

𝟎Observation
Analysis𝟎𝒊 𝒏

×: Smoother●: Filter××××𝒊 𝒏𝒊time

Point 1: the observation should not
be assimilated more than oncecorrect×incorrect

assimilatedin this cycle

𝒊 𝒏not assimilated
in this cycle

Point 2: Usually obs at t=0
is DA’ed by the last cycle𝒊

Point 3: Smoother is more accurate than filter
Analysis RMSEs of filter solutions should becompared w/ other filters (e.g. KF)

RMSEs in DA window
DA window(2 days)

Larger RMSE for later timein 4DVAR window

Smaller RMSE for earlier timein 4DVAR window

Smoother solutions
is more accurate
than filter solutions

DA window: 2 days, Bii=0.15, BFGS of Scipy

Sensitivity to B and window
when B is diagonal matrix.

 b110Bbii0bnn 
from KF’s ave

Analysis RMSE of filters are compared

KotsukiT. Saito

F. Kawasaki

Sensitivity to Obs. Network
Case 1: Num. Obs. = X
Case 2: Num. Obs. = X

Case 1

homogeneous
dense

Case 2

Iterations of 4DVAR (20 DAs)
By F. Kawasaki

terminated when

20 DAs

Scipy BFGS

Hand-write BFGS

Hand-write BFGS
terminated when

DA window: 2 days, Bii=0.15

Hand-write BFGS used Armijo (=0.1)

Visualization ofCost Function
F. Kawasaki

Experimental Setting‣ Focus on one DA by 4DVAR w/ Lorenz 96
NL Experiment

QUO: 4DVAR

Nonlinear Model

Nonlinear Model(N1) employ L96 model fcst
(N2) departure:(N3) get TLM and ADJ models(N1) employ L96 model fcst
(N2) departure:(N3) get TLM and ADJ models||

Linearized Model(L1) Cost and its gradient||(L2) BGGS to update(L2) BGGS to update(N4) update analysis(N4) update analysis
QUO: Quadratic unconstrained optimization

iteration

Linearized Model(L1) Cost and its gradient

iteration|

Estimating only two vars
‣ we estimate only x０ and x1 among 40 vars

‣ For x2-39, truth is used for computing J and grad J‣ While analysis increment is obtained for 40 vars,
only x０ and x1 are updated over the iteration‣ メモ:‣ X2=x39はtruthを使って、コスト＆勾配を計算‣ Grad(J)は40変数だけど、x0とx1だけを更新

Cost function & Iteration

AN of NL is closer to True than AN of QUO because(1) Updating TLM and ADJ around the trajectory(2) Nonlinear model is used for departure

QUONLComparison of solvers
Scipy BFGS (default)

Hand-write BFGS w/ Armijo•••))

Sensitivity to Armijo params

Larger reduces increments of an iteration••)Sensitivity to Armijo params

Larger increases increments of an iteration••)Tips forimplementations

For hand-write 4DVAR‣ (1) to test BFGS w/ a simple problem

Solve the minimization problem
from different initial states‣ (2) to test BFGS w/ 3DVAR‣ (3) to test BFGS w/ 4DVAR

(2) BFGS for 3DVAR (L96)
3DVAR can be solved by deterministically or iteratively by BFGS=+=BFGS
deterministicThank you for your attention!
Presented by Shunji Kotsuki(shunji.kotsuki@chiba-u.jp)
Further information is available athttps://kotsuki-lab.com/