Title: Robin Hogan
1Variational methods for retrieving cloud, rain
and hail properties combining radar, lidar and
radiometers
- Robin Hogan
- Julien Delanoe
- Department of Meteorology, University of Reading,
UK
2Outline
- Increasingly in active remote sensing, many
instruments are being deployed together, and
individual instruments may measure many variables - We want to retrieve an optimum estimate of the
state of the atmosphere that is consistent with
all the measurements - But most algorithms use at most only two
instruments/variables and dont take proper
account of instrumental errors - The variational approach (a.k.a. optimal
estimation theory) is standard in data
assimilation and passive sounding, but has only
recently been applied to radar retrieval problems - It is mathematically rigorous and takes full
account of errors - Straightforward to add extra constraints and
extra instruments - In this talk, two applications will be
demonstrated - Polarization radar retrieval of rain rate and
hail intensity - Retrieving cloud microphysical profiles from the
A-train of satellites (the CloudSat radar, the
Calipso lidar and the MODIS radiometer)
3Polarization radar Z
- We need to retrieve rain rate for accurate flood
forecasts - Conventional radar estimates rain-rate R from
radar reflectivity factor Z using ZaRb - Around a factor of 2 error in retrievals due to
variations in raindrop size and number
concentration - Attenuation through heavy rain must be corrected
for, but gate-by-gate methods are intrinsically
unstable - Hail contamination can lead to large
overestimates in rain rate
4Polarization radar Zdr
- Differential reflectivity Zdr is a measure of
drop shape, and hence drop size Zdr 10 log10
(ZH /ZV) - In principle allows rain rate to be retrieved to
25 - Can assist in correction for attenuation
- But
- Too noisy to use at each range-gate
- Needs to be accurately calibrated
- Degraded by hail
Drop
1 mm
ZV
3 mm
ZH
4.5 mm
ZDR 0 dB (ZH ZV)
- Drop shape is directly related to drop size
larger drops are less spherical - Hence the combination of Z and ZDR can provide
rain rate to 25.
ZDR 1.5 dB (ZH gt ZV)
ZDR 3 dB (ZH gtgt ZV)
5Polarization radar fdp
- Differential phase shift fdp is a propagation
effect caused by the difference in speed of the H
and V waves through oblate drops - Can use to estimate attenuation
- Calibration not required
- Low sensitivity to hail
- But
- Need high rain rate
- Low resolution information need to take
derivative but far too noisy to use at each
gate derivative can be negative! - How can we make the best use of the Zdr and fdp
information?
6Using Zdr and fdp for rain
- Useful at low and high R
- Differential attenuation allows accurate
attenuation correction but difficult to implement
- Need accurate calibration
- Too noisy at each gate
- Degraded by hail
Zdr
- Calibration not required
- Low sensitivity to hail
- Stable but inaccurate attenuation correction
- Need high R to use
- Must take derivative far too noisy at each gate
fdp
7Variational method
- Start with a first guess of coefficient a in
ZaR1.5 - Z/a implies a drop size use this in a forward
model to predict the observations of Zdr and fdp - Include all the relevant physics, such as
attenuation etc. - Compare observations with forward-model values,
and refine a by minimizing a cost function
Smoothness constraints
Observational errors are explicitly included, and
the solution is weighted accordingly
For a sensible solution at low rainrate, add an a
priori constraint on coefficient a
8How do we solve this?
- The best estimate of x minimizes a cost function
- At minimum of J, dJ/dx0, which leads to
- The least-squares solution is simply a weighted
average of m and b, weighting each by the inverse
of its error variance - Can also be written in terms of difference of m
and b from initial guess xi
- Generalize suppose I have two estimates of
variable x - m with error sm (from measurements)
- b with error sb (background or a priori
knowledge of the PDF of x)
9The Gauss-Newton method
- We often dont directly observe the variable we
want to retrieve, but instead some related
quantity y (e.g. we observe Zdr and fdp but not
a) so the cost function becomes - H(x) is the forward model predicting the
observations y from state x and may be complex
and non-analytic difficult to minimize J - Solution linearize forward model about a first
guess xi - The x corresponding to yH(x), is equivalent
to a direct measurement m - with error
y
Observation y
x
xi
xi1
xi2
(or m)
10- Substitute into prev. equation
- If it is straightforward to calculate ?y/?x then
iterate this formula to find the optimum x - If we have many observations and many variables
to retrieve then write this in matrix form - The matrices and vectors are defined by
Where the Hessian matrix is
State vector, a priori vector and observation
vector
Error covariance matrices of observations and
background
The Jacobian
11Finding the solution
New ray of data First guess of x
- In this problem, the observation vector y and
state vector x are
Forward model Predict measurements y and Jacobian
H from state vector x using forward model H(x)
Compare measurements to forward model Has the
solution converged? ?2 convergence test
No
Gauss-Newton iteration step Predict new state
vector xi1 xiA-1HTR-1y-H(xi)
B-1(b-xi) where the Hessian is AHTR-1HB-1
Yes
Calculate error in retrieval The solution error
covariance matrix is SA-1
Proceed to next ray
12Chilbolton example 3-GHz radar25-m dish
Forward-model values at final iteration are
essentially least-squares fits to the
observations, but without instrument noise
13A ray of data
- Zdr and fdp are well fitted by the forward model
at the final iteration of the minimization of the
cost function - The scheme also reports the error in the
retrieved values - Retrieved coefficient a is forced to vary
smoothly - Represented by cubic spline basis functions
14Enforcing smoothness 1
- Cubic-spline basis functions
- Let state vector x contain the amplitudes of a
set of basis functions - Cubic splines ensure that the solution is
continuous in itself and its first and second
derivatives - Fewer elements in x ? more efficient!
Forward model Convert state vector to high
resolution xhrWx Predict measurements y and
high-resolution Jacobian Hhr from xhr using
forward model H(xhr) Convert Jacobian to low
resolution HHhrW
Representing a 50-point function by 10 control
points
The weighting matrix
15Enforcing smoothness 2
- Background error covariance matrix
- To smooth beyond the range of individual basis
functions, recognise that errors in the a priori
estimate are correlated - Add off-diagonal elements to B assuming an
exponential decay of the correlations with range - The retrieved a now doesnt return immediately to
the a priori value in low rain rates - Kalman smoother in azimuth
- Each ray is retrieved separately, so how do we
ensure smoothness in azimuth as well? - Two-pass solution
- First pass use one ray as a constraint on the
retrieval at the next (a bit like an a priori) - Second pass repeat in the reverse direction,
constraining each ray both by the retrieval at
the previous ray, and by the first-pass retrieval
from the ray on the other side
16Response to observational errors
- Nominal Zdr error of 0.2 dB Additional
random error of 1 dB
17What if we use only Zdr or fdp ?
Retrieved a
Retrieval error
Zdr and fdp
- Very similar retrievals in moderate rain rates,
much more useful information obtained from Zdr
than fdp
Zdr only
Where observations provide no information,
retrieval tends to a priori value (and its error)
fdp only
fdp only useful where there is appreciable
gradient with range
18Heavy rain andhail
Difficult case differential attenuation of 1 dB
and differential phase shift of 80º
19How is hail retrieved?
- Hail is nearly spherical
- High Z but much lower Zdr than would get for rain
- Forward model cannot match both Zdr and fdp
- First pass of the algorithm
- Increase error on Zdr so that rain information
comes from fdp - Hail is where Zdrfwd-Zdr gt 1.5 dB and Z gt 35 dBZ
- Second pass of algorithm
- Use original Zdr error
- At each hail gate, retrieve the fraction of the
measured Z that is due to hail, as well as a. - Now the retrieval can match both Zdr and fdp
20Distribution of hail
Retrieved a
Retrieval error
Retrieved hail fraction
- Retrieved rain rate much lower in hail regions
high Z no longer attributed to rain - Can avoid false-alarm flood warnings
- Use Twomey method for smoothness of hail retrieval
21Enforcing smoothness 3
- Twomey matrix, for when we have no useful a
priori information - Add a term to the cost function to penalize
curvature in the solution ld2x/dr2 (where r is
range and l is a smoothing coefficient) - Implemented by adding Twomey matrix T to the
matrix equations
22Summary
- New scheme achieves a seamless transition between
the following separate algorithms - Drizzle. Zdr and fdp are both zero use a-priori
a coefficient - Light rain. Useful information in Zdr only
retrieve a smoothly varying a field (Illingworth
and Thompson 2005) - Heavy rain. Use fdp as well (e.g. Testud et al.
2000), but weight the Zdr and fdp information
according to their errors - Weak attenuation. Use fdp to estimate attenuation
(Holt 1988) - Strong attenuation. Use differential attenuation,
measured by negative Zdr at far end of ray
(Smyth and Illingworth 1998) - Hail occurrence. Identify by inconsistency
between Zdr and fdp measurements (Smyth et al.
1999) - Rain coexisting with hail. Estimate rain-rate in
hail regions from fdp alone (Sachidananda and
Zrnic 1987)
Hogan (2006, submitted to J. Appl. Meteorol.)
23TheA-train
- The CloudSat radar and the Calipso lidar were
launched on 28th April 2006 - They join Aqua, hosting the MODIS, CERES, AIRS
and AMSU radiometers - An opportunity to tackle questions concerning
role of clouds in climate - Need to combine all these observations to get an
optimum estimate of global cloud properties
2413.10 UTC June 18th
MODIS RGB composite
25MODIS Infrared window
13.10 UTC June 18th
Lake district
Isle of Wight
France
England
Scotland
26Met Office rain radar network
13.10 UTC June 18th
Lake district
Isle of Wight
France
England
Scotland
27(No Transcript)
287 June 2006
- Calipso lidar
- CloudSat radar
5500 km
Eastern Russia
Japan
Sea of Japan
East China Sea
29Motivation
- Why combine radar, lidar and radiometers?
- Radar Z?D6, lidar b?D2 so the combination
provides particle size - Radiances ensure that the retrieved profiles can
be used for radiative transfer studies - Some limitations of existing radar/lidar ice
retrieval schemes (Donovan et al. 2000, Tinel et
al. 2005, Mitrescu et al. 2005) - They only work in regions of cloud detected by
both radar and lidar - Noise in measurements results in noise in the
retrieved variables - Elorantas lidar multiple-scattering model is too
slow to take to greater than 3rd or 4th order
scattering - Other clouds in the profile are not included,
e.g. liquid water clouds - Difficult to make use of other measurements, e.g.
passive radiances - Difficult to also make use of lidar molecular
scattering beyond the cloud as an optical depth
constraint - Some methods need the unknown lidar ratio to be
specified - A unified variational scheme can solve all of
these problems
30Why not to invert the lidar separately
- Standard method assume a value for the
extinction-to-backscatter ratio, S, and use a
gate-by-gate correction - Problem for optical depth dgt2 is excessively
sensitive to choice of S - Exactly the same instability identified for radar
in 1954 - Better method (e.g. Donovan et al. 2000)
retrieve the S that is most consistent with the
radar and other constraints - For example, when combined with radar, it should
produce a profile of particle size or number
concentration that varies least with range
Implied optical depth is infinite
31Formulation of variational scheme
- Observation vector State vector
- Elements may be missing
- Logarithms prevent unphysical negative
values
32Solution method
- Find x that minimizes a cost function J of the
form - J deviation of x from a-priori
- deviation of observations from
forward model - curvature of extinction profile
New ray of data Locate cloud with radar
lidar Define elements of x First guess of x
Forward model Predict measurements y from state
vector x using forward model H(x) Also predict
the Jacobian H
Gauss-Newton iteration step Predict new state
vector xi1 xiA-1HTR-1y-H(xi)
-B-1(xi-xa)-Txi where the Hessian
is AHTR-1HB-1T
No
Has solution converged? ?2 convergence test
Yes
Calculate error in retrieval
Proceed to next ray
33Radar forward model and a priori
- Create lookup tables
- Gamma size distributions
- Choose mass-area-size relationships
- Mie theory for 94-GHz reflectivity
- Define normalized number concentration parameter
- The N0 that an exponential distribution would
have with same IWC and D0 as actual distribution
- Forward model predicts Z from extinction and N0
- Effective radius from lookup table
- N0 has strong T dependence
- Use Field et al. power-law as a-priori
- When no lidar signal, retrieval relaxes to one
based on Z and T (Liu and Illingworth 2000,
Hogan et al. 2006)
Field et al. (2005)
34Lidar forward model multiple scattering
- 90-m footprint of Calipso means that multiple
scattering is a problem - Elorantas (1998) model
- O (N m/m !) efficient for N points in profile and
m-order scattering - Too expensive to take to more than 3rd or 4th
order in retrieval (not enough) - New method treats third and higher orders
together - O (N 2) efficient
- As accurate as Eloranta when taken to 6th order
- 3-4 orders of magnitude faster for N 50 ( 0.1
ms)
Wide field-of-view forward scattered
photons may be returned
Narrow field-of-view forward scattered
photons escape
Ice cloud
Molecules
Liquid cloud
Aerosol
Hogan (2006, Applied Optics, in press). Code
www.met.rdg.ac.uk/clouds
35(No Transcript)
36Radiance forward model
- MODIS solar channels provide an estimate of
optical depth - Only very weakly dependent on vertical location
of cloud so we simply use the MODIS optical depth
product as a constraint - Only available in daylight
- MODIS, Calipso and SEVIRI each have 3 thermal
infrared channels in atmospheric window region - Radiance depends on vertical distribution of
microphysical properties - Single channel information on extinction near
cloud top - Pair of channels ice particle size information
near cloud top - Radiance model uses the 2-stream source function
method - Efficient yet sufficiently accurate method that
includes scattering - Provides important constraint for ice clouds
detected only by lidar - Ice single-scatter properties from Anthony
Barans aggregate model - Correlated-k-distribution for gaseous absorption
(from David Donovan and Seiji Kato)
37Ice cloud non-variational retrieval
Donovan et al. (2000)
Aircraft-simulated profiles with noise (from
Hogan et al. 2006)
Observations State variables Derived
variables
Retrieval is accurate but not perfectly stable
where lidar loses signal
- Donovan et al. (2000) algorithm can only be
applied where both lidar and radar have signal
38Variational radar/lidar retrieval
Observations State variables Derived
variables
Lidar noise matched by retrieval
Noise feeds through to other variables
- Noise in lidar backscatter feeds through to
retrieved extinction
39add smoothness constraint
Observations State variables Derived
variables
Retrieval reverts to a-priori N0
Extinction and IWC too low in radar-only region
- Smoothness constraint add a term to cost
function to penalize curvature in the solution
(J l Si d2ai/dz2)
40add a-priori error correlation
Observations State variables Derived
variables
Vertical correlation of error in N0
Extinction and IWC now more accurate
- Use B (the a priori error covariance matrix) to
smooth the N0 information in the vertical
41add visible optical depth constraint
Observations State variables Derived
variables
Slight refinement to extinction and IWC
- Integrated extinction now constrained by the
MODIS-derived visible optical depth
42add infrared radiances
Observations State variables Derived
variables
Poorer fit to Z at cloud top information here
now from radiances
- Better fit to IWC and re at cloud top
43Convergence
- The solution generally converges after two or
three iterations - When formulated in terms of ln(a), ln(b) rather
than a, b, the forward model is much more linear
so the minimum of the cost function is reached
rapidly
44Radar-only retrieval
Observations State variables Derived
variables
Use a priori as no other information on N0
Profile poor near cloud top no lidar for the
small crystals
- Retrieval is poorer if the lidar is not used
45Radar plus optical depth
Observations State variables Derived
variables
Optical depth constraint distributed evenly
through the cloud profile
- Note that often radar will not see all the way to
cloud top
46Radar, optical depth and IR radiances
Observations State variables Derived
variables
47Ground-based example
Observed 94-GHz radar reflectivity Observed
905-nm lidar backscatter Forward model radar
reflectivity Forward model lidar backscatter
Lidar fails to penetrate deep ice cloud
48Retrieved extinction coefficient a Retrieved
effective radius re Retrieved normalized number
conc. parameter N0 Error in retrieved extinction
Da
Radar only retrieval tends towards a-priori
Lower error in regions with both radar and lidar
49Conclusions and ongoing work
- Variational methods have been described for
retrieving cloud, rain and hail, from combined
active and passive sensors - Appropriate choice of state vector and smoothness
constraints ensures the retrievals are accurate
and efficient - Could provide the basis for cloud/rain data
assimilation - Ongoing work cloud
- Test radiance part of cloud retrieval using
geostationary-satellite radiances from
Meteosat/SEVIRI above ground-based radar lidar - Retrieve properties of liquid-water layers,
drizzle and aerosol - Apply to A-train data and validate using in-situ
underflights - Use to evaluate forecast/climate models
- Quantify radiative errors in representation of
different sorts of cloud - Ongoing work rain
- Validate the retrieved drop-size information,
e.g. using a distrometer - Apply to operational C-band (5.6 GHz) radars
more attenuation! - Apply to other problems, e.g. the radar
refractivity method
50Sd
An island of Indonesia
Banda Sea
51Antarctic ice sheet
Southern Ocean