VARIATIONAL DATA ASSIMILATION E'G' FOR ATMOSPHERIC MODELLING

1
VARIATIONAL DATA ASSIMILATION (E.G. FOR
ATMOSPHERIC MODELLING)

• Theory and basic computation
• Practical and physical aspects
• Examples and applications
• Related useful adjoint techniques
• Open questions

F. Bouttier (Météo-France), LNCC meeting,
Petropolis, Aug 2004
2
THEORY motivation
Compute the analysis values on a model grid
xa(A) Given observed values y1, y2, y3,. xb
background a priori estimate e.g. a previous
forecast xa(A)-xb(A) Si wi(yi-xbi) i.e.
problem is to compute wi Variational
analysis Define a real cost function J(w) of
vector wwi Such that wArg min J yields an
optimal analysis xa
Analysis problem
Y
y2
y1
y3
A
O
X
3
THEORY methodology

• Steps for variational analysis
• Define what to analyse control variable x
• Define an optimality criterion for xa (given the
  available data observation, background,
  physical constraints)
• Derive a cost-function J(x) which satisfies the
  optimality at its minimum
• Compute xaArg min J(x) as efficiently as
  possible
• Remark numerical efficiency will depend on
  choices made in all 4 steps need for compromises

4
THEORY a trivial example

• Compute Aub where u is the unknown vector
• Method 1 direct computation (Gauss elimination,
  etc.), expensive in a large space
• Method 2 define J(u) Au b 2 and
  minimize J
• Iterative solution is cheap if problem is
  well-posed and an approximation is good enough
  with just a few iterations of e.g. conjugate
  gradient minimizer
• Minimization of a well-posed quadratic
  cost-function is equivalent to solving a linear
  system
• J(u)uTAu2bTu
• when A symmetric positive definite, u minimizes J
  if and only if
• ?J0 i.e. Aub

5
Organization of sequential assimilation for
dynamical systems
time
obs
obs
3D analysis
xb
xb
xa
xa
etc
analysis
analysis
forecast
forecast
obs
obs
obs
obs
obs
obs
obs
4D analysis
xb
xa
etc
analysis
analysis
xaxb
6
THEORY optimal statistical analysis

• well suited to sequential problems combining
  imperfect data
• If xb is the a priori model state (background), B
  its estimation error covariance matrix,
• If y is the vector of observations, H the
  observation operator such that H(x) are
  model-simulated observations, R is the obs error
  covariance matrix (covs. of y-Hxtrue),
• If H can be linearized with respect to errors in
  x,
• Theorem the following equations define the same
  analysis xa
• xaxbBHT(HBHTR)-1(y-Hxb)
• xaArg min J with J(x)(x-xb)TB-1(x-xb)(y-Hx)TR-
  1(y-Hx)
• Jb(x)Jo(x)
• and xa is optimal in the sense that it has the
  smallest analysis error variance E(xa-xtrue)2

7
Impact of 2 terms on cost-function solution
8
Proof of statistical optimality

• Assuming xb-xtrue and y-H(xtrue) are unbiased and
  H can be linearized
• Criterion 1 xa is the most likely state given
  the knowledge of Gaussian pdfs (xb,B) and (y,R)
  (Bayesian maximum likelihood)
• Criterion 2 xa has the smallest expectation of
  quadratic error E(xa-xtrue)2 given the error
  covariances B and R
• Either criterion leads to the same optimum xa,
  the Best Linear Unbiased Estimator (BLUE)
• xaxbBHT(HBHTR)-1(y-Hxb)
• (B-1HR-1HT)-1HTR-1(y-Hxb)

9
THEORY variational optimal analysis

• Minimizing J(x)(x-xb)TB-1(x-xb)(y-Hx)TR-1(y-Hx)
• Is mathematically equivalent to solving for xa in
• xaxbBHT(HBHTR)-1(y-Hxb)
• With a small number of minimizer iterations, the
  variational form is approximate but much more
  efficient for large models (modern meteorological
  models have dim x gt108)
• Also works for a weakly non-linear H with some
  loss of optimality
• There exists a more general derivation of the
  equations without a need to distinguish xb and y
  (both are data sources, extra data can be added)

10
Glossary

• In J(x) (x-xb)TB-1(x-xb) (y-Hx)TR-1(y-Hx),
• background term Jb (x-xb)TB-1(x-xb)
• observation term Jo (y-Hx)TR-1(y-Hx)
• Iterative mimizers such as the conjugate gradient
  require the operator x ? J(x),?J(x) called the
  simulator
• In the gradient ?J2B-1(x-xb)2HTR-1(y-Hx) the
  term HT depends on the inner product and is
  called the adjoint of H
• The matrix of second derivatives J2B-12HTR-1H
  is called the Hessian of J
• preconditioner a change of variable ?Lx such
  that
• J(?)(J o L-1) is numerically easier to minimize
  than J
• penalty term an additionnal term Jc added to Jb
  and Jo in order to enforce constraints
  (initialization, coupling with another model,
  etc.)

11
Iterative minimization of cost-function
12
Basic computation

• On input, (xb,B,y,R,H) define the simulator
  software, xb the initial point of the
  minimization for the descent algorithm.
• Need to stop the minimization using a sensible
  criterion, e.g. ?J(xi) lt0.01 ?J (xini)
• The CPU-expensive parts are the computations of J
  and its gradient
• Minimization speed i.e. cost is mainly sensitive
  to the conditioning ellipticity of the Hessian.
  It must be improved by a good preconditioning
  e.g. using B yields a solution in about 100
  iterations.
• Minimization in a smaller space (e.g. PSAS in
  observation space) does not imply the analysis
  will be cheaper.

13
The importance of good preconditioningto reduce
cost-function ellipticity
14
2. Practical aspects 3DVar, 4DVar

• Variational analysis can solve the analysis
  problem at a given time its 3D-Var, to include
  in a data assimilation cycle (then xb, xa and y
  are valid at the same time)
• Or it can solve the whole data assimilation
  problem over a time interval (t1,t2) its
  4D-Var. Then,
• xb and xa are valid at t1
• y contains all the obs. between t1 and t2
• H contains a linearized forecast model from
  t1 to all observation times H ?i Hi?M1?i
• Mathematically 4D-Var is the same algorithm as
  3D-Var. It is much more expensive, the model
  adjoint software HT has to be coded, and it has
  interesting physical properties.

15
Principle of 4D-VAR assimilation
obs
Jo
previous forecast
analysis
Jo
obs
xb
corrected forecast
Jb
Jo
xa
obs
9h
12h
15h
Assimilation window
16
Variational assimilation 4D-Var

• 3D-Var is a static analysis algorithm uses obs
  at the time of analysis
• 4D-Var performs a complete data assimilation in a
  finite time window produces an optimal sequence
  of analysed states between t1 and t2
• 4D-Var is optimal if the model is perfect and can
  be linearized. Then it is mathematically
  equivalent to the Kalman filter, and much
  cheaper.
• In practice, 4D-Var is limited by model error and
  nonlinearities and its cost !
• 3DVar, 4DVar, KF are all limited by errors in B,
  R, H !
• 3D-FGAT is a cheap compromise between 3DVar and
  4DVar

17
2. Practical aspect observation processing

• Cost function is quadratic bad data may have
  very large weight
• Observation Quality Control reject unbelievable
  ones, thin too dense sets (i.e. with correlated
  errors) e.g. with variational QC
• Direct assimilation or preliminary retrieval
  (e.g. 1DVar on atmospheric columns) if difficult
  observation operator
• Accurate observation processing is essential for
  satellite data e.g. cloud-clearing of IR radiances

18
The weight of observations can be tuned
19
2. Practical aspect incremental technique

• A very useful trick to reduce the cost of 3DVar
  and 4DVar
• Use an approximate (low-resolution) cost-function
  in the vicinity of xb dxx-xb
• J(dx)dxTB-1x(y-Hdx-Hxb)TR-1(y-Hdx-Hxb)
• Similar to a Taylor expansion of H(xbdx)
• Makes J quadratic even if H is non-linear
• Ok if small scales are driven by large scales in
  dynamical system
• It will always be necessary for 4D-Var because
  4D-Var costs 100x more than the forecast model

20
Organization of 4D-Var incremental analysis
obs
obs
obs
obs
background forecast
xb
Compute y-HtMt(xb)
Minimize with y-HtMt(xb)-HtMt?x in Jo
xb
(3DVar FGAT assumes Midentity)
?x
updated forecast
Compute xa xb ?x
xa
21
Example of incremental technique storm surface
pressure corrections in 4D-Var
High-resolution model field 4D-Var field
22
2. Practical aspect Jb and the structure
functions

• Structure function shape of the increment xa-xb
  for a single observation
• Extremely important, determines the physical
  properties of the analysis
• For 1 obs, (HBHTR)-1(y-Hxb) is a scalar, HT is a
  column vector, and
• xa-xbBHT? where ? is a scalar
• the 3D structure of the increment is implied by B
• In 4DVar, HT includes the adjoint linearized
  model the structure functions are
  flow-dependent.
• Good design of B is the most important and
  difficult job when developing a 3D-Var and 4D-Var.

23
The impact of two versions of Jb on ECMWF
forecast performance
24
Typical 3DVar and 4DVar structure
functionscomputed using 1 geopotential
observation
3D-Var
24-h 4D-Var
25
2. How to build a good Jb term

• B is a symmetric, positive definite operator
• Its diagonal defines background error variances,
  which determine how closely the observations are
  fit
• The cross-covariances for each 3D field define
  the smoothing properties of the analysis, and the
  preferred increment structures ( most likely
  errors in xb)
• The inter-parameter covariances define physical
  balance properties geostrophy, convective
  balance
• Most of B is computed from forecast error
  statistics (Monte-Carlo methods and
  cross-validation)
• In nature, B is time- and location-dependent, the
  Kalman Filter aims to improve B
• The tropical atmosphere requires a specific B

26
Smoothing and geostrophy in Jb term (Northern
Hemisphere !)
27
Example in real 3DVar 1 pressure obs, structure
of wind correction
28
Interpolation of corrections between several
observations
? wind profiler station
29
3. Other applications of variational analysis

• 3D-Var and 4D-Var analyses of the atmosphere and
  the ocean for real-time forecasts (now in nearly
  all global weather centres !) and reanalysis
• 1D-Var analysis (retrievals) of atmospheric
  columns from satellite incomplete observations
• 2D-Var (4D-Var on a single column) analysis of
  soil moisture
• 3D-Var and 4D-Var analyses of atmospheric
  chemistry
• Applications in biology, fluid mechanics, etc.
  (systems with incomplete observations)

30
4. Related useful adjoint techniques

• variational thinking assuming the response of
  model and observations to errors is linear.
• Allows efficient linear algebra techniques
• Singular vectors (Lanczos algorithm) understand
  short-term error growth in forecasts, run
  ensemble predictions
• Adjoint sensitivity diagnose what was the cause
  of a forecast error (in analysis, observation,
  model tuning parameters)
• Optimality diagnostics check whether matrices B
  and R in J are consistent with errors in
  forecasts and observations
• Observation targeting on any single day, find
  where observations are most needed to make a good
  forecast
• Correction of a forecast by a human expert
• Efficient approximations to the Kalman filter,
  etc.

31
Conclusion open questions for scientists

• Variational analysis is a mathematically simple
  concept, successful because it is efficient and
  flexible for
• modern computers and remote-sensed data.
• Few fundamental scientific questions, but some
  interesting problems at the interface between
  maths and physics
• model non-linearities (e.g. clouds in the
  tropics)
• Choice of control variable and observation when
  errors are non-Gaussian (e.g. radars, humidity)
• Design and computation of B (ensemble techniques)
• Balance properties cloud structures, mountains
• Handling complex boundary conditions surface
  fluxes, lateral boundary coupling with a larger
  model
• (topics for visiting scientists in Météo-France)