Recent Developments in Data Assimilation at NCAR MM5WRF 3DVAR

1
Recent Developments in Data Assimilation at
NCAR(MM5/WRF 3DVAR)

• S. R. H. Rizvi
• National Center For Atmospheric Research
• NCAR/MMM, Bolder, CO-80307, USA
• Email rizvi_at_ucar.edu

2
Outline of talk

• Overview of MM5/WRF 3DVAR
• Algorithm
• Control variables
• Balance
• Statistical parameters
• Observations
• Computational efficiency
• Latest development results
• Future plan
• Conclusion

3
An overview of WRF/MM5 3DVAR
Namelist File
Xb
BE
Yo
3DVAR START
Setup Observations
Setup Background Errors
Read Namelist
Setup Background
Setup MPP
Compute Analysis
Calculate (O B)
Minimise Cost Function
Outer Loop
Output Analysis
Calculate Diagnostics
Tidy up
Diagnostic File
Xa
3DVAR END
4
3DVAR Algorithm

• Incremental 3DVAR approach
• Courtier et al. (1994) Veerse F. J.-N.
  Thepaut (1998)
• The cost function (J) is defined as
• J(X?) 1/2 X?T B1X? H?(X?) - dT
  (O F)1 H?(X?) d
• Cost 1/2 Background
  Observations
• X? Analysis increments (X - Xb )
• d Innovation vector (Yo- H(Xb) )
  
• Yo Observation vector
• H Forward (Non-linear) Observation
  Operator (FOO)
• H? Tangent linear operator of the forward
  operator, H
• B Background (previous forecast) errors
• O Observation (instrumental) errors
• F Representivity (observation operator)
  errors

5
3DVAR Algorithm Contd.

• In terms of control variable ( V , X? UV
  where B UUT ) the cost function (J) can
  written as
• J(V) 1/2 VVT H? (UV
  )- dT (O F) 1 H? (UV)- d
• Cost 1/2 Background
  Observatios
• Thus minimisation ( ?J/?v 0 ) of cost function
  (J) leads to,
• V - UT H?T (O F)1 (d H?UV) 0
  or
• AV R
  Analysis equation
• Where,
• A I UT H?T (O F)1 H? U
  and 
• R UT H?T (O F)1 d
• Analysis equation is solved for V and thus X? is
  determined

6
3DVAR Algorithm Contd.

• Practical implementation of 3DVAR requires
  simplifications
• Simplified error covariances.
• Linearized observation operators, balance
  equation.
• Thinning of observations.
• Suitable choice of analysis control variables
• etc.

7
Control variables

• In MM5/WRF there are three choices for the
  control variable
• cv_option 1
• U-component of wind
• V-component of wind
• Temperature
• Pressure
• Moisture variable as specific
  humidity or relative humidity
• cv_option 2
• Stream function (?)
• Velocity potential (?)
•    Unbalanced part of pressure (Pu)
• Moisture variable as specific humidity
  or relative humidity
• cv_option 3
• Stream function (?)
• Unbalanced part of velocity potential (?u)
• Unbalanced part of temperature (Tu)
• Log of surface pressure
• Pseudo relative humidity

8
Control variables Contd.

• Control variable (V) is defined as,
• X? UV , U Up
  Uv Uh
• Where, B E L ET UUT , U E L1/2
• E and L are the eigenvectors and eigenvalues of
  B.
• Uh Horizontal transform, using Recursive
  filters
• 
  Purser and Hyden (1998)
• Uv Vertical transform, using eigenvectors (E
  and L)
• Up Physical transform, using Physical/Dynamical
  laws

9
Recursive Filter
Right moving Bi a Bi-1 (1- a) Ai
(Input A output B Left moving Ci a Ci1
(1- a) Bi (Input B output C ) a
smoothing factor (0 lt a lt 1) Characteristic
Scalelength (R)
R, N and d are fixed to get a
10
Recursive Filter
Number of Passes
11
Recursive Filter
Response of Scalelength
12
Horizontal background errors

• Uh Isotropic/homogeneous recursive filter
  algorithm.

Correlation lengthscale
Single T obs (O-B1K, p500hPa)
13
Vertical background errors

• Uv Vertical EOFs transform

Eigenvectors of PSI
Single U obs (O-B1m/s, p200hPa)
14
Wind and Mass Balance
Linearized geostrophic, cyclostrophic balance
equation
Unbalanced pressure
C Regression coefficient
(Statistically determined)
15
Statistical parameters

• Observation error
• Background error
• Balance rgeression coefficients
• Characteristic lengthscale

16
Observations

• Conventional
• Upper air (TEMP, PIBAL, AIREP, ACARS,PROFILER).
• Surface (SYNOP, METAR, SHIP,BUOY)
• Remotely sensed retrievals
• Cloud-track winds (SATOBS).
• ATOVS thicknesses (SATEMs).
• Ground-based GPS TPW/ZTD.
• SSM/I oceanic surface wind speed and TPW.
• SSM/T1 temperature retrievals.
• SSM/T2 relative humidity retrievals.
• Scatterometer (Quikscat) oceanic surface winds.
• Radiances
• SSM/I brightness temperatures.

17
Statistical approximations

• Climatological background errors
• Estimated via tuned NMC-method
  statistics
• Simplified horizontal background error
  covariances represented by simple recursive
  filters
• Uncorrelated observation errors
• Neglect error correlations between analysis
  variables (streamfunction, potential,
  unbalanced pressure and humidity variable (q or
  RH)
• Approximate balance relationship
• Geostrophic, Cyclostrophic, Hydrostatic
  increments.

18
Computational efficiency

• 3DVAR MPP Domain Decompositions

1
Recursive Filters and FFTs
2
1
2
3
4
3
4
Minimisation
2
1
1
3
Obs. Operators
4
3
4
3
19
Computational efficiency Contd.
IBM-SP, Domain size 140x150x41
20
Computational Efficiency Contd.

• Data Compression Via Truncation Of Vertical EOFS
• Cost for 100x100x31 (CAA domain) 45km 3DVAR
• with conventional observations

Conclusion Halve 3DVAR cost with data, with
negligible lose of accuracy
21
3DVAR/MM5 AFWA Global Theaters
22
3DVAR/MM5 AFWA Tropical Theaters
23
Recent Developments

• Conjugate Gradient Minimisation
• Implementation of outer loop
• Background error computation
• Surface data assimilation
• Improved vertical interpolation
• Assimilation of Radar data

24
Conjugate Gradient (CG) Method

• Conjugate gradient method is an efficient way of
  solving simultaneous system of linear equations,
  
• A X B
• It is an iterative method and converges very fast
  if A is positive definite matrix
• Convergence is accelerated by selecting the
  search direction to make sure that the search is
  always made in a direction perpendicular to the
  direction already searched.
• Following Golub and Van Loan (1990) and Chandra
  (1978), CG algorithm have been coded in MM5/WRF
  3DVAR as an additional minimisation option.

25
CG Algorithm

• Equation to be solved (Analysis equation) AV
  R , where
• A I UTH?T (O F)1 H?U
  and
• R UTH?T (O F )1 d
• Set, r0 - R and P1 - r0
• where R UTH?T(O F)1d
• For k ? 1
• fk A Pk where A I UTH?T (O
  F)1 H?U
•   S ?rk-1 , rk-1? / ?Pk , fk?
• Vk Vk-1 S ? Pk
• rk rk-1 S ? fk
• Pk1 - rk ?rk , rk? / ?rk-1 , rk-1? ?
  Pk
• Iterate the above sequence till the desired
  convergence is achieved

26
CG Performance
27
CG Efficiency
28
Implementation of outer loop

• The analysis equation (AV B) is solved using
  double
• iteration loop as follows
• Set X Xb , so that V 0
•   Start of outer iteration
• X Xb U V
• Start of inner iteration
• Compute R UTH?T(O F)1d
• Solve A d R for d
• Update V V d
• End of inner loop
• End of outer loop

29
Outer loop Performance
30
Outer loop advantages

• Use of additional data
• Non-linearity of the forward observation operator
• Multiple background Quality control
• Efficient utilisation of PBL information for
  assimilation of surface data
• Effective utilisation of meso-scale data

31
Background errors

• The background errors are computed using
  NMC-method (Parrish and Derber,1992).
• This is the most popular method used at various
  operational NWP centers
• Generally one month forecast data is used for
  computing the various statistical parameters
  including the background errors, used in 3DVAR
• It is highly compute intensive and requires huge
  amount of computing resources.
• Currently MM5/WRF 3DVAR uses interpolated
  statistics generated from AVN forecast.
• The original statistics is at 101x181x21 with 210
  Km. resolution.

32
Global Domain For Calculation Of AVN
Background Errors101 x 181 x 21 at 210 Km.
resolution
33
New Background errors

• For understanding the issues involved with using
  the interpolated statistics new Background errors
  were generated for the following three regions
• Indian 75 x 101 x 23 at
  90 Km.
• AMPS 181 x 101 x 29 at 90
  Km.
• T4B 226 x 289 x 41 at
  15 Km.

34
Lengthscale (psi)
35
CPU requirements
36
Single-obs (Temperature) test (T4B)

• New BE

Old BE
37
Wind response (T4B)

• N
• E
• W
• BE

O L D BE




38
Single-obs (Temperature) test (AMPS)


New BE


Old BE
39
Wind response (AMPS)

• N
• E
• W
• BE

O L D BE
40
New approach for Surface data assimilation

• Why new approach?
• It was observed that too many surface observation
  reports, specially over the complex terrain, were
  getting rejected
• While using surface obs by reducing it to the
  lowest model sigma level some background error
  information also gets blended into obs

41
Surface data assimilation - Contd.

• Forward observation operator and its adjoint are
  developed to compute 10-m wind and 2-m
  temperature, moisture based on Monin-Obukhov
  similarity theory Cardinali et al. (1994),
  Courtier et al. (1998) and
• Guo et al. (2002)
• Suitable correction is applied to Surface
  Pressure corresponding to the difference between
  the actual and model terrain height
• Conclusion 50-60 additional surface data
  utilisation

42
Improved vertical interpolation
h, p, T
k1
Vertical Interpolation Old 1) If ho not
observed, derive ho from po . 2) Interpolate in
h. New 1) Interpolate in h or p (depends on
observation).
To, po
k
h, p, T
k-1
h, p, T
43
Radar Data Assimilation

• Additional control variables for assimilating
• Reflectivity Radial velocity
• - vertical velocity (w)
• - cloud rain water (qr)
• - cloud liquid water/ice (qc)
• - cloud water vapor (qv)
• Full Micro-Physics have been added as a part of
  FOO and its tangent linear code is developed and
  integrated with the minimisation scheme.

44
Radar data assimilation Contd.

• Observation operator (Sun and Crook, 1998)
• ri - distance between radar and the observation
• qr - rainwater

45
Radar data assimilation Contd.

• Richardson equation which is the combination of
  continuity, thermodynamic and hydrostatic
  relationship is used to diagnose vertical
  velocity (W)
• Linearized Adiabatic Richardson equation for W'

46
Future plan

• a) Community WRF 3DVAR
• Testing WRF/WRF 3DVAR cycling.
• Release new version (V2.0) in June 2004 and
  support.
• b) New observations
• IR/MW radiances (ATOVS, etc).
• Radar reflectivity.
• GPS refractivity.
• c) Ground temperature as new control variable
• d) NCAR real-time applications
• Antarctica Mesoscale Prediction System (AMPS).
• CONUS (impact of GPS data).
• e) Variational algorithm
• Synoptically-dependent background errors.
• First-guess at analysis time (FGAT), incremental
  4DVAR?

47
Conclusion

• For 3DVAR, CG method is ideally suited. The
  efficiency comes from faster convergence and
  memory saving.
• Implementation of Outer loop is very useful
  for effective data utilization and quality
  control point of view.
• In 3DVAR, it is not good idea to use interpolated
  statistical parameters from the statistics which
  is computed at much coarse resolution. However,
  there is no problem when it is used the other
  way.
• For computation of 3DVAR input statistics, it
  will always be good to have an initial estimate
  of the lengthscale parameter with limited
  dataset. It helps a lot in saving the
  computational cost.
• It will always be desirable to use flow dependent
  BE.

48
Acknowledgments

• Y.-H Kuo
• Dale M Barker
• Yong -R Guo
• Wei Huang
• Xiao Qingnong
• Mi-Seon Lee
• etc. at NCAR

49
Thank you !