Robin Hogan

1
Synergistic cloud, aerosol and precipitation
productsProgress so far in RATEC

• Robin Hogan
• Julien Delanoë
• Nicola Pounder
• University of Reading

2
Synergy products as defined in CASPER
3
Overview of our role in RATEC

• Radiative transfer for EarthCARE (RATEC)
• Lead contractor Howard Barker
• Started August 23, 2009
• 15 month project
• Our role is to start the development of synergy
  algorithms
• Also define the required imager and target
  classification products
• WP120 Review - first 3 months
• Products and Algorithms Requirements Document
  (PARD)
• Algorithm Development Plan (ADP)
• WP 220 Design - next 5 months
• Product Definition Document (PDD)
• WP 320 Prototype and test - final 7 months
• Algorithm Theoretical Basis Document (ATBD)

4
Overview of talk

• Laying the foundations for the algorithm...
• Summary of retrieval framework
• State variables (describing ice, liquid, rain and
  aerosol)
• Forward models (for radar, HSRL lidar, infrared
  solar radiometers)
• Cost function
• Minimization techniques
• Gauss-Newton Levenberg-Marquardt
• Gradient Descent (adjoint method)
• Ensemble methods
• Coding progress
• C library
• Scattering code
• Calculation and storage of error descriptors
• Solution error covariance matrix
• Averaging kernel
• Progress for specific atmospheric constituents
• Liquid cloud gradient constraint

5
Retrieval framework

• Ingredients developed before
• Not yet developed

6
Proposed state variables
State variable Representation with height / constraint A-priori
Ice clouds and snow
Visible extinction coefficient One variable per pixel with smoothness constraint None
Number conc. parameter Cubic spline basis functions with vertical correlation Temperature dependent
Lidar extinction-to-backscatter ratio Cubic spline basis functions 20 sr
Riming factor Likely a single value per profile 1
Liquid clouds
Liquid water content One variable per pixel but with gradient constraint None
Droplet number concentration One value per liquid layer Temperature dependent
Rain
Rain rate Cubic spline basis functions with flatness constraint None
Normalized number conc. Nw One value per profile Dependent on whether from melting ice or coallescence
Melting-layer thickness scaling factor One value per profile 1
Aerosols
Extinction coefficient One variable per pixel with smoothness constraint None
Lidar extinction-to-backscatter ratio One value per aerosol layer identified Climatological type depending on region
7
Forward model components

• From state vector x to forward modelled
  observations H(x)...

x
Ice snow
Liquid cloud
Rain
Aerosol
Lookup tables to obtain profiles of extinction,
scattering backscatter coefficients, asymmetry
factor
Lots of expensive matrix multiplications likely
to be the most expensive part of the entire
algorithm
Sum the contributions from each constituent
Radiative transfer models
8
Solution methods (1/3)

• We want to minimize this cost function
• Possibly using the gradient (a vector)
• and the second derivative (a matrix)
• 1. Gauss-Newton method (Rodgers p85)
• Advantage rapid convergence (instant convergence
  for linear problems)
• Another advantage get the error covariance of
  the solution for free
• Disadvantage need the Jacobian of every forward
  model can be expensive

9
Solution methods (2/3)

• 2. Gradient descent method (e.g. ECMWF data
  assimilation system)
• Advantage we dont need to calculate the
  Jacobian so forward model is cheaper!
• Quasi-Newton method to get search direction
  (conjugate gradient, BGFS etc)
• Disadvantage more iterations needed since we
  dont know curvature of J(x)
• Disadvantage dont get the error for free at
  the end
• Why dont we need the Jacobian H?
• The adjoint of a forward model takes as input
  the vector . and outputs the vector ??Jobs
  without needing to calculate the matrix H on the
  way
• Adjoint can be coded to be only 3 times slower
  than original forward model

Simple steepest descent requires many steps
Conjugate gradient method more intelligent
10
Solution methods (3/3)

• 3. Levenberg-Marquardt (hybrid of Gauss-Newton
  steepest descent)
• If cost function reduces, keep h small rapid
  convergence
• If cost function increases, increase h so that
  always go downhill
• More robust than Gauss-Newton, but still need the
  Jacobian
• 4. Ensemble method (e.g. Ensemble Kalman Filter)
• Create ensemble of trial solutions or
  particles in state space
• Spread of results allows Jacobian to be
  calculated automatically to create better set
  of particles repeat until convergence
• Applied to GPM retrievals by Mircea Grecu (NASA)
• Advantage no Jacobian required
• Estimate of retrieval error from final spread
• Disadvantage may need many ensemble
  members so probably slower than adjoint
  method
• 5. Simulated annealing (e.g. Donovan in ECSIM)
• Suitable for very non-linear problems
• Too slow for this application?

.
.
.
.
.
.
11
Radiative transfer forward models

• Computational cost can scale with number of
  points describing vertical profile N we can cope
  with an N2 dependence but not N3

Radar/lidar model Applications Speed Jacobian Adjoint
Single scattering bb exp(-2t) Radar lidar, no multiple scattering N N2 N
Platts approximation bb exp(-2ht) Lidar, ice only, crude multiple scattering N N2
Photon Variance-Covariance (PVC) method (Hogan 2006, 2008) Lidar, ice only, small-angle multiple scattering N or N2 N2 N
Time-Dependent Two-Stream (TDTS) method (Hogan and Battaglia 2008) Lidar radar, wide-angle multiple scattering N2 N3 N2
Depolarization capability for TDTS Lidar radar depol with multiple scattering N2 N2

• Lidar will use PVCTDTS, radar will use single
  scatteringTDTS
• Jacobian of TDTS is too expensive need to
  develop reduced-resolution Jacobian or
  alternatively adjoint code
• Also need depolarization forward model with
  multiple scattering

Radiometer model Applications Speed Jacobian Adjoint
RTTOV (used at ECMWF Met Office) Infrared and microwave radiances N N
Two-stream source function technique (e.g. Delanoe Hogan 2008) Infrared radiances N N2
RADIANT (from Graeme Stephens) Solar radiances N N2 N

• Infrared will probably use RTTOV, solar radiances
  will use RADIANT

12
Coding progress

• Coding has started on the underpinning libraries
• Use C because of key language features
• Matrix and Vector Expression Library
• Operator overloading and expression templates
  allows complex operations on matrices to be
  written on one line, with similar efficiency as
  Fortran-90
• Important feature is to allow a matrix object to
  behave as part of another matrix important when
  different parts of the Jacobian refer to
  different instruments and need to be passed to
  functions in bits
• Optimal Synergy Library
• Object orientation and polymorphism Radar and
  Lidar classes inherit from generic Observation
  class IceCloud and Aerosol inherit from generic
  Constituent class
• Enables code to be written to be completely
  flexible observations can be added and removed
  without needing to keep track of indices to
  matrices, so same code can be applied to many
  platforms
• Also written code to provide scattering libraries
  in NetCDF files

13
Storage of error descriptors
Measurements
Measurements

• If ice is present at N points and M variables are
  retrieved
• MN retrieved variables to report with MN errors
• Error covariance and averaging kernel matrices
  also needed
• A standard output from a variational algorithm
• But (MN)2/2 values in each matrix, much larger
  than the retrieved variables and their errors
• Matrix is a different size each ray!
• Need a suitable compression strategy...

Retrievals
Retrievals
14
Retrieval error covariance matrix (Extinction,
lidar ratio, N)
Top -gt Height -gt bottom
N
S
a
Top -gt Height -gt bottom
a
N
Top -gt Height -gt bottom
Top -gt Height -gt bottom
15
Retrieval error correlation matrix (Extinction,
lidar ratio, N)
Top -gt Height -gt bottom
N
S
a
Top -gt Height -gt bottom
a
N
Top -gt Height -gt bottom
Top -gt Height -gt bottom
16
For autocorrelation of extinction...

• 3 Gaussians requires 8N or fewer values to
  approximate the error covariance matrix
• Need to work on correlations between different
  variables

17
Averaging kernel

• The averaging kernel matrix expresses how the
  retrieval at a point depends on the truth at
  every other point
• The sum of each column expresses how much the
  retrieval at that point derives from the
  measurements (rather than the prior) 1 is good!
• The shape of each column expresses the effective
  vertical resolution of the retrieval more peaked
  is good!
• Easy to parameterize
• Store integral and peak
• Or store integral and the width for each column

Height (m)
Information content
18
Gradient constraint
A. Slingo, S. Nichols and J. Schmetz, Q. J. R.
Met. Soc. 1982

• We have a good constraint on the gradient of the
  state variables with height for
• LWC in stratocu (adiabatic profile, particularly
  near cloud base)
• Rain rate (fast falling so little variation with
  height expected)
• Not suitable for the usual a priori constraint
• Solution add an extra term to the cost function
  to penalize deviations from gradient c

19
Example in liquid clouds

• Using simulated observations
• Triangular cloud observed by a 1- or
  2-field-of-view lidar
• Retrieval uses Levenberg-Marquardt minimization
  with Hogan and Battaglia (2008) model for lidar
  multiple scattering

Optical depth30 Footprint100m Footprint600m
20
Ongoing work

• Implement various instrument forward models
• Not too difficult interfacing different codes
  and languages
• Write adjoint for slowest models
• Solver
• Test gradient descent algorithms (conjugate
  gradient, BGFS etc)
• Implementation of different atmospheric
  constituents
• Ice cloud we know how to do this
• Liquid cloud to be done in RATEC
• Aerosol and precipitation basic implementation
  in RATEC longer-term development needs
  collaboration with DAME, IRMA etc...
• Test datasets
• Several airborne comparisons for validation (TC4,
  POLARCAT etc)
• A-Train
• ECSIM (after RATEC)
• Two-instrument synergy products
• To be tested using same code removing one
  instrument (after RATEC)

21
(No Transcript)
22
Aerosol refractive index
355-nm lidar

• Values from Ellie Highwood
• Strong dependence on aerosol type
• Varies with hydration
• Can we use one continuous variable to represent
  different aerosol properties?
• Need to make simplifications if information from
  multiple wavelengths is to be combined

Solar MSI channels
Thermal MSI channels
23
Matrix and Vector Expression Library

• Intuitive interface and ability to do complex
  expressions
• Matrix M(3,3)
• Vector x(3) 1.0, 2.0, 3.0 // Simple
  initialization
• M N R exp(P) // invokes element
  multiplication
• y transpose(x) () inverse(R) // () invokes
  matrix multiplication
• x solve(A, b) // Linear algebra
• Important aspects not in any other library (to my
  knowledge)
• S.subset(M, 3, 30, 2, 20) // Matrix subsetting
• M(find(M lt 0.0)) 0.0 // Matlab-style indexing
• BitfieldMatrix b b.bit(4) 1 // Handling of
  classification data
• Efficient and robust implementation
• Expression templates to avoid temporary objects
  in long expressions
• Reference counting for allocated objects protects
  against memory leaks

24
Key classes

• State
• Methods
• solve()
• add_obs(Observation)
• Data
• List of Observations
• List of Constituents
• Vector x, y
• Matrix Jacobian etc.

• Constituent
• Methods
• calc_scat_props(wavelength)
• Data
• Vector x, x_prior

• Observation
• Methods
• calc_scat_props()
• virtual forward_model()
• Data
• Real wavelength
• Vector y, y_error
• Matrix Jacobian

25
Forward model part 1
State vector x
xice . . . . . xliq . . .
26
Forward model part 2
Radiative transfer
y'radar
Hradar
?yradar / ?xradar
Z . . .