Transcript Document

Retrieving consistent
profiles of clouds and
rain from instrument
synergy
Robin Hogan, Nicola Pounder,
Chris Westbrook University of Reading, UK
Julien Delanoë LATMOS, France
Alessandro Battaglia University of Leicester, UK
Overview
• Why a “unified” algorithm?
• Some new components
– New fast model for depolarization due to multiple scattering
– Automatic adjoints
• Ice, rain and melting-ice retrieval
– Testing on simulated profiles
• Demonstration on A-Train data
• Skill of global cloud forecasts
• Outlook
Why a “unified” algorithm?
• Combine all measurements available (radar, lidar, radiometers)
– Forms the observation vector y
• Retrieve cloud, precipitation and aerosol properties simultaneously
– Ensures integral measurements can be used when affected by more
than one species (e.g. radiances affected by ice and liquid clouds)
– Forms the state vector x
• Variational approach
– Minimizing a cost function J(x) allows for rigorous treatment of errors
• Aim to be completely flexible
– Applicable to ground-based, airborne and space-borne platforms
• Behaviour should tend towards existing two-instrument synergy algos
– Radar+lidar for ice clouds: Donovan et al. (2001), Delanoe & H (2008)
– CloudSat+MODIS for liquid clouds: Austin & Stephens (2001)
– Calipso+MODIS for aerosol: Kaufman et al. (2003)
– CloudSat surface return for rainfall: L’Ecuyer & Stephens (2002)
• This algorithm will provide one of the standard EarthCARE products
1. New ray of data: define state vector
Use classification to specify variables describing each species at each gate
Ice: extinction coefficient, N0’, lidar extinction-to-backscatter ratio
Liquid: extinction coefficient and number concentration
Rain: rain rate, drop diameter and melting ice
Aerosol: extinction coefficient, particle size and lidar ratio
Unified
retrieval
Ingredients developed
Implement previous work
Not yet developed
2. Convert state vector to radar-lidar resolution
Often the state vector will contain a low resolution description of the profile
3. Forward model
3a. Radar model
Including surface return
and multiple scattering
3b. Lidar model
Including HSRL channels
and multiple scattering
4. Compare to observations
Check for convergence
3c. Radiance model
Solar and IR channels
6. Iteration method
Derive a new state vector
Adjoint of full forward model
Quasi-Newton scheme
Not converged
Converged
7. Calculate retrieval error
Error covariances and averaging kernel
Proceed to next ray of data
Forward models
Observation
Radar reflectivity
factor
Radar reflectivity
factor in deep
convection
Radar Doppler velocity
HSRL lidar in ice and
aerosol
HSRL lidar in liquid
cloud
Lidar depolarization
Infrared radiances
Solar radiances
Model
Multiscatter: single scattering option
Speed
N
Status
OK
Multiscatter: single scattering plus TDTS MS
model (Hogan and Battaglia 2008)
N2
OK
Single scattering OK if no NUBF; fast MS model N2
with Doppler does not exist
Multiscatter: PVC model (Hogan 2008)
N
Not available
for MS
OK
Multiscatter: PVC plus TDTS models
N2
OK
Multiscatter: under development
Delanoe and Hogan (2008) two-stream source
function method
LIDORT (Robert Spurr)
N2
N
In progress
Needs to be
plugged in
Testing
N
• Multiscatter combines two fast multiple scattering models, PVC & TDTS
– Includes a Fortran-90 interface, adjoint model, HSRL capability...
– For lidar, much more accurate than Platt’s approximation with mu=0.7
– Can be used in retrievals and in instrument simulators
– Fast: One profile can cost the same as a single Monte Carlo photon!
• Freely available from http://www.met.rdg.ac.uk/clouds/multiscatter
• Can we model effect of multiple
scattering on depolarization?
• Potentially very useful information
on extinction (e.g. Sassen &
Petrilla 1986)
Battaglia et al. (2007)
Scattering regimes
• Regime 1: Single scattering
– Apparent backscatter b’ is easy to
calculate
– Zero depolarization from droplets
• Regime 2: Small-angle
multiple scattering
Footprint x
– Only for wavelength much less than
particle size, e.g. lidar & ice clouds
– Fast Photon Variance-Covariance
(PVC) model of Hogan (2008)
– Depolarization due to backscatter
slightly away from 180 degrees
• Regime 3: Wide-angle multiple
scattering
– Fast Time Dependent Two Stream
(TDTS) method of Hogan & Battaglia
– Depolarization increases with
number of scattering events
A typical Mie
phase function
for a distribution
of droplets
Fraction of cross-polar
rather than co-polar
scattered radiation
Forward
scattering is
unpolarized
The glory is
polarized
Compare new model to Monte Carlo using I3RC case
Wide-angle
scattering dominates
Small-angle
scattering
dominates
• Small-angle scattering: convolve cross-polar phase function with
modelled distribution of near-backscatter scattering angles
• Wide-angle scattering: assume that each scattering event randomizes
the polarization by a certain fraction = 0.6f + 0.85(1–f), where f is
the fraction of energy remaining in the field-of-view of the lidar
(coefficients derived by comparing to Alessandro’s Monte Carlo)
• New model appears to perform well for different fields of view
Unified retrieval: Forward model
• From state vector x to forward modelled observations H(x)...
Ice & snow
Liquid cloud
Rain
Aerosol
x
Gradient of cost function (vector)
xJ=HTR-1[y–H(x)]
Lookup tables to obtain profiles of extinction, scattering
& backscatter coefficients, asymmetry factor
Ice/radar
Ice/lidar
Ice/radiometer
Liquid/radar
Liquid/lidar
Liquid/radiometer
Rain/radar
Rain/lidar
Rain/radiometer
Aerosol/lidar
Aerosol/radiometer
Vector-matrix multiplications: around
the same cost as the original forward
operations
Sum the contributions from each constituent
Radar scattering
profile
Lidar scattering
profile
Radiative transfer models
Radar forward
modelled obs
Lidar forward
modelled obs
Radiometer
scattering profile
H(x)
Radiometer fwd
modelled obs
Adjoint of radar
model (vector)
Adjoint of lidar
model (vector)
Adjoint of radiometer
model
Adjoint of radiative transfer models
yJ=R-1[y–H(x)]
Minimization methods - in 1D
Quasi-Newton method (e.g. L-BFGS)
J
J/x
Gauss-Newton method
J
2J/x2
x1
x1
x2
x
x4 3
x
x6 5
xx8 7
J/x
x
x5
x4
x3
x2
x
• Requires the curvature 2J/x2
• Rolling a ball down a hill
– Smart choice of direction in many
– A matrix
dimensions aids convergence
– More expensive to calculate
• Requires the gradient J/x
• Fewer iterations to converge
– A vector (efficient to store)
– Assume J is quadratic and
– Efficient to calculate using
jump to the minimum
adjoint method
• Limited to smaller retrieval
• Used in data assimilation
problems
Coding up adjoints and Jacobian matrices is very time
consuming and error prone – is there a better way?
Automatic adjoints
• Algorithm y(x) in C++:
• Quite fiddly and error-prone to codeup dJ/dx given dJ/dy
adouble algorithm(const adouble x[2]) {
adouble
y = 4.0;
double
algorithm_AD(const
double x[2],
adouble
s = y_AD[1],
2.0*x[0] double
+ 3.0*x[1]*x[1];
double
x_AD[2]) {
y *= sin(s);
double
y = 4.0;
return
y;
double s = 2.0*x[0] + 3.0*x[1]*x[1];
}
y *= sin(s);
// Main
code
/* Adjoint
part: */
adept::Stack
stack;
double s_AD
= 0.0;
y = algorithm(x);
y_AD[0] += sin(s) * y_AD[0];
stack.compute_adjoint();
Do the hard work
s_AD += y * cos(s) *//
y_AD[0];
x_AD[0] += 3.0 * s_AD;
x_AD[1] += 6.0 * x[0] * s_AD;
s_AD = 0.0;
y_AD[0]
= 0.0;
This can be done automatically in C++
using
a special active type
return y;
“adouble” and overloading mathematical
operators and functions
}
double algorithm(const double x[2]) {
double y = 4.0;
double s = 2.0*x[0] + 3.0*x[1]*x[1];
y *= sin(s);
return y;
}
•
• Existing libraries CppAD and ADOL-C provide this capability but
they’re quite slow
• New library “Adept” uses expression templates in C++ to do this
much faster
• Freely available at http://www.met.rdg.ac.uk/clouds/adept
• Hogan (Mathematics & Computers in Simulation, in review)
Benchmark results
• Tested PVC and TDTS multiple scattering algorithms (here for PVC)
• Time relative to original code for profile with N=50 cloudy points:
Adjoint
Jacobian (50x350)
Hand-coded
3.0
New C++ library: Adept
3.5
20
ADOL-C
25
83
CppAD
29
352
• Adjoint calculation is:
– Only 5-20% slower than hand-coded adjoint
– 5-9 times faster than leading alternative libraries ADOL-C and CppAD
• Jacobian calculation is:
– 4-20 times faster than ADOL-C/CppAD for a matrix of size 50x350
• Now used for entire unified algorithm
• Sorry, it won’t work for Fortran until Fortran has template capability!
Ice retrieval: status
Same as Delanoe & Hogan (2008, 2010) ice-only retrieval
• State variables
– Extinction coefficient in geometric optics approximation
– Normalized number concentration N0*
– Lidar ratio
• Notable aspects
– Molecular signal beyond cloud is used when detectable
– Oblate spheroids with aspect ratio 0.6 for radar: good match to
aircraft data (Hogan et al. 2012)
– Doppler model using Heymsfield & Westbrook (2010) fall-speeds
• Remaining steps needed
– Add density as a retrieved variable to exploit Doppler in riming graupel
conditions
Liquid cloud retrieval: status
Discussed in Nicola Pounder’s talk
• State variables
– Liquid water content (actually ln LWC)
– Total number concentration (one number per liquid layer)
• Notable aspects
– Wide-angle lidar multiple scattering is included and provides useful
constraint on optical depth
– “One-sided gradient constraint” ensures retrieval near cloud base tends
towards the known adiabatic profile
• Remaining steps needed
– Forward model shortwave radiances and radar PIA
Rain and melting-layer: status
• State variables
– Rain rate
– Number concentration param (NL; currently fixed at a-priori value)
• Notable aspects
– Radar multiple-scattering included
– “Flatness constraint” on rain rate: extra terms in cost function penalize
deviations from a constant rain rate with height
– Different a-priori values for stratocu drizzle and rain from melting ice
– Melting layer radar attenuation uses Matrosov’s empirical relationship:
2-way attenuation [dB] = 2.2 x rain rate [mm h-1]
– Additional term in cost function penalizes difference in snow flux above
melting layer and rain rate below
• Remaining steps needed
– Use radar PIA information over the ocean
– Do we need a more complex melting-layer model?
Idealized nimbostratus retrieval
Constraint of constant flux
across melting layer and smooth
rain profile, but we have the
classic ill-constrained
attenuation retrieval problem
Idealized nimbostratus retrieval
Much better
behaviour with
flatness constraint
on rain rate
A-Train case
• Forward modelled radar and
lidar match observations
well, indicating good
convergence
• Can also simulate the
Doppler velocity that would
be observed by EarthCARE
– Currently omits multiple
scattering and air motion
effects on Doppler
Retrievals
• Ice cloud properties
retrieved similarly to
Delanoe and Hogan (2008,
2010) algorithm
• Water flux is
approximately conserved
across the melting layer
• Rain rate is relatively
constant with range
A mixed-phase case
• Observations
• Retrievals
Further work
• Forward models and observations
– Implement LIDORT solar radiance model (has adjoint/Jacobian)
– Implement Delanoe & Hogan infrared radiance code
– Implement multiple scattering model with depolarization (but are
measurements too noisy?)
– Incorporate radar path integrated attenuation
– Incorporate Doppler velocity
• Retrieved quantities
– Add “riming” factor
– Refine primitive aerosol retrieval
• Verification
– Consistency of different sources of information using A-Train data
– Aircraft data with in-situ sampling from NASA ER-2 and French
aircraft
– EarthCARE simulator (ECSIM) scenes using EarthCARE specification
Minimizing the cost function
1
1
T
T
1
J  y  H (x) R y  H (x)  x  a  B 1 x  a 
2
2
Gradient of cost function (a vector)
x J  HT R 1y  H (x)  B1 x  a
Gauss-Newton method

x i 1  x i   J
–
–
–
–
2
x

1
x J
Rapid convergence (instant for linear
problems)
Get solution error covariance “for
free” at the end
Levenberg-Marquardt is a small
modification to ensure convergence
Need the Jacobian matrix H of every
forward model: can be expensive for
larger problems as forward model may
need to be rerun with each element of
the state vector perturbed
and 2nd derivative (the Hessian matrix):
2x J  HT R 1H  B1
Gradient Descent methods
xi 1  xi  Ax J
– Fast adjoint method to calculate xJ
means don’t need to calculate Jacobian
– Disadvantage: more iterations needed
since we don’t know curvature of J(x)
– Quasi-Newton method to get the search
direction (e.g. L-BFGS used by ECMWF):
builds up an approximate inverse Hessian
A for improved convergence
– Scales well for large x
– Poorer estimate of the error at the end
Terminal fall-speed (m s-1)
Ice fall speeds
• Heymsfield & Westbrook
(2010) expression
predicts fall speed as a
function of particle mass,
maximum dimension and
“area ratio”
• Currently we assume
Brown and Francis (1995)
mass-size relationship, so
fall speed is a function of
size alone
• In convective clouds, intend to add a multiplication factor (or similar) to
allow denser particles (e.g. rimed aggregates, graupel and hail) to be
retrieved using the Doppler measurements
Simple melting-layer model
• Minimalist approach:
– 2-way radar
attenuation in dB is
2.2 times rain rate
(Matrosov 2008)
– No effect on other
variables
– Add term to cost
function penalizing
difference between
ice flux above and
rain flux below
melting layer
Matrosov (IEEE Trans. Geosci. Rem. Sens. 2008)
Model skill
• Use “DARDAR” CloudSatCALIPSO cloud mask
• How well is mean cloud
fraction modelled?
– Tend to underestimate
mid & low cloud fraction
• How good are models at
forecasting cloud at right
time? (SEDI skill score)
– Winter mid to upper
troposphere: excellent
– Tropical mid to upper
troposphere: fair
– Tropical and sub-tropical
boundary layer: virtually
no skill!
• Hogan, Stein, Garcon &
Delanoe (in prep)