A variational data assimilation system for the global ocean

1
A variational data assimilation system for the
global ocean
Anthony Weaver CERFACS Toulouse,
France Acknowledgements N. Daget, E. Machu,
A. Piacentini, S. Ricci, P. Rogel (CERFACS) C.
Deltel, J. Vialard (LODYC, Paris) D. Anderson
and the ECMWF Seasonal Forecasting Group ENACT
project consortium (EC Framework 5)
2
Outline

• Scientific objectives and development strategy.
• General formulation and key characteristics of
  the variational system.
• Some results from tropical Pacific and global
  ocean applications.
• Summary and future directions.

3
Scientific objectives

• Develop a global ocean data assimilation system
  that can satisfy two purposes simultaneously
• Provide estimates of the ocean state over
  multi-annual to multi-decadal periods (currently
  up to 40 years ERA40).
• Provide ocean initial conditions for seasonal to
  multi-annual range forecasts.
• Much of this work has been coordinated through
  the EC-FP5 project ENACT (2002-2004).

4
Basic considerations in designing a practical
data assimilation system

• The assimilation method should have a solid
  theoretical foundation.
• It must be practical for large-dimensional
  systems involving GCMs
• State vector O(106) to O(107) elements.
• Non-differentiable parameterizations and
  algorithms.
• There should be a clear pathway to more advanced
  data assimilation systems.

5
CERFACS assimilation system development strategy
for OPA

• Develop a system based on variational data
  assimilation.
• Use an incremental approach (Courtier et al.
  1994, QJRMS).
• Provide a clear development pathway from
• 3D-Var
• 4D-Var short window, strong model constraint
• 4D-Var long window, weak model constraint

6
Development strategy cont.

• Why 3D-Var?
• An effective 3D-Var system provides a solid
  foundation for 4D-Var (they share most
  components!).
• 3D-Var is a simpler and cheaper alternative to
  4D-Var.
• 3D-Var provides a valuable reference for
  evaluating the cost benefits of 4D-Var.
• Some of the flow-dependent features implicit in
  4D-Var can be built into 3D-Var.
• 3D-Var requires significantly less maintenance
  and development than 4D-Var (the tangent-linear
  and adjoint of the forecast model are not
  needed).
• What can 4D-Var do that 3D-Var cant?
• 4D-Var can exploit tendency information in the
  observations.
• 4D-Var computes implicitly flow-dependent,
  time-evolving covariances within the assimilation
  window.

7
The OPA-VAR assimilation system

• OPA version 8.2 (Madec et al. 1999)
• Configurations available for assimilation
• Tropical Pacific (TDH) 1o x 0.5o at eq., 25
  levels (rigid lid)
• (Weaver et al. 2003, MWR Vialard et al. 2003,
  MWR Vossepoel et al. 2004, MWR Ricci et al.
  2004, MWR)
• Global (ORCA) 2o x 0.5o at eq., 31 levels (free
  surface)
• Variational assimilation environment
• 160 Fortran routines ( 43,000 lines) for the
  OPA tangent-linear and adjoint models, and
  associated validation routines.
• 190 Fortran routines ( 54,000 lines) for the
  rest (observation operators, covariance
  operators, minimization routine,).
• cf. 230 Fortran routines (76,000 lines) for
  the global OPA forecast model.

8
General formulation of the variational problem

• Let denote the vector of prognostic model
  state variables.
• Let denote the vector of analysis control
  variables where
• Find that minimizes
  where
• background term
• observation term
• where

9
Incremental formulation

• Let be an increment
  to the state
• Let be an increment
  to the control where
• Find that minimizes
  where

background term
quadratic obs. term
where
10
Choice of analysis control variables

• In the ocean model ( T, S, ?, u, v)
• As analysis variables we take ( T, Su,
  ?u, uu, vu)
• and assume these variables are mutually
  uncorrelated (so is block diagonal).
• The transformation is a
  balance constraint (Derber and Bouttier, 1999,
  Tellus)
• strong constraint if Su ?u uu vu 0
• weak constraint if Su ? ?u ? uu ? vu ? 0

unbalanced variables
11
Interpretation of the balance operator

• If is linear then
• defines the multivariate covariances in
• When dim( ) lt dim( ), has a null
  space.
• E.g., with applied as a strong
  constraint, the observations will project only
  onto the balanced modes.

12
Choice of balance operator

• We construct as a lower
  triangular matrix (and hence easily invertible)
  transformation using the following constraints
• Linearized local T-S relationships ?
  balanced S
• (Ricci et al. 2004, MWR)
• Dynamic height (baroclinic) ?
  balanced ?
• Geostrophy, ß-plane approx. near eq. ? balanced
  (u, v)
• We can interpret
• Su ? S(T) ?
  unbalanced S
• ?u barotropic component ?
  unbalanced ?
• (uu, vu) ageostrophic velocity ?
  unbalanced (u, v)

13
Multivariate 3D-Var covariances
Ex covariance relative to a SSH point at
(0o,144oW)
(surface)
(surface)
14
Choice of linear propagator

• involves integrating the nonlinear forward
  model
• from initial time to the
  observation times .
• involves integrating a linear forward
  model
• In 3D-Var (FGAT)
  persistence
• In 4D-Var
  approx. TL model
• where

15
Linear approximation in the tropical
Pacific (from Weaver et al. 2003, MWR)
Latitude
Latitude
16
TL approximation in the tropical Pacific (from
Weaver et al. 2003, MWR)
October start date
TIWs
17
Interpretation of the linear propagator

• The linear propagator defines how the background
  error covariances evolve within the assimilation
  window .
• E.g., for observations located only at time ,
  the effective background-error covariance matrix
  at is
• (cf. Extended Kalman filter)

18
Diagnosing implicit background temperature error
standard deviations ( ) in 4D-Var (Weaver et
al. 2003 MWR)
In 3D-Var
In 4D-Var (cf. EKF)
4D-Var (ti 30 days)
3D-Var
19
Impact of a single SSH observation in 4D-Var
SSH innovation 10 cm at (0o,160oW) at t 30
days
SSH analysis increment
Amplitude (cm)
Temperature analysis increment
Depth (m)
20
Preconditioning

• The minimization is preconditioned via a change
  of variables
• so that and
• where
• For a single observation, the minimization
  converges in a single iteration.

21
Specifying background error covariances general
remarks

• There is not enough information (and never will
  be) to determine all the elements of
  (typically gt O(1010)).
• must be approximated by a statistical
  model (e.g., prescribed covariance functions)
  with a limited number of tunable parameters.
• In 3D-Var/4D-Var, is implemented as an
  operator (a matrix-vector product).
• For the preconditioning transformation we require
  access to a square-root operator (and its
  adjoint ).
• Constructing an effective operator
  requires substantial development and tuning!

22
Modelling univariate background error covariances

• We solve a generalized diffusion equation (GDE)
  to perform the smoothing action of the
  square-root of the correlation operator (
  ).
• (Weaver and Courtier 2001, QJRMS Weaver and
  Ricci 2004, ECMWF Sem. Proc.)
• Simple parameterizations for the standard
  deviations of background error ( )
• (Balanced) T background vertical T-gradient
  dependent
• Unbalanced S background mixed-layer depth
  dependent
• Unbalanced SSH function of latitude
• Unbalanced (u,v) function of depth

23
Univariate correlation modelling using a
diffusion equation (Derber Rosati 1989 -
JPO Egbert et al. 1994 - JGR Weaver Courtier
2001 - QJRMS)

• A simple 1D example
• Consider with
  constant .
  
• on with
  as
• Integrate from and with
  as IC

24

• Solution
• This integral solution defines, after
  normalization, a correlation operator
• The kernel of is a Gaussian correlation
  function
• where is the length scale.
• Basic idea To compute the action of on a
  discrete grid we can iterate a diffusion
  operator.
• This is much cheaper than solving an integral
  equation directly.

25
Constructing a family of correlation functions
on the sphere using a GDE (Weaver Courtier
2001, QJRMS Weaver Ricci 2004 ECMWF Sem.
Procs.)
shape spectrum
L 500 km
Gaussian
Gaussian
26
Some remarks on numerical implementation

• The full correlation operator is formulated in
  grid-point space as a sequence of operators
• is the diffusion operator and is formulated
  in 3D as a product
• of a 2D (horizontal)
  and 1D (vertical) operator.
• is a diagonal matrix of volume elements,
  and appears in because of the
  self-adjointness of .
• The factor means iterations of
  the diffusion operator.

27
Some remarks on numerical implementation

• We can let where is
  a diffusion tensor that can be used to stretch
  and/or rotate the coordinates in the correlation
  model to account for anisotropic or
  flow-dependent structures.
• BCs are imposed directly within the discrete
  expression for using a land-ocean mask.
• contains normalization factors to ensure the
  variances of are equal to one.
• The diffusion approach to correlation modelling
  has many similarities to spline smoothing (Wahba
  1982) and recursive filtering (Purser et al. 2003
  - MWR).

28
GDE-generated correlation functions using
time-implicit scheme
Example T-T correlations at the equator
29
GDE-generated correlation functions
Example flow-dependent correlations (Weaver
Courtier 2001-QJRMS cf. Riishojgaard
1998-Tellus Daley Barker 2001-MWR)
Background isothermals
T-T correlations
Depth
15oN
15oS
15oS
15oN
30
Variational formulation main point

• The main scientific component of the algorithm is
  the transformation from control space to
  observation space in the Jo term

31
Incremental variational formulation

• And for incremental Var we need the linearized
  transformation (and its adjoint)

32
Impact of improved covariances on the mean zonal
velocity in the tropical Pacific
1993-96 climatology
eastward current bias
33
Impact of in situ T (GTSPP) data assim. on the
mean salinity state in the global model
Pacific
Atlantic
Indian
Equator
0
Control (no d.a.)
Depth (m)
Pacific
Atlantic
Indian
Equator
0
500
Depth (m)
Longitude
3D-Var univariate (T)
500
Longitude
34
Impact of in situ T (GTSPP) data assim. on the
mean salinity state in the global model
Pacific
Atlantic
Indian
Equator
0
Control (no d.a.)
Depth (m)
Pacific
Atlantic
Indian
Equator
0
500
Depth (m)
Longitude
3D-Var multivariate (T, S, u, v, SSH)
500
Longitude
35
Global reanalysis set-up and control

• Experimental set-up for ENACT
• Stream 1 1987 2001
• Stream 2 1962 2001
• Daily mean ERA-40 surface fluxes
• Weak 3D relaxation to Levitus T and S
• Strong relaxation to Reynolds SST (-200 W/m2/K)
• Control no data assimilation (streams 1 and 2)
• Getting a satisfactory control run was not
  straightforward!
• Post-correction to ERA-40 precipitation to remove
  a tropical bias.
• Stronger relaxation needed at high latitudes to
  avoid numerical instabilities.
• Daily correction to global mean E-P to remove sea
  level drift.

36
Global reanalysis experiments(completed or
currently running)

• 3D-Var (streams 1 and 2)
• In situ T data from ENACT QC data-set
• 10-day window
• Multivariate B (with balance)
• Incremental Analysis Updating (IAU) (Bloom et
  al. 1994, MWR)
• 4D-Var (stream 1)
• In situ T data from ENACT QC data-set
• 30-day window
• Univariate B (no balance)
• Instantaneous update

37
Cycling of 3D-Var and 4D-Var
using IAU
10-day window
observations
Analysis
30-day window
Background
Background trajectory
Analysed trajectory
38
ENACT QC historical in situ dataset (Met. Office)

• 200,000 300,000 in situ T observations / month
• 50,000 100,000 in situ S observations / month

Example of T data distribution on a 10 day window
Jan. 1987
Jan. 1995
39
Box regions for ENACT diagnostics
40
Assimilation diagnostics
1987-2001 global temperature statistics
Depth (m)
Mean (oC)
Standard deviation (oC)
41
Assimilation diagnostics
1987-2001 regional temperature statistics
NW extra-trop Pacific
NW extra-trop Pacific
Depth (m)
Mean (oC)
Standard deviation (oC)
42
Assimilation diagnostics
1987-2001 regional temperature statistics
Nino3
Nino3
Depth (m)
Mean (oC)
Standard deviation (oC)
43
Summary

• 3D-Var (FGAT) and 4D-Var incremental systems
  developed for a global version of OPA.
• Major coding and validation effort required.
• Clear development path towards more advanced
  systems.
• Substantial effort devoted to developing
  covariance models and balance operators.
• Balance constraints have a significant positive
  impact in 3D-Var.
• And a positive impact in 4D-Var with single
  observations (but has not yet been evaluated in
  real-data experiments).
• Production and assessment of global ocean
  reanalyses is ongoing (ENACT).
• Preliminary results indicate that the
  assimilation is correcting for a large model bias
  in the upper ocean.
• But assimilation is introducing a bias of its own
  below the thermocline.
• Further improvements to the assimilation system
  are needed

44
Future directions

• Background error modelling and estimation
• Observation error modelling and quality control
• Combined in situ T, S, altimeter and SST
  assimilation
• Model bias detection/correction
• Improving the computational efficiency of 4D-Var
• Ongoing reanalysis production and evaluation