Robin Hogan

1
Variational methods for retrieving cloud, rain
and hail properties combining radar, lidar and
radiometers

• Robin Hogan
• Julien Delanoe
• Department of Meteorology, University of Reading,
  UK

2
Outline

• 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)

3
Polarization 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

4
Polarization 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)
5
Polarization 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?

6
Using 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
7
Variational 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
8
How 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)

9
The 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
11
Finding 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
12
Chilbolton example 3-GHz radar25-m dish

• Observations
• Retrieval

Forward-model values at final iteration are
essentially least-squares fits to the
observations, but without instrument noise
13
A 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

14
Enforcing 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
15
Enforcing 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

16
Response to observational errors

• Nominal Zdr error of 0.2 dB Additional
  random error of 1 dB

17
What 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
18
Heavy rain andhail
Difficult case differential attenuation of 1 dB
and differential phase shift of 80º

• Observations
• Retrieval

19
How 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

20
Distribution 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

21
Enforcing 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

22
Summary

• 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.)
23
TheA-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

24
13.10 UTC June 18th
MODIS RGB composite
25
MODIS Infrared window
13.10 UTC June 18th
Lake district
Isle of Wight
France
England
Scotland
26
Met Office rain radar network
13.10 UTC June 18th
Lake district
Isle of Wight
France
England
Scotland
27
(No Transcript)
28
7 June 2006

• Calipso lidar
• CloudSat radar

5500 km
Eastern Russia
Japan
Sea of Japan
East China Sea
29
Motivation

• 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

30
Why 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
31
Formulation of variational scheme

• Observation vector State vector
• Elements may be missing
• Logarithms prevent unphysical negative
  values

32
Solution 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
33
Radar 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)
34
Lidar 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)
36
Radiance 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)

37
Ice 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

38
Variational 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

39
add 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)

40
add 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

41
add 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

42
add 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

43
Convergence

• 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

44
Radar-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

45
Radar 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

46
Radar, optical depth and IR radiances
Observations State variables Derived
variables
47
Ground-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
48
Retrieved 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
49
Conclusions 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

50
Sd

• sdf

An island of Indonesia
Banda Sea
51
Antarctic ice sheet
Southern Ocean