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

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

4. Polarization 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)

5. Polarization 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?

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

7. How do we solve this? • 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) • The best estimate of x minimizes a cost function: • At minimum of J, dJ/dx=0, 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:

8. The Gauss-Newton method • We often don’t 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 y=H(x), is equivalent to a direct measurementm: …with error: y Observation y x xi+2 xi+1 xi (or m)

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

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

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

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

13. Enforcing smoothness 1 Forward model Convert state vector to high resolution:xhr=Wx Predict measurements y and high-resolution Jacobian Hhrfrom xhrusing forward modelH(xhr) Convert Jacobian to low resolution:H=HhrW • 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! Representing a 50-point function by 10 control points The weighting matrix

14. 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 doesn’t 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

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

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

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

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

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

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

21. 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) Hogan (2006, submitted to J. Appl. Meteorol.)

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

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

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

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

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

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

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

30. 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 Aerosol visible extinction coefficient profile • Observation vector • State vector • Elements may be missing • Logarithms prevent unphysical negative values

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

32. 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 (2006, Applied Optics, in press). Code: www.met.rdg.ac.uk/clouds

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

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

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

36. …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

37. …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

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

39. …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

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

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

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

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

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

45. 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 • 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

46. Sd • sdf An island of Indonesia Banda Sea

47. Antarctic ice sheet Southern Ocean