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

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 (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 don’t 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. Active sensing Passive sensing No attenuation With attenuation Isolated weighting functions (or Jacobians) so don’t need to bother with variational methods? With attenuation (e.g. spaceborne lidar) weighting functions are broader: variational method required • 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 factorZ using Z=aRb • 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= 10log10(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 > ZV) ZDR = 3 dB (ZH >> ZV)

6. Chilbolton 3GHz radar: fdp phase shift • 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. Simple Zdr method Retrieval Noisy or no retrieval Rainrate Lookup table Observations • Use Zdr at each gate to infer a in Z=aR1.5 • Measurement noise feeds through to retrieval • Noise much worse in operational radars Noisy or Negative Zdr

8. Variational method • Start with a first guess of coefficient a in Z=aR1.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

9. 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 modelH(x) Compare measurements to forward model Has the solution converged? 2 convergence test No Gauss-Newton iteration step Predict new state vector: xi+1= xi+A-1{HTR-1[y-H(xi)] +B-1(b-xi)} where the Hessian is A=HTR-1H+B-1 Yes Calculate error in retrieval The solution error covariance matrix is S=A-1 Proceed to next ray

10. First guess of a • Use difference between the observations and forward model to predict new state vector (i.e. values of a), and iterate First guess: a =200 everywhere Rainrate

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

12. 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

13. 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

14. Response to observational errors Nominal Zdr error of ±0.2 dB Additional random error of ±1 dB

15. 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

16. Heavy rain andhail Difficult case: differential attenuation of 1 dB and differential phase shift of 80º • Observations • Retrieval

17. How is hail retrieved? • Hail is nearly spherical • High Z but much lower Zdrthan would get for rain • Forward model cannot match both Zdr andfdp • First pass of the algorithm • Increase error on Zdrso that rain information comes from fdp • Hail is where Zdrfwd-Zdr> 1.5 dB and Z > 35 dBZ • Second pass of algorithm • Use original Zdrerror • 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 andfdp

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

19. Summary • New scheme achieves a seamless transition between the following separate algorithms: • Drizzle.Zdr andfdp are both zero: use a-prioria 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 (2006, submitted to J. Appl. Meteorol.)

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

21. 13.10 UTC June 18th Lake district Isle of Wight France England Scotland MODIS RGB composite

22. MODIS Infrared window 13.10 UTC June 18th Lake district Isle of Wight France England Scotland

23. Met Office rain radar network 13.10 UTC June 18th Lake district Isle of Wight France England Scotland

25. 7 June 2006 Aerosol from China? Molecular scattering Non-drizzling stratocumulus Cirrus Mixed-phase altocumulus Drizzling stratocumulus Rain • Calipso lidar • CloudSat radar 5500 km Japan Eastern Russia East China Sea Sea of Japan

26. 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 • Eloranta’s 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

27. Why not 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 d>2 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

28. Formulation of variational scheme Ice visible extinction coefficient profile Attenuated lidar backscatter profile Ice normalized number conc. profile Radar reflectivity factor profile (on different grid) Extinction/backscatter ratio for ice Liquid water path and number conc. for each liquid layer Visible optical depth Infrared radiance Radiance difference Microwave radiances for precipitation? Aerosol visible extinction coefficient profile • Observation vector • State vector • Elements may be missing • Logarithms prevent unphysical negative values

29. 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)

30. Lidar forward model: multiple scattering Wide field-of-view: forward scattered photons may be returned Narrow field-of-view: forward scattered photons escape • 90-m footprint of Calipso means that multiple scattering is a problem • Eloranta’s (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) Ice cloud Molecules Liquid cloud Aerosol Hogan (Applied Optics, 2006). Code: www.met.rdg.ac.uk/clouds

31. 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 • Likely to be degraded by 3D cloud effects • MODIS, Calipso and SEVIRI each have 3 thermalinfrared 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 Baran’s aggregate model • Correlated-k-distribution for gaseous absorption (from David Donovan and Seiji Kato)

32. Ice cloud: non-variational retrieval Donovan et al. (2000) Aircraft-simulated profiles with noise (from Hogan et al. 2006) • Donovan et al. (2000) algorithm can only be applied where both lidar and radar have signal Observations State variables Derived variables Retrieval is accurate but not perfectly stable where lidar loses signal

33. Variational radar/lidar retrieval • Noise in lidar backscatter feeds through to retrieved extinction Observations State variables Derived variables Lidar noise matched by retrieval Noise feeds through to other variables

34. …add smoothness constraint • Smoothness constraint: add a term to cost function to penalize curvature in the solution (J’ = l Sid2ai/dz2) Observations State variables Derived variables Retrieval reverts to a-priori N0 Extinction and IWC too low in radar-only region

35. …add a-priori error correlation • Use B (the a priori error covariance matrix) to smooth the N0 information in the vertical Observations State variables Derived variables Vertical correlation of error in N0 Extinction and IWC now more accurate

36. …add visible optical depth constraint • Integrated extinction now constrained by the MODIS-derived visible optical depth Observations State variables Derived variables Slight refinement to extinction and IWC

37. …add infrared radiances • Better fit to IWC and re at cloud top Observations State variables Derived variables Poorer fit to Z at cloud top: information here now from radiances

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

39. Radar-only retrieval • Retrieval is poorer if the lidar is not used 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

40. Radar plus optical depth • Note that often radar will not see all the way to cloud top Observations State variables Derived variables Optical depth constraint distributed evenly through the cloud profile

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

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

43. Conclusions and ongoing work Lake district Isle of Wight France England Scotland • 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

44. Sd • sdf An island of Indonesia Banda Sea

45. Antarctic ice sheet Southern Ocean