Title: Robin Hogan
1Synergistic cloud, aerosol and precipitation
productsProgress so far in RATEC
- Robin Hogan
- Julien Delanoë
- Nicola Pounder
- University of Reading
2Synergy products as defined in CASPER
3Overview 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)
4Overview 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
5Retrieval framework
- Ingredients developed before
- Not yet developed
6Proposed 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
7Forward 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
8Solution 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
9Solution 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
10Solution 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?
.
.
.
.
.
.
11Radiative 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
12Coding 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
13Storage 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
14Retrieval 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
15Retrieval 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
16For autocorrelation of extinction...
- 3 Gaussians requires 8N or fewer values to
approximate the error covariance matrix - Need to work on correlations between different
variables
17Averaging 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
18Gradient 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
19Example 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
20Ongoing 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)
22Aerosol 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
23Matrix 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
24Key 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
25Forward model part 1
State vector x
xice . . . . . xliq . . .
26Forward model part 2
Radiative transfer
y'radar
Hradar
?yradar / ?xradar
Z . . .