Robin Hogan

1
How to distinguish rain from hail using radarA
cunning, variational method

• Robin Hogan
• Last Minute Productions Inc.

2
Outline

• Increasingly in active remote sensing (radar and
  lidar), 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, it will be applied to polarization
  radar measurements of rain rate and hail
  intensity
• Met Office recently commissioned new polarization
  radar
• A variational retrieval is a very useful step
  towards assimilation of polarization data

3
Passive sensing

• Radiance at a particular wavelength has
  contributions from large range of heights
• A variational method is used to retrieve the
  temperature profile

4
Chilbolton 3GHz 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

5
Chilbolton 3GHz 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)
6
Chilbolton 3GHz 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?

7
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
8
Simple Zdr method
Observations

• Use Zdr at each gate to infer a in ZaR1.5
• Measurement noise feeds through to retrieval
• Noise much worse in operational radars

Noisy or Negative Zdr
9
Variational method

• Start with a first guess of coefficient a in
  ZaR1.5
• Z/R 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
10
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)

11
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)
12

• 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
13
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
14
First guess of a
First guess a 200 everywhere
Rainrate

• Use difference between the observations and
  forward model to predict new state vector (i.e.
  values of a), and iterate

15
Final iteration

• Zdr and fdp are well fitted by forward model at
  final iteration of minimization of cost function

Rainrate

• Retrieved coefficient a is forced to vary
  smoothly
• Prevents random noise in measurements feeding
  through into retrieval (which occurs in the
  simple Zdr method)

16
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
• Prevents random noise in the measurements feeding
  through into the retrieval

17
Enforcing smoothness

• In range cubic-spline basis functions
• Rather than state vector x containing a at
  every range gate, it is the amplitude of smaller
  number of basis functions
• Cubic splines ? solution is continuous in itself,
  its first and second derivatives
• Fewer elements in x ? more efficient!

Representing a 50-point function by 10 control
points

• In azimuth Two-pass Kalman smoother
• 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

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

20
Full scan from Chilbolton

• Observations
• Retrieval
• Note validation required!

Forward-model values at final iteration are
essentially least-squares fits to the
observations, but without instrument noise
21
Response to observational errors

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

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

• Observations
• Retrieval

24
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

25
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

26
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

27
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)
• Could be applied to new Met Office polarization
  radars
• Testing required higher frequency ? higher
  attenuation!

Hogan (2007, J. Appl. Meteorol. Climatology)
28
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
• Incorporate microwave radiances for deep
  precipitating clouds
• 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 radar problems, e.g. the radar
  refractivity method