Cloud and Precipitation Best Estimate ATBD
|EarthCARE Level 2 Documentation |
| |
|ACM-CAP |
|Cloud and Precipitation Best Estimate |
|Algorithm Theoretical Basis Document (ATBD) |
| |
|VARSY Project |
|Code: |L2b-ACM-BE-ATBD |
|Issue: |01-draft |
|Date: |04/04/2012 |
| |Name |Function |Signature |
|Prepared by |Robin Hogan |Project Scientist | |
|Reviewed by |Pavlos Kollias |Project Scientist | |
|Approved by |Pavlos Kollias |Project Manager | |
|Signatures and approvals on original |
| |
|[pic] |
| |
This page intentionally left blank
Document Information
|Contract Data |
|Contract Number: |4000104528/11/NL/CT |
|Contract Issuer: |ESA-ESTEC |
| | |
|Internal Distribution |
|Name |Unit |Copies |
|Pavlos Kollias |McGill University |1 |
|Aleksandra Tatarevic |McGill University |1 |
|Wanda Szyrmer |McGill University |1 |
|Internal Confidentiality Level |
|Unclassified |( |Restricted |( |Confidential |( |
|External Distribution |
|Name |Organisation |Copies |
|Tobias Wehr |ESA-ESTEC |1 |
|Michael Eisinger |ESA-ESTEC |1 |
|Dulce Lajas |ESA-ESTEC |1 |
|Robin Hoganogan |University of Reading |1 |
|Julien Delanoë |LATMOS |1 |
|Gerd-Jan van Zadelhof |KNMI |1 |
|David Donovan |KNMI |1 |
|Alessandro Battaglia |University of Leicester |1 |
|Archiving |
|Word Processor: |MS Word 2003 |
|File Name: |VARSY-ATBD-01-draft |
Document Status Log
|Issue |Change description |Date |Approved |
|RATEC 2.0 |Initial version: ATBD from RATEC project |01/11/2011 | |
|1.0 |First version for VARSY |10/4/2012 | |
| | | | |
| | | | |
Table of Contents
1. Purpose and Scope 10
2. Documents 11
2.1. Applicable Documents 11
2.2. Reference Documents 12
2.3. List of Abbreviations 12
2.4. Scientific references 13
3. Scientific Background of the algorithm 17
3.1. Ice cloud 17
3.2. Liquid cloud 18
3.3. Rain 18
3.4. Variational retrieval theory 19
3.4.1. Definition of the cost function 19
3.4.2. Minimization of the cost function 21
4. Justification for the selection of the algorithm 25
5. Mathematical algorithm Description 26
5.1. Input variables 26
5.1.1. AMB-MO Merged observations 26
5.1.1.1. Product description 26
5.1.1.2. Variable list 27
5.1.1.3. Notes 28
5.1.2. ACM-TC Target Classification 29
5.1.2.1. Product description 29
5.1.3. C-FMR Feature Mask and Radar Reflectivity factor 33
5.1.3.1. Product description 33
5.1.3.2. Variable list 33
5.1.4. C-CD Corrected Doppler 34
5.1.4.1. Product description 34
5.1.4.2. Variable list 34
5.1.5. X-MET Meteorological and surface data 34
5.1.5.1. Product description 34
5.1.5.2. Variable list 35
5.2. Configuration Parameters 36
5.3. Output variables 41
5.3.1. Product description 41
5.3.2. Variable List 42
5.4. Algorithm overview and flow chart 43
5.5. Algorithm definition 44
5.5.1. Class hierarchy 44
5.5.2. Global initialization 46
5.5.2.1. Global configuration and building object lists 46
5.5.2.2. Scattering look-up tables 46
5.5.2.3. Creating scattering look-up tables 49
5.5.3. Profile initialization 53
5.5.3.1. Contents of vectors and matrices describing retrieved quantities 53
5.5.3.2. Contents of vectors and matrices describing observed quantities 54
5.5.3.3. Populating the vectors and matrices 55
5.5.4. Minimizing the cost function 57
5.5.5. Calculation of the cost function 58
5.5.5.1. Different approaches to representing profiles of variables 60
5.5.5.2. Merging the scattering profiles 61
5.5.6. Calculation of the retrieval error 63
5.5.6.1. Calculation of the solution error covariance matrix 63
5.5.6.2. Transformation of error covariance matrix 63
5.5.6.3. One-dimensional error descriptors derived from the error covariance matrix 64
5.5.6.4. Error descriptors derived from the averaging kernel 64
5.5.7. Retrieved atmospheric constituents 65
5.5.7.1. Ice cloud 65
5.5.7.1.1. State variables 65
5.5.7.1.2. Microphysical and scattering assumptions 67
5.5.7.2. Liquid cloud 70
5.5.7.2.1. State variables 70
5.5.7.2.2. Microphysical and scattering assumptions, and additional constraints 71
5.5.7.3. Rain 72
5.5.7.3.1. State variables 72
5.5.7.3.2. Microphysical and scattering assumptions, and additional constraints 73
5.5.7.4. Melting layer 73
5.5.7.5. Aerosol 74
5.5.7.5.1. State variables 75
5.5.7.5.2. Microphysical and scattering assumptions 76
5.5.8. Radiative transfer forward models 77
5.5.8.1. CPR Reflectivity factor 78
5.5.8.2. CPR Doppler velocity 78
5.5.8.3. CPR Surface echo 79
5.5.8.4. ATLID HSRL Backscatter 79
5.5.8.5. ATLID Depolarization 80
5.5.8.6. MSI Thermal infrared radiances 80
5.5.8.7. MSI Solar radiances 81
6. Algorithm performance, sensitivity studies, limitations 82
6.1. Computational cost 82
6.2. Limitations 83
7. Validation status 85
7.1. Maturity 85
7.2. Liquid cloud example 89
8. Annex A: Implementation details and supporting software 91
8.1. Automatic differentiation 91
8.1.1. The problem 91
8.1.2. A solution 92
8.2. Scattering and terminal fall-speed calculations 93
8.3. Efficient storage of inverse of background error covariance matrix 94
List of Tables
Table 1: Applicable Documents 11
Table 2: Reference Documents 12
Table 3: List of abbreviations 13
List of Figures
Figure 1: Schematic of the operation of a quasi-Newton scheme for minimizing the cost function. 23
Figure 2: Top-level flowchart of the operations carried out by the algorithm. 43
Figure 3: Main class hierarchy. The arrows with solid heads indicate that a class contains pointers to one or more of the subclasses that the arrow points to. A class that points to another class with an open-headed arrow inherits from it. White boxes denote abstract classes that define a generic interface but cannot be instantiated, while the grey boxes denote concrete classes that can be instantiated. 45
Figure 4: Schematic of the operation of the forward model, which takes as input the state vector x and outputs the forward modelled observations yf=H(x). 62
Figure 5. (a) The temperature dependence of N0* for each size distribution within the large in situ database of Delanoë et al. (2005) (dots), superimposed by the mean N0* in 5[pic]C temperature ranges and various ranges of ice water content IWC (lines and symbols). (b) The same but for the variable N0’ = N0*/αv0.6. From Delanoë and Hogan (2008). 66
Figure 6. The colours represent the fall speed of ice particles in m s-1, as a function of their diameter (maximum dimension) and the fraction of the volume of the smallest sphere that encloses the particle that is ice, according to the model of Heymsfield and Westbrook (2010). The black line depicts the relationship between ice fraction and diameter of Brown and Francis (1995). 68
Figure 7. Aircraft measurements of the vertical profile of liquid water content and stratocumulus clouds (solid lines), together with the adiabatic profile calculated from the cloud-base temperature and pressure (Slingo et al. 1982). 71
Figure 8. Summary of the mean time per profile spent in each part of the algorithm, averaged over the processing of 1400 rays of simulated EarthCARE data (“test2” in the Algorithm Test Plan). In this scene the mean total time spent on each profile was 0.3 s. 83
Figure 9. Demonstration of the retrieval for a deep cloud observed by CloudSat and Calipso. (Top panel) Forward-modelled attenuated lidar backscatter coefficient at the final iteration of the algorithm, below which are the corresponding observed values. (Third panel) Forward modelled radar reflectivity factor, below which are the corresponding observed values. (Bottom panel) The forward modelled Doppler velocity (positive downwards) that would be measured in the absence of vertical wind and multiple scattering. The retrieved variables are shown in Figure 10. 86
Figure 10. Main retrieved atmospheric constituents for the case shown in Figure 9. 88
Figure 11. (Left first two panels) Calipso lidar backscatter and CloudSat radar reflectivity factor for a mixed-phase cloud with an ice cloud above. (Bottom-left panel) Forward modelled Doppler velocity. (Right panels) Retrieved variables. 89
Figure 12. Demonstration of retrieval of liquid water content profile from ATLID using simulated profiles and the Best Estimate algorithm. The top row shows the Mie and Rayleigh channels while the bottom row shows the true and retrieved profiles in each case. The left profile is of an adiabatic cloud retrieved using a smoothness constraint, the middle profile considers the same profile retrieved using a one-sided gradient constraint, while the right panel considers a cloud that is adiabatic only in the lower half again retrieved using the one-sided gradient constraint. 90
Purpose and Scope
This document describes the theoretical basis of the algorithm deriving the ACM-CAP data product.
The algorithm retrieves the properties of clouds, aerosols and precipitation from the combination of nadir-pointing cloud radar, lidar and radiometer. For each profile of data analyzed, the algorithm takes as input the ATLID-MSI “merged observations” product (AMB-MO) containing the lidar and imager observations interpolated or averaged to the same horizontal and vertical grid, the two CPR products C-FMR and C-CD containing the radar observations, the ACM-TC “target classification” product, which indicates what atmospheric constituents are present at each range gate (e.g. ice, liquid cloud, rain, aerosols and insects), and the X-MET product containing auxiliary meteorological, surface and other background field required to perform the retrieval. The target classification is used to define what variables to retrieve for each profile. Within this algorithm, a variational framework is employed to perform the retrieval, which enables errors and prior constraints to be incorporated rigorously. A key component is a sophisticated forward model, enabling a description of the atmosphere (specifically the detailed properties of clouds and precipitation) to be used to predict the observations; by then comparing these predictions with the actual observations, the description of the atmosphere can be refined. The retrieved properties of each constituent are then reported together with their estimated error and other diagnostic information.
The algorithm has been developed with the EarthCARE suite of instruments in mind: a Doppler cloud profiling radar (CPR), a high spectral resolution polarization atmospheric lidar (ATLID) and a multi-wavelength imager (MSI). The satellite will also carry a broad-band radiometer, which is not to be used for retrievals but would be available for independent evaluation). Despite being tailored for EarthCARE, the algorithm is designed to be very flexible, accepting arbitrary combinations of instrument types, and being applicable to both ground-based and space-borne platforms. This is essential for it to be tested on real data, since no current platform (e.g. airborne) has exactly the same specification as the EarthCARE instrument fit. It can be regarded as a “unification” of many previous algorithms, since when applied to one atmospheric constituent (e.g. ice clouds) observed with just one or two instruments (e.g. radar and lidar), it will tend towards similar behaviour to previous algorithms reported in the literature.
Documents
1 Applicable Documents
Table 1: Applicable Documents
|Reference |Code |Title |Issue |
|[SOW] |EC-SW-ESA-SY-0310 |Statement of Work: VARSY - 1-Dimensional VARiational Retrieval of|1.0 |
| | |SYnergistic EarthCARE Products | |
|[CC] |Appendix 2 to |Draft Contract (attachment to SOW) |1.0 |
| |AO/1-6823/11/NL/CT | | |
|[AD 1] |EC-SW-ESA-SY-0152 |EarthCARE Level 2 Processor Development |1.0 |
| | |General Requirements Baseline | |
|[AD 2] |EC.ICD.ASD.SY.00004 |EarthCARE Product Definitions. Vol. 0: | |
| | |Introduction | |
|[AD 3] |EC.ICD.ASD.SY.00005 |EarthCARE Product Definitions. Vol. 1: |1.0 |
| | |Common Product Definitions | |
|[AD 4] |EC.ICD.ASD.ATL.00021 |EarthCARE Product Definitions. Vol. 2b: |1.0 |
| | |ATLID level 1 | |
|[AD 5] |EC.ICD.ASD.BBR.00022 |EarthCARE Product Definitions. Vol. 3b: |1.0 |
| | |BBR level 1 | |
|[AD 6] |EC.ICD.ASD.MSI.00023 |EarthCARE Product Definitions. Vol. 4b: |1.0 |
| | |MSI level 1 | |
|[AD 7] |ECSIM-DMS-TEC-ICD01-R |ECSIM Simulator Interface Control Document | |
|[AD 8] |PE-TN-ESA-GS-0001 |Ground Segment: File Format Standard |1.0 |
|[AD 9] |EC-TN-ESA-GS-0218 |Tailoring of the Earth Explorer File Format Standard for the |2.0 |
| | |EarthCARE Ground Segment | |
2 Reference Documents
Table 2: Reference Documents
|Reference |Code |Title |Issue |
|[RD1] |ECSIM-DMS-TEC-SUM-01-R |ECSIM System User Manual | |
|[RD2] |ECSIM-KNMI-MAD01-R |ECSIM Model and Algorithms Document | |
|[RD3] |EE-MA-DMS-GS-0001 |Earth Explorer Mission CFI Software: | |
| | |General Software User Manual | |
|[RD4] |EOP-SM/1567/TW |EarthCARE Mission Requirements Document | |
|[CASPER-PARD] |CASPER-DMS-PARD-001 |CASPER Products and Algorithms Requirement Document (PARD) |2.0 |
|[CASPER-ATBD] |CASPER-DMS-ATBD-URSW1-10 |Algorithm Theoretical Basis Document (ATBD) for ACM-Ice-Reading |1.0 |
| | |met.reading.ac.uk/clouds/other_reports.html | |
|[RATEC-ATBD] |RATEC-ATBD-READING-2.0 |Algorithm Theoretical Basis Document (ATBD) for Cloud, Aerosol and |2.0 |
| | |Precipitation Best Estimate | |
|[RATEC-FR] |RATEC-FR-READING-1 |RATEC Final Report |1.0, April |
| | | |2011 |
|[VARSY-PARD] |VARSY-UR-PARD-1.0 |VARSY Products and Algorithms Requirement Document – Reading contribution |1.0 |
| | |(PARD) | |
|[VARSY-PDD] |VARSY-UR-PDD-1.0 |VARSY Products Definition Document – Reading contribution (PARD) |1.0 |
3 List of Abbreviations
Table 3: List of abbreviations
|Abbreviation |Name |
|1D-VAR RS |1-dimensional variational retrieval scheme |
|ATLID |Atmospheric Lidar (The EarthCARE lidar) |
|CASPER |Cloud and Aerosol Synergetic Products from EarthCARE retrievals |
|CPR |Cloud Precipitation Radar (The EarthCARE radar) |
|EarthCARE |The Earth Clouds, Aerosols and Radiation Explorer |
|ECSIM |EarthCARE Simulator |
|HSRL |High-Spectral Resolution Lidar |
|MSI |Multi-spectral Imager (The EarthCARE imager ) |
4 Scientific references
Austin, R. T., and G. L. Stephens, 2001: Retrieval of stratus cloud microphysical parameters using millimeter-wave radar and visible optical depth in preparation for CloudSat - 1. Algorithm formulation. J. Geophys. Res., 106, 28233–28242.
Baran, A. J., 2004: On the scattering and absorption properties of cirrus cloud. J. Quant. Spectroscopy Rad. Transfer, 89, 17–36.
Battaglia, A., S. Tanelli, S. Kobayashi, D. Zrnic, R. J. Hogan and C. Simmer, 2010: Multiple-scattering in radar systems: a review. J. Quant. Spectroscopy, 111, 917-947.
Battaglia, A., T. Augustynek, S. Tanelli and P. Kollias, 2011: Multiple scattering identification in spaceborne W-band radar measurements of deep convective cores. Submitted to J. Geophys. Res.
Beard, K. V., 1976: Terminal velocity and shape of cloud and precipitation drops aloft. J. Atmos. Sci. 33, 851–864.
Brandes, E. A., G. Zhang and J. Vivekanandan, 2002: Experiments in rainfall estimation with a polarimetric radar in a subtropical environment. J. Appl. Meteorol., 41, 674-685.
Brown, P.R., and P.N. Francis, 1995: Improved Measurements of the Ice Water Content in Cirrus Using a Total-Water Probe. J. Atmos. Oceanic Technol., 12, 410–414.
Delanoë, J., and R. J. Hogan, 2008: A variational method for retrieving ice cloud properties from combined radar, lidar, and infrared radiometer. J. Geophys. Res., 113, D07204, doi:10.1029/2007JD009000.
Delanoë, J. M. E., and R. J. Hogan, 2008b: Algorithm Theoretical Basis Document for ACM-Ice-Reading (variational radar-lidar-radiometer ice cloud retrieval). Produced for the European Space Agency as part of the project "CASPER" (Cloud and Aerosol Synergetic Products from EarthCARE Retrievals). Available from:
Delanoë, J., and R. J. Hogan, 2010: Combined CloudSat-CALIPSO-MODIS retrievals of the properties of ice clouds. J. Geophys. Res., D00H29, doi:10.1029/2009JD012346.
Delanoe, J., A. Protat, J. Testud, D. Bouniol, A. J. Heymsfield, A. Bansemer, P. R. A. Brown and R.M. Forbes, 2005: Statistical properties of the normalized ice particle size distribution. J. Geophys. Res., 110, D10201, doi:10.1029/2004JD005405.
Donovan, D. P., and A. C. A. P. van Lammeren, 2001: Cloud effective particle size and water content profile retrievals using combined lidar and radar observations - 1. Theory and examples. J. Geophys. Res., 106(D21), 27425-27448.
Fabry, F., and I. Zawadzki, 1995: Long-term radar observations of the melting layer of precipitation and their interpretation. J. Atmos. Sci., 52, 838-851.
Field, P. R., R. J. Hogan, P. R. A. Brown, A. J. Illingworth, T. W. Choularton and R. J. Cotton, 2005: Parametrization of ice particle size distributions for mid-latitude stratiform cloud. Q. J. R. Meteorol. Soc., 131, 1997-2017.
Fisher, M., and P. Courtier, 1995: Estimating the covariance matrices of analysis and forecast error in variational data assimilation. ECMWF Tech. Memo., 220, 28 pp.
Fox, N. I., and A. J. Illingworth, 1997: The potential of a spaceborne cloud radar for the detection of stratocumulus clouds. J. Appl. Meteorol, 36, 676-687.
Gilbert, J. C., and C. Lemaréchal, 1989: Some numerical experiments with variable-storage quasi-Newton algorithms. Mathematical Programming, 45, 407-435.
Giering, R., and T. Kaminski, 1998: Recipes for adjoint code construction. ACM Trans. on Mathematical Software (TOMS), 24, 437-474.
Grecu, M., and W. S. Olsen, 2006: Bayesian estimation of precipitation from satellite passive microwave observations using combined radar-radiometer retrievals. J. Appl. Meteorol. Climatol., 45, 416-433
Hansen, P. C., 1987: Rank-deficient and discrete ill-posed problems: numerical aspects of linear inversion. SIAM Philadelphia.
Hawkness-Smith, 2010: A novel retrieval of liquid water path and an evaluation of the representation of drizzle in numerical models. PhD thesis, University of Reading.
Hesse, M., P. Koepke and I. Schult, 1998: Optical properties of aerosols and clouds: the software package OPAC. Bull. Am. Meteorol. Soc., 79, 831-844.
Heymsfield, A. J., and C. D. Westbrook, 2010: Advancements in the estimation of ice particle fall speeds using laboratory and field measurements. J. Atmos. Sci., 67, 2469-2482.
Hogan, R. J., 2006: Fast approximate calculation of multiply scattered lidar returns. Appl. Opt., 45, 5984-5992.
Hogan, R. J., 2007: A variational scheme for retrieving rainfall rate and hail intensity from polarization radar. J. Appl. Meteorol. Climatology, 46, 1544-1564.
Hogan, R. J., 2008: Fast lidar and radar multiple-scattering models – 1. Small-angle scattering using the photon variance-covariance method. J. Atmos. Sci., 65, 3621-3635.
Hogan, R. J., 2013: Fast reverse-mode automatic differentiation using expression templates in C++. Submitted to ACM Trans. Mathematical Softw.
Hogan, R. J., and A. Battaglia, 2008: Fast lidar and radar multiple-scattering models – 2. Wide-angle scattering using the time-dependent two-stream approximation. J. Atmos. Sci., 65, 3636-3651.
Hogan, R. J., and C. D. Westbrook, 2013: Equation for the microwave backscatter cross section of aggregate snowflakes using the Self-Similar Rayleigh-Gans approximation. To be submitted to J. Atmos. Sci.
Hogan, R. J., C. Jakob and A. J. Illingworth, 2001: Comparison of ECMWF winter-season cloud fraction with radar derived values. J. Appl. Meteorol., 40, 513-525.
Hogan, R. J., M. P. Mittermaier and A. J. Illingworth, 2006a: The retrieval of ice water content from radar reflectivity factor and temperature and its use in the evaluation of a mesoscale model. J. Appl. Meteorol. Climatology, 45, 301-317.
Hogan, R. J., L. Tian, P. R. A. Brown, C. D. Westbrook, A. J. Heymsfield and J. D. Eastment, 2012: Radar scattering from ice aggregates using the horizontally aligned oblate spheroid approximation. J. Appl. Meteorol. Climatology., 51, 655-671.
Illingworth, A. J., and T. M. Blackman, 2002: The need to represent raindrop size spectra as normalized gamma distributions for the interpretation of polarization radar observations. J. Appl. Meteorol., 41, 286-297.
Josset, D., J. Pelon, A. Protat and C. Flamant, 2008: New approach to determine aerosol optical depth from combined CALIPSO and Cloud-Sat ocean surface echoes. Geophys. Res. Lett., 35, L10805, doi:10.1029/2008GL033442.
Kudryavtsev, V., 2003: A semiempirical model of the normalized radar cross-section of the sea surface. J. Geophys. Res., 108, doi:10.1029/2001JC001003.
L’Ecuyer, T.S., and G.L. Stephens (2002), An estimation-based precipitation retrieval algorithm for attenuating radars, J. Appl. Meteorol., 41, 272-285.
Li, L., G. M. Heymsfield, L. Tian and P. E. Racette, 2005: Backscattering measurements of ocean surface backscattering using an airborne 94-GHz cloud radar - implication for calibration of airborne and spaceborne W-band radars. J. Atmos. Oceanic Technol, 22, 1033-1045.
Liebe, H., 1985: An updated model for millimetre-wave propagation in moist air. Radio Sci., 20, 1069-1089.
Liu, Y., and J. Hallett, 1997: The ‘1/3’ power law between effective radius and liquid-water content. Q. J. R. Meteorol. Soc., 123, 1789-1795.
Liu, D. C., and J. Nocedal, 1989: On the limited memory method for large scale optimization. Mathematical Programming B, 45, 503-528.
Martin, G. M., D. W. Johnson and A. Spice, 1994: The measurement and parameterization of effective radius of droplets in warm stratocumulus clouds. J. Atmos. Sci., 51, 1823-1842.
Matrosov, S. 2007; Potential for attenuation-based estimates of rainfall rate from CloudSat, Geophys. Res. Lett., 34, L05817, doi: 10.1029/2006GL029161.
Matrosov, S. Y., 2008: Assessment of radar signal attenuation caused by the melting hydrometeor layer. IEEE Trans. Geosci. Rem. Sens., 46, 1039-1047.
Miles, N. L., J. Verlinde and E. E. Clothiaux, 2000: Cloud droplet size distributions in low-level stratiform clouds. J. Atmos. Sci., 57, 295–311.
Nocedal, J., 1980: Updating quasi-Newton matrices with limited storage. Mathematics of Computation, 35, 773-782.
O'Connor, E. J., A. J. Illingworth and R. J. Hogan, 2004: A technique for auto-calibration of cloud lidar. J. Atmos. Oceanic Technol., 21, 777–786.
Polonsky, I. N., and A. B. Davis, 2004: Lateral photon transport in dense scattering and weakly absorbing media of finite thickness: asymp-totic analysis of the space-time Green function. J. Opt. Soc. Am. A, 21, 1018-1025.
Polonsky, I. N., L. C. Labonnote and S. Cooper, 2008: Level 2 cloud optical depth product process description and interface control document (version 5). Available from
Pounder, N. L., R. J. Hogan, T. Varnai, A. Battaglia and R. F. Cahalan, 2012: A variational method to retrieve the extinction profile in liquid clouds using multiple field-of-view lidar. J. Appl. Meteorol. Climatology, 51, 350-365.
Saunders, R., P. Rayer, T. Blackmore, M. Matricardi, P. Bauer and D. Salmond, 2007: A new fast radiative transfer model – RTTOV-9. Available from: eumetsat.int/Home/ Main/Publications/Conference_and_Workshop_Proceedings/groups/cps/documents/ document/pdf_conf_p50_s2_10_saunders_p.pdf
Slingo, A., S. Nicholls and J. Schmeitz, 1982: Aircraft observations of marine stratocumulus during JASIN, Quart. J. Roy. Meteorol. Soc., 108, 833-856.
Toon, O. B., C. P. McKay, T. P. Ackerman and K. Santhanam, 1989: Rapid calculation of radiative heating rates and photodissociation rates in inhomogeneous multiple scattering atmospheres. J. Geophys. Res., 94, 16287–16301.
Wu, Z. S., and Y. P. Wang, 1991: Electromagnetic scattering for multi-layered sphere – Recursive algorithm. Radio Sci., 26, 1393-1401.
Scientific Background of the algorithm
In terms of inputs and outputs, the heritage of this algorithm comes from a number of synergy algorithms published in the last decade or two in which two instruments are combined to produce an estimate of two properties of a particular cloud type, typically water content and particle size. Since the Best Estimate algorithm is flexible enough to be run with any instrument combination, we would expect it to produce similar outputs as the best of these previous algorithms would, given the same inputs and the same microphysical assumptions. In terms of the inner workings of the Best Estimate algorithm, we draw on the narrower subset of these previous algorithms that use a variational (or “optimal estimation theory”) formulation. Two key novelties of our algorithm are (1) that it can incorporate any number of different instruments (the previous maximum in a cloud retrieval was three, e.g. radar, lidar and infrared radiometer by Delanoe and Hogan 2008), and (2) that it retrieves multiple atmospheric particle types (ice cloud, liquid cloud, rain and aerosol) simultaneously. We are not aware of an existing variational algorithm that has this capability. We now briefly review the relevant previous algorithms for each type of particle; mostly we focus on synergy algorithms, but in the case of liquid cloud we draw on a broader range of previous algorithms. A more detailed review was provided in [CASPER-PARD]. We then briefly review the variational approach in which a cost function is defined, and the quasi-Newton method of minimizing it.
1 Ice cloud
The synergy of radar and lidar in ice clouds has great potential for retrieving particle size due to the very different size dependence at the two wavelengths and the fact that spaceborne lidar often penetrates much deeper into ice clouds than liquid clouds. The main problems to overcome are (1) lidar attenuation, (2) lidar multiple scattering, and (3) consistent retrievals between regions where both radar and lidar detect the cloud and regions where only one instrument detects the cloud.
The earliest radar-lidar algorithm for ice clouds by Intrieri et al. (1993) neglected all these factors. Wang and Sassen (2002) and Okamoto et al. (2003) went further and corrected the lidar for attenuation as an initial step, before combining with radar to estimate particle size and water content. Donovan and van Lammeren (2000) showed that lidar attenuation could be corrected more reliably by using the radar information, retrieving the lidar extinction-to-backscatter ratio (constant with range) that produced a profile of effective particle size had the least variation with range. Tinel et al. (2005) took an analogous approach but seeking the extinction-to-backscatter ratio that produced a profile of effective particle number concentration that varied least with range. The properties of these algorithms were compared carefully by Hogan et al. (2006), and both have the capability for lidar multiple scattering to be incorporated.
More recently, Delanoe and Hogan (2008) developed a variational algorithm that combined these advantages with the capability to incorporate infrared radiometers, and also retrieve seamlessly between regions detected by both active instruments and regions detected by just one of these instruments. This has since been applied to CloudSat and Calipso (Delanoe and Hogan 2010). The capability to use HSRL has since been added [CASPER-ATBD], and it is this algorithm that has been the starting point scientifically for the algorithm described in this document, although it has been recoded from scratch to maximise flexibility.
2 Liquid cloud
There is a much lower level of maturity in synergistic liquid-cloud retrievals due to the problem that the lidar is rapidly attenuated and the radar return is often dominated by occasional drizzle drops (Fox and Illingworth 1997). The difficulty in using the radar or lidar backscatter directly means that often we must also make use of path constraints (such as from the radar sea-surface return beneath liquid clouds or from solar radiances), when available.
The official CloudSat product is based on Austin and Stephens (2001), who combined the radar and visible optical depth from a shortwave radiometer in a variational framework. This algorithm assumes a single-mode size distribution, i.e. neglecting larger drizzle droplets. For those drizzle-free clouds that the radar detects we would expect their properties to be retrieved quite accurately, but when drizzle is present we would expect both liquid water content (LWC) and effective radius to be overestimated. Obviously this method only works in the day when solar radiances can be used. Hawkness-Smith (2008) demonstrated that the vertically integrated LWC could be estimated from the radar path integrated attenuation (PIA) over the ocean, although since the PIA due to liquid clouds is much lower than for rain, it was necessary to pay extra care in characterizing the backscatter of the ocean surface given its dependence on wind speed. Polonsky and Davis (2004) proposed to exploit the dependence of lidar multiple scattering on optical depth to estimate optical depth. Pounder et al. (2012) have demonstrated the potential of this using a variational algorithm for retrieving optical depth and the extinction profile in liquid clouds using both multiple field-of-view lidar, and more recently with the single field-of-view Calipso lidar with its footprint on the cloud of 90 m.
Other possibilities exist that have been less explored in the literature, for example using additional constraints on the shape of the profile, taking advantage that liquid clouds often have an “adiabatic” LWC profile (e.g. Slingo et al. 1982), or using HSRL to obtain extinction coefficient near cloud top. We anticipate that an algorithm that forward models all the observations mentioned in this section, with intelligent prior constraints on the shape of the LWC profile, should be able to draw these ideas together into a robust retrieval.
3 Rain
Invariably in raining situations the lidar will have been attenuated by the cloud above, so we will be reliant on the radar return. The frequency of 94 GHz is far from ideal for rain retrievals due to non-Rayleigh scattering of the rain and attenuation by both rain and the melting layer. Nonetheless, several previous authors have proposed and implemented methods to retrieve rain rate from spaceborne cloud radar in situations where the rainfall was not too heavy.
Matrosov (2007) showed that even though reflectivity factor was not well correlated to rain rate, the decrease in reflectivity with penetration depth into the rain due to attenuation was quite well correlated and could be used to provide an estimate of rain rate. The operational CloudSat algorithm (L’Ecuyer and Stephens 2002) uses the path-integrated attenuation from the radar surface echo over the ocean to estimate the total attenuation, attributing most of it to the rain attenuation from which rain rate may be estimated. O’Connor et al. (2005) showed that in light rain beneath stratocumulus cloud, rain rate could be estimated from a combination of radar reflectivity factor and Doppler velocity, although this ignored possible non-uniform beam filling effects. Optimal estimation theory has been used in the development of algorithms intended for the Global Precipitation Mission (e.g. Grecu and Olsen 2006).
In building these ideas into a variational retrieval, there are also useful constraints that can be used but which have not been incorporated into a retrieval until now. These include the constraint that the rain rate should not vary rapidly with range (since raindrops fall fast and evaporate slowly) and that the snow flux above the melting layer should be approximately equal to the rain rate below the melting layer.
4 Variational retrieval theory
A variational approach, sometimes called “optimal estimation theory” (Rodgers 2000), defines the state vector x to contain all the variables to be retrieved, and makes use of a forward model that predicts the observations yf=H(x) that would correspond to the current contents of x. A scalar cost function, J, is defined that expresses the norm of the difference between the forward modelled observations yf and the actual observations yo, plus the norm of the difference between the state variables and any prior estimates of them. The x that minimizes the cost function is denoted the solution. The cost function is defined in section 3.4.1and methods to minimize it in section 3.4.2, in particular the quasi-Newton scheme used by the Best Estimate algorithm
1 Definition of the cost function
The cost function for this problem may be written in either summation form or matrix form as
[pic] (1)
The first term on the right-hand-side expresses the sum of the squared deviations between the observations and their forward-modelled counterparts, normalized by the error variance of the observations. The term δy is shorthand for yo – yf. In the summation form, σio is the error standard deviation of the ith observation, while in the matrix form, R represents the error covariance matrix of the observations. In practice we assume that the observational errors are uncorrelated, in which case R is diagonal with the ith element of the diagonal given by (σio)2.
The second term on the right-hand-side expresses the sum of the squared deviations between the state variables and their prior estimates, normalized by the error variance of the prior estimate. A prior estimate can be thought of as the climatological mean of that variable before the information from the observations has been taken into account. In practice, the prior error variance of many retrieved variables is large, which ensures that it does not provide a strong constraint on the solution, but can still help to stabilize the early iterations towards that solution. The term δx is shorthand for xp – x, where xp is the prior estimate of the state vector (referred to as the background in data assimilation) with its jth element denoted xjp. In the summation form, σjp is the error standard deviation of the jth prior estimate, while in the matrix form, B represents the error covariance matrix of the prior estimate. The matrix form offers more flexibility, since often it is convenient to represent error covariances between nearby prior estimates via the off-diagonal elements in B.
The third term on the right-hand-side represents a smoothness constraint, and is applied to only certain parts of the state vector (the coefficient λj is often zero). It is used to prevent noise in the observations (particularly the lidar backscatter profile) from feeding into the retrievals; see Delanoe and Hogan (2008) for a fuller discussion. It works by penalizing the curvature of x via a mechanism often referred to as Tikhonov regularization; the term inside the summation is simply the square of the finite difference form of the second derivative of x. In matrix form this is represented by the “Twomey-Tikhonov” matrix T (see Rodgers 2000 or Ansmann and Muller 2005), which has the form
[pic], (2)
where coefficient λ controls the degree of smoothness and may be optimized for a given problem using L-curve analysis (Hansen 1987, Pounder et al. 2012). For some variables it is more desirable to apply a “flatness” constraint, in which we penalize squared deviations of the gradient of the variable with height from a zero gradient. An example is rain rate, which is often close to constant with range due to the fast fall speed of raindrops compared to the time it takes them to evaporate. In this case the matrix T in (1) is given by:
[pic] . (3)
Also in this family of constraints is the one-sided gradient constraint for liquid clouds, to prevent the retrieval of profiles of liquid water content that have a gradient steeper than what would be expected for an adiabatic cloud. The details of this are provided in section 5.5.7.2.2.
2 Minimization of the cost function
At the minimum of the cost function, its gradient with respect to each element of the state vector, (xJ = (J/(x is zero. Taking the derivative of the cost function with respect to x yields this gradient (a vector):
[pic], (4)
where H is the Jacobian matrix, consisting of the rate of change of every forward-modelled observation with respect to every element of the state vector such that Hji=(xi/(yjo. In principle, this matrix may be computed at the same time as the forward model yf = H(x). Unfortunately we cannot simply set (xJ = 0 and solve for x because of the presence of the non-linear function H(x) in the definition of δy.
Newton’s method for finding the minimum is to iteratively apply the formula
[pic], (5)
where i denotes the iteration count, x1 is the first guess of the state vector, and provided that the cost function is well behaved, successive applications of this formula should yield better and better approximations to x. The term (x2J is the second derivative of the cost function, often referred to as the Hessian matrix. Note that in practice for efficiency the Hessian is not inverted but rather kept on the left-hand-side and the symmetric matrix problem is solved by a method such as Cholesky decomposition.
Since the cost function has a quadratic dependence on x, in this case the Hessian is given by
[pic]. (6)
Applying this in (5) is known as the Gauss-Newton method (e.g. Rodgers 2000). This method was used by Delanoe and Hogan (2010) in variational retrievals of ice-cloud properties from the CloudSat and Calipso instruments, and is also used in the standard CloudSat retrieval algorithms. Because it makes use of second-derivative information, the Gauss-Newton method tends to converge rapidly. In the same family of methods is the Levenberg-Marquardt method (e.g. Rodgers 2000), which can be thought of as a blend between the Gauss-Newton method and a gradient-descent method: if after a Gauss-Newton iteration the cost function does not decrease then another search direction is tried that is closer to the direction that leads to the most rapid reduction of the cost function with respect to a change to the state variables. At the final iteration of either the Gauss-Newton or Levenberg-Marquardt methods, an estimate of the error covariance matrix of the solution is given simply by the inverse of the Hessian. A disadvantage is that at each iteration we need to compute the full Jacobian matrix, which is considerably more expensive than simply running the forward model. It is usually possible (although laborious) to code up an analytic Jacobian, rather than simply perturbing each of the state variables by a small amount, and another option is to use automatic differentiation (e.g. Hogan 2013). For a forward model whose computational cost scales as the number of state variables, i.e. O(n), the computational cost of filling the corresponding Jacobian will typically be O(n2), where we have assumed that the number of state variables is of the same order as the number of observed variables. In the case of the Hogan and Battaglia (2008) wide-angle multiple scattering model (applicable to both radar and lidar), the basic forward model has a cost O(n2), and there appears to be no way to compute the corresponding Jacobian more efficiently than O(n3). For the modest values of n used here (typically no larger than 100), computing the Jacobian matrix is still viable, and it should also be pointed out that if automatic differentiation is used, computing this matrix is readily parallelized. However, as n gets larger, methods for minimizing the cost function in which the Jacobian does not need to be calculated become competitive. For very large problems, such as atmospheric data assimilation where n may be in excess of 107, neither calculation nor storage of the Jacobian matrix is feasible, and other methods need to be used.
Quasi-Newton methods use an alternative to (5) where the inverse Hessian is approximated by matrix Q:
[pic]. (7)
|[pic] |
|Figure 1: Schematic of the operation of a quasi-Newton scheme for minimizing the cost function. |
The order of operations is given in Figure 1. In the first iteration, Q is set equal to the identity matrix and thus the gradient (xJ simply provides the search direction. A line search is then performed in state space to find the minimum of the cost function along that line. This involves repeated calls to the forward model to calculate J for different values of x. Some line-search algorithms will also request the gradient of the cost function. Typically a large tolerance is put on the line search so that only a few calls to the forward model are made.
The quasi-Newton algorithm then chooses a new search direction and performs a new line search. The gradient information obtained in each outer loop for different values of x is used to build up a better and better approximation to the inverse Hessian, Q. Thus the first few iterations converge at a similar rate to a simple steepest descent scheme, after which the convergence becomes closer to a full Newton scheme (and long line searches are no longer required). We use the limited memory Broyden-Fletcher-Goldfarb-Shanno” (L-BFGS) scheme (Nocedal 1980, Liu and Nocedal 1989). This was found by Gilbert and Lemaréchal (1989) to be superior to several of its competitors for large-scale problems, and their implementation of L-BFGS is currently used in the data assimilation system of the European Centre for Medium Range Weather Forecasts.
It appears from (4) that calculating the gradient (xJ requires H to first be calculated, which is expensive. However, this can be avoided by using the adjoint method. We write the first term on the right-hand-side of (4) as
[pic]. (8)
The vector [pic] is the gradient of the observational part of the cost function with respect to the forward-modelled observations, and is very efficient to calculate since R is diagonal. The adjoint method provides a way to transform the vector (yJo into the vector (xJo without requiring the intermediate matrix H. It involves coding the adjoint of the forward model (see Giering and Kaminski 1998), which takes as input both x and (yJo, and returns (xJo. It is typically three times more expensive than the forward model itself if hand-coded adjoints are used, and up to typically four times more expensive with the Adept automatic differentiation package (Hogan 2013), so much more efficient to calculate than the full Jacobian. Thus, even though more iterations are required, this is often more than made up for by the fact that each iteration is much faster, particularly for large n.
In practice, the typical value of n in this retrieval is such that the Levenberg-Marquardt and quasi-Newton methods are of a similar order of magnitude, so both methods are available in ACM-CAP. In the case of the Levenberg-Marquardt method, this has been coded up from scratch. In the case of the quasi-Newton method, a package is used to perform the L-BFGS minimization, and in fact two are available in ACM-CAP. The first is the implementation in the GNU Scientific Library[1] available under the GNU General Public License, and the second is the liblbfgs package[2] available under the permissive MIT license. In practice we find liblbfgs to be faster due to due to its line search algorithm, so this is the default, but note that it has been modified to improve its robustness. Another option is the M1QN3[3] algorithm of Gilbert and Lemaréchel (1989), released under the GNU General Public License and used by ECMWF. All three of these libraries require the user to provide a function that calculates the cost function for a given state vector, and another that calculates its gradient for a given state vector. Section 5.5.4 describes how these functions are formulated for the current algorithm.
Justification for the selection of the algorithm
This algorithm uses all the information available to try to obtain the best possible estimate of cloud, aerosol and precipitation properties in any situation. A combined approach is essential if integral measurements (e.g. path-integrated attenuation and solar radiances) are to be used when multiple species are present in the profile. With a variational methodology, the retrievals have the prospect of being more accurate and with the most robustly derived error statistics than any alternative approach, making it attractive for use in scene construction. Moreover, this has the potential to be a flagship product for EarthCARE, exploiting its “synergy by design” ethos with the three key instruments mounted on the same platform for the first time.
A further justification is provided by considering that the ultimate goal of the EarthCARE mission is to retrieve cloud and aerosol properties along the vertical cross-section defined by the active sensors such that top of atmosphere (TOA) radiative fluxes for each (10 km)2 broad-band radiometer (BBR) cell surrounding the cross-section can be estimated to within 10 W m-2. This radiative consistency is much more likely to be achieved if the retrievals have also been forced to be consistent with as many narrowband instrument wavelengths as possible (including both active radar and lidar, and passive solar and infrared imager) as part of the formulation of the retrieval algorithm. Ensuring consistency with multiple instruments is central to the variational approach of this algorithm, so it can be considered to be a central part of achieving EarthCARE’s ultimate science goal.
This algorithm is therefore Mandatory.
Mathematical algorithm Description
1 Input variables
In addition to the configuration file described in section 5.2, the best-estimate algorithm reads in the following products, all required to be on the Joint Standard Grid (JSG) in time, but when not on the JSG in height, interpolation is performed as soon as the variables are read into memory:
• AMB-MO: ATLID-MSI-BBR Merged Observations (lidar HSRL channels averaged on to the same vertical and horizontal grid, together with the averaged MSI radiances for the same column and error estimates for all observed variables)
• ACM-TC: ATLID-CPR-MSI Target Classification (ice particles, liquid cloud droplets, rain or drizzle, melting ice, insects, aerosols, the surface of the earth, and combinations thereof)
• C-FMR: CPR Feature mask and reflectivity
• C-CD: CPR Corrected Doppler
• X-MET: meteorological and surface data
The following subsections outline the contents of each that are required by the Best Estimate algorithm (thus variables in these products that are not required will not be listed).
1 AMB-MO Merged observations
1 Product description
In addition to standard geolocation variables, this product contains:
• The lidar Mie, Rayleigh and depolarization channels, averaged to the same horizontal and vertical grid (around 100 m in height and 1 km along-track). As the lidar is at 2 degrees there might be a need for twice the lidar data, once with the pure binning and once binned relative to a reference height, e.g. 5 or 8 km, to make sure that the profiles in lidar and radar describe the exact same column of air
• “Instrument detection status” fields for the lidar, indicating whether the instrument can see some particulate target or something else. It is important to note that this contains what the instrument actually saw, not what it would be expected to see. So at the ground we would normally see a few range gates containing “ground detected” with “totally extinguished” below. If the ground signal is obscured by multiple scattering or total extinction by cloud, then this is reported as a different number, so users who wish to use the ground return only use actual observations of it.
• The radiances from each MSI channel averaged or interpolated to the same along-track spacing as the radar and lidar profiles above.
2 Variable list
The following table lists the variables contained in the product; some notes to assist in interpretation are provided in the next section.
|Variable |
|time |UTC time |s |time |real |0,1010 |
|latitude |Latitude of co-located CPR+ATLID |degree |time |real |-90, 90 |
| |footprints at the ground | | | | |
|longitude |Longitude of co-located CPR+ATLID |degree |time |real |-180, 180 or 0,|
| |footprints at the ground | | | |360 |
|height |Height above mean sea level; an |m |height |real |-1000, 30000 |
| |attribute would indicate which | | | | |
| |reference height was used in correcting| | | | |
| |for the lidar 2-degree off-nadir | | | | |
| |viewing | | | | |
|wavelength |Wavelength of centre of radiance |m |wavelength |real |0.1e-6, 15e-6 |
| |channel (solar or infrared) | | | | |
|Geographic information |
|surface_altitude |Height of surface above mean sea level |m |time |real |-300, 9000 |
|satellite_altitude |Height of satellite above mean sea |m |time |real |105, 106 |
| |level | | | | |
|day_night_flag |Day/night flag: 1 indicates the solar |none |time |bool |0, 1 |
| |zenith angle is less than 90 degrees, 0| | | | |
| |indicates that it is equal or more than| | | | |
| |90 | | | | |
|solar_zenith_angle |Angle between local zenith and the sun |degree |time |real |0, 180 |
|Measured variables |
|bscat_mie |Lidar co-polar attenuated backscatter, |m-1 sr-1 |time, height |real |1e-8, 1e-2 |
| |HSRL Mie channel | | | | |
|bscat_ray |Lidar co-polar attenuated backscatter, |m-1 sr-1 |time, height |real |1e-8, 1e-2 |
| |HSRL Rayleigh channel | | | | |
|sigma_0_bscat |Lidar surface echo |dB |time |real |TBD |
|bscat_cross |Lidar cross-polar attenuated |m-1 sr-1 |time, height |real |1e-8, 1e-2 |
| |backscatter | | | | |
|radiance |Radiances for each MSI channel |W m-2 sr-1 |time, |real |0, 20 |
| | |um-1 |wavelength | | |
| Random errors in measured variables |
|bscat_mie_error |Random measurement error in lidar |m-1 sr-1 |time, height |real |1e-8, 1e-2 |
| |backscatter for HSRL Mie channel | | | | |
|bscat_ray_error |Random measurement error in lidar |m-1 sr-1 |time, height |real |1e-8, 1e-2 |
| |backscatter, HSRL Rayleigh channel | | | | |
|sigma_0_bscat_error |Random measurement error in lidar |dB |time |real |TBD |
| |surface echo | | | | |
|bscat_cross_error |Random measurement error in lidar |none |time, height |real |0, 1 |
| |depolarisation ratio | | | | |
|radiance_error |Random measurement error in radiances |W m-2 sr-1 |time, |real |0, 20 |
| |for each MSI channel |um-1 |wavelength | | |
|Instrument detection status (8-bit unsigned integers) |
|lidar_detection_status |0 = lidar not working |none |time, height |byte |0, 6 |
| |1 = ground detected | | | | |
| |2 = totally extinguished | | | | |
| |3 = clear | | | | |
| |4 = target detected | | | | |
| |5 = molecular | | | | |
| |6 = don’t know | | | | |
3 Notes
The following points should be noted about these tables (some points are applicable to other tables in this document):
• All error variables are 1-sigma random errors.
• The variables are almost all functions of three possible dimensions: “time”, “height” and “wavelength” (the latter corresponding to the MSI channels).
• Variable types are “real” (4-byte floating point number), “int” (4-byte signed integer), “byte” (unsigned 1-byte integer) or “bool” (boolean value of 0 or 1, typically stored as a 1-byte integer).
2 ACM-TC Target Classification
This product facilitates identifies the nature of the targets in each pixel and highlights those that are bad or ambiguous, thereby informing subsequent algorithms where they should perform a retrieval (e.g. an ice cloud algorithm would only be applied to pixels containing ice cloud) and in some cases the confidence that they should assign to the observations at each pixel. In addition to facilitating application of other algorithms, it can also be used to derive cloud fraction and cloud overlap on arbitrary model-type grids.
1 Product description
In addition to standard geo-location variables, this product contains target classification fields indicating the occurrence of the various types of target that are possible: ice particles, liquid cloud droplets, rain or drizzle, melting ice, insects, aerosols, the surface of the earth, or combinations of the above. In some cases it is possible to further sub-classify this information, for example presenting the different types of aerosol present, and distinguishing “warm rain” (that originating from collision and coalescence of droplets within liquid clouds) from “cold rain” (that originating from melting ice). Note that the format that this information is stored is important: it needs to be flexible enough that combinations of types can be reported for the same pixel, while being easy to interpret by the user. Furthermore, it is important to communicate not only whether a particular type is present or not, but when it is impossible to tell because of attenuation of a particular instrument.
The classification of the atmosphere by type is described by the following variables. Since all these variables are unsigned bytes with no units that are a function of time and height, the extra table columns have been removed for clearer presentation. The rationale for this format of storage is given in the [PDD]. Note that the configuration file for the Best Estimate algorithm specifies which numbers indicate the presence of which atmospheric constituents, so if any of the following numbers is changed then it is a trivial task to update the configuration file in order to correctly interpret the new numbers.
|Classifications by type (8-bit unsigned integers) |
|aerosol_classification |0 = ground |
| |1 = none |
| |2 to 8 = types (derived from lidar only mask) |
| |9 = don’t know |
|ice_classification |0 = ground |
| |1 = none |
| |2 = ice or snow (implicit assumption that observationally they are acontinuum) |
| |3 = melting ice |
| |4 = stratospheric particles detected only by lidar (could be aerosols but impossible to tell; this could |
| |be part of aerosol_classification instead) |
| |9 = don’t know |
|rain_classification |0 = ground |
| |1 = none |
| |2 = warm rain, originating from collision and coalescence in liquid-only clouds |
| |3 = rain originating from melting ice |
| |9 = don’t know |
| |An additional class of “rain from an unknown source” may be needed, depending on the algorithm. |
|liquid_classification |0 = ground |
| |1 = none |
| |2 = warm liquid cloud (wet-bulb temperature > 0°C) |
| |3 = warm but only inferred from radar (lidar extinguished, wet-bulb temperature > 0°C) |
| |4 = supercooled |
| |9 = don’t know (lidar attenuated and T > –40°C) |
|insect_classification |0 = ground |
| |1 = none |
| |2 = yes |
| |9 = don’t know |
|convective_classification |0 = ground |
| |1 = none |
| |2 = warm convective core (wet-bulb temperature > 0°C) |
| |3 = cold convective core |
| |9 = don’t know |
|Summary classifications (8-bit unsigned integers) |
|simplified_classification |0 = ground |
| |1 = clear sky |
| |2 = liquid cloud |
| |3 = ice only |
| |4 = ice and supercooled liquid |
| |5 = ice but can’t tell if liquid too (lidar extinguished) |
| |6 = melting ice |
| |7 = rain |
| |8 = rain and liquid cloud |
| |9 = aerosol |
| |10 = insects |
| |11 = stratospheric |
| |12 = convective core |
| |13 = don’t know |
|extended_classification |0 = no targets |
| |10 = clear sky |
| |11 = liquid cloud only |
| |12 = probably liquid cloud based on radar alone |
| |13 = aerosol only (13-17 available to distinguish some basic types) |
| |18 = stratospheric feature not detected by radar |
| |19 = clear sky but lidar attenuated so might contain aerosol or liquid |
| |20 = ice only |
| |21 = ice and liquid cloud |
| |29 = ice but lidar attenuated so might also contain liquid |
| |30 = melting ice only |
| |31 = melting ice and liquid cloud |
| |39 = melting ice but lidar attenuated so might also contain liquid |
| |40 = warm rain only |
| |41 = warm rain and liquid cloud |
| |42 = warm rain, and probably liquid based on radar alone |
| |49 = warm rain but lidar attenuated so might also contain liquid |
| |50 = rain from melting ice only |
| |51 = rain from melting ice and liquid cloud |
| |59 = rain from melting ice, but lidar attenuated so liquid possible too |
| |60 = warm convective core with no liquid droplets |
| |61 = warm convective core and liquid droplets observed by lidar |
| |69 = warm convective core and lidar attenuated |
| |70 = cold convective core with no liquid droplets (e.g. T < –40(C) |
| |71 = cold convective core with liquid droplets observed by lidar |
| |79 = cold convective core and lidar attenuated |
| |80 = insects only |
| |81 = insects and liquid cloud |
| |83 = insects and aerosol (83-87 available to distinguish aerosol type) |
| |89 = insects but lidar attenuated so might also contain liquid |
| |99 = don’t know |
3 C-FMR Feature Mask and Radar Reflectivity factor
1 Product description
This product contains:
• Radar reflectivity factor
• An “Instrument detection status” fields the radar, indicating whether the instrument can see some particulate target or something else.
2 Variable list
The relevant variables are listed in the following table:
|Variable |
|Z |Radar reflectivity factor |dBZ |time, height |real |-50,40 |
|sigma_0_Z |Radar surface echo |dB |time |real |-20, 40 |
|velocity |Radar reflectivity weighted Doppler |m s-1 |time, height |real |-10, 10 |
| |velocity | | | | |
|PIA |Two-way radar path-integrated |dB |time |real |0, 60 |
| |attenuation | | | | |
| Random errors in measured variables |
|Z_error |Random measurement error in radar |dB |time, height |real |0, 2 |
| |reflectivity factor | | | | |
|sigma_0_Z_error |Random measurement error in surface |dB |time |real |0, 10 |
| |echo | | | | |
|PIA_error |Random error in two-way radar |dB |time |real |0, 20 |
| |path-integrated attenuation | | | | |
|Instrument detection status (8-bit unsigned integers) |
| |0 = radar not working |none |time, height |byte |0, 9 |
|radar_detection_status |1 = ground detected | | | | |
| |2 = totally extinguished | | | | |
| |3 = clear | | | | |
| |4 = target detected | | | | |
| |5 = dominated by multiple-scattering | | | | |
| |9 = don’t know | | | | |
4 C-CD Corrected Doppler
1 Product description
This product contains the Doppler velocity corrected for spacecraft motion, non-uniform beam filling and folding.
2 Variable list
|Variable |
|velocity |Radar reflectivity weighted Doppler |m s-1 |time, height |real |-10, 10 |
| |velocity | | | | |
| Random errors in measured variables |
|velocity_error |Random measurement error in |m s-1 |time, height |real |0, 1 |
| |reflectivity-weighted Doppler velocity | | | | |
5 X-MET Meteorological and surface data
1 Product description
This product contains:
• ECMWF profiles of temperature, pressure, humidity and ozone, for use in the radiative-transfer forward models of subsequent retrieval algorithms.
• Estimates of the surface emissivity and spectral albedo in each MSI channel beneath each profile, also for use in radiative forward models.
• Gaseous radar attenuation predicted from the ECMWF thermodynamic profiles. This would probably not be used directly by the variational retrieval, which would forward model the radar reflectivity including gaseous attenuation, but would be useful for any simpler algorithms such as an empirical expression for ice water content as a function of radar reflectivity and temperature.
• Wet-bulb temperature, Tw; this is the temperature that falling particles attain when falling through sub-saturated air. This is cooler than the actual temperature due to evaporation, and increases as the humidity drops further below saturation. This is included because ice melts at Tw = 0°C, so is useful as a first-guess for the location of the melting layer.
2 Variable list
|Variable |
|Temperature |Air temperature |K |time, height |real |180, 340 |
|Pressure |Air pressure |Pa |time, height |real |0, 110000 |
|Ozone |Ozone concentration (mass mixing ratio)|kg/kg |time, height |real |0, 0.2e-6 |
|humidity |Specific humidity |kg/kg |time, height |real |0, 0.05 |
|wet_bulb_temperature |Wet-bulb temperature |K |time, height |real |180, 340 |
|skin_temperature |Skin temperature |K |time |real |180, 340 |
|surface_pressure |Surface pressure |Pa |time |real |80000, 110000 |
|tropopause_altitude |Tropopause altitude |m |time |real |5000, 20000 |
|surface_albedo |Surface albedo in each MSI channel |none |time, |real |0, 1 |
| |(emissivity = 1 – albedo) | |wavelength | | |
|gaseous_radar |Two-way radar attenuation from space |dB |time, height |real |0, 15 |
|_attenuation |due to atmospheric gases | | | | |
|Random errors in model variables |
|temperature_error |Air temperature random error |K |time, height |real |0, 20 |
|pressure_error |Air pressure random error |Pa |time, height |real |0, 5000 |
|ozone_error |Ozone concentration random error |kg/kg |time, height |real |0, 0.1e-6 |
|humidity_error |Specific humidity random error |kg/kg |time, height |real |0, 0.02 |
|wet_bulb_temperature |Error in wet-bulb temperature |K |time, height |real |180, 340 |
|_error | | | | | |
|skin_temperature_error |Skin temperature random error |K |time |real |0, 20 |
|surface_pressure_error |Surface pressure random error |Pa |time |real |0, 5000 |
|tropopause_altitude |Tropopause altitude random error |m |time |real |0, 2000 |
|_error | | | | | |
|surface_albedo_error |Surface albedo random error |none |time, |real |0, 0.2 |
| | | |wavelength | | |
|gaseous_radar |Random error in two-way radar |dB |time, height |real |0, 5 |
|_attenuation_error |attenuation from space due to | | | | |
| |atmospheric gases | | | | |
2 Configuration Parameters
The algorithm is designed to be fully flexible with regard to what observations are used, what atmospheric constituents are retrieved, and the assumptions and constraints used in the calculations. The algorithm is even quite flexible with regard to the format and variable names in the input data file. In this section we describe how the behaviour is configured, which can be done without recompiling.
The code may be configured via a single configuration file that essentially consists of a set of key-value pairs grouped into named sections. This could be implemented in XML, or as illustrated in the SUM, a simpler ASCII format. The global parameters to be set are as shown in Table 1.
|Table 1. Valid entries in the global part of the algorithm configuration file. All are character strings (represented by the courier font) except|
|where indicated. |
|Key name |Value |
|minimizer |Normally LBFGS (the limited memory Broyden-Fletcher-Goldfarb-Shanno quasi-Newton scheme using the |
| |liblbfgs library), but other valid options are LBFGS_GSL (the GNU Scientific Library implementation|
| |of this algorithm), simplex (the GNU Scientific Library implementation of the gradient-free |
| |“Simplex” method) and levenberg_marquardt for the Levenberg-Marquardt method. |
|max_iterations |(Integer) The maximum number of iterations to perform per retrieved profile. |
|converged_gradient_norm |(Real) When the L2 norm of the gradient of the cost function with respect to the state vector falls|
| |below this number, convergence is deemed to have occurred. |
|max_step_size |(Real) Maximum step-size in state space between one iteration and the next (Levenberg-Marquardt |
| |only) |
|wolfe_coefficients |(Two real numbers) The values of coefficients c1 and c2 in the Wolfe conditions governing |
| |convergence of the line search used by the liblbfgs library (see the Wikipedia entry “Wolfe |
| |conditions” for the exact meaning and suitable values of these two coefficients). |
|observations |A list of strings naming the observations to be assimilated, each referring to a named section of |
| |the file where the details of that observation type are given. For EarthCARE the strings are likely|
| |to be CPR ATLID MSI. |
|constituents |A list of strings naming the constituents to be retrieved, each referring to a named section of the|
| |file where the details of that constituent type are given. For EarthCARE the strings are |
| |LiquidCloud IceCloud Rain Aerosol MeltingIce. |
|time_name, altitude_name, height_name, |The names of the variables in the input file that contain the time, altitude of the satellite, |
|pressure_name, temperature_name, |height of each range-gate, pressure at each range gate, temperature, specific humidity, ozone mass |
|specific_humidity_name, ozone_mmr_name, |mixing ratio and surface temperature. An additional set of optional variables altitude_scaling, |
|surface_temperature_name |height_scaling, etc., provide the scalings that need to apply to convert the variables into the SI |
| |units expected by the algorithm. |
|input, output |The name of the input and output data files (which for maximum flexibility could be either NetCDF |
| |or HDF format). |
|environment_input |The name of the input file containing the meteorological and surface variables, since this may be |
| |in a different file from the one containing observations. |
The global keys observations and constituents refer to later sections of the configuration where specific instruments and retrieved constituents are defined. The sections describing each constituent contain keys listed in Table 2.
Within each constituent section are subsections defining each type of state variable needed to describe that constituent, as declared in the state_variables field. The valid entries for these subsections are given in Table 3.
|Table 2. Valid entries in a section of the configuration describing a retrieved constituent. See Table 3 for a description of how the |
|individual variables describing a particular retrieved constituent are represented in height. |
|Key name |Value |
|type |One of ice, liquid, rain, aerosol or simple_melting_layer. |
|minimum_size, maximum_size |(Real numbers) A range of particle size distributions are created, composed of different |
| |concentrations of particles in the single-particle scattering files (described later). These keys |
| |define the minimum and maximum median sizes of these distributions. |
|num_sizes |(Integer) The number of median sizes (i.e. the number of different size distributions) to create. |
|distribution_shape |The shape of the size distribution: can be exponential, gamma, lognormal, field or delanoe. |
|shape_factor or gamma_order |(Real number) A number describing some property (typically the relative width) of the distribution;|
| |for gamma distributions this is the “mu” parameter (zero is exponential and larger values yield |
| |narrower distributions), while for lognormal distributions it is the standard deviation of the |
| |natural logarithm of particle sizes. |
|classification_name |The name of the variable in the input file that contains a target classification. |
|classification_codes |(Integer list) A list of codes that if contained in the classification variable would indicate the |
| |presence of this consituent. |
|state_variables |A list of names of subsections within the section for this particular constituent that describe the|
| |details of each of the types of state variable to use to describe the constituent |
|derived_variables |A list of names of time-height variables to be derived from the state variables and reported in the|
| |output file, along with their errors. Valid values are water_content (or mass_content), extinction |
| |(in the geometric optics approximation), mass_flux, effective_radius, D0 |
|integrated_variables |A list of names of vertical integrals to be derived from the state variables and reported in the |
| |output file, along with their errors. Valid values are optical_depth (in the geometric optics |
| |approximation) and water_path. |
|Table 3. Each constituent is described by up to four different types of state variable as listed in Table 8. The present table provides the valid |
|entries in a subsection of a constituent section of the configuration describing one of these types of state variables. |
|Key name |Value |
|type |The type of state variable, which may be one of mass_content (or water_content), |
| |extinction, mass_flux (or precipitation_rate), median_volume_diameter, |
| |number_concentration, normalized_number_concentration, primed_number_concentration, |
| |fraction_scaling (for riming) or lidar_backscatter_extinction_ratio |
|representation |One of predefined (the variable is not to be retrieved but fixed at the prior value), |
| |constant (a single value will be retrieved constant with height), direct (separate values |
| |will be retrieved at each range gate), cubic_spline (the profile will be represented by a |
| |cubic spline, the resolution of which is determined by basis_function_spacing)or linear |
| |(the variable will be assumed to vary linearly with height, and two variables will be |
| |retrieved describing its mean and its gradient with height). |
|prior |(Real number or numbers) The a-priori value for the variable. If a vector of values are |
| |provided then the prior is assumed to be temperature dependent, and the temperature_interp |
| |variable is needed to provide the corresponding temperature coordinate. |
|prior_error |(Real number) The 1-sigma error in the a-priori values. |
|temperature_interp |(Real numbers) The temperature in Kelvin at which the prior values are defined for |
| |interpolation. If this entry is omitted then the prior is assumed to be constant with |
| |temperature. |
|basis_function_spacing |(Real number) Only used if representation=cubic_spline; this entry controls how many cubic |
| |spline control points are used by setting the target number of range gates between each |
| |control point. A larger number leads to a smaller number of state variables and so more |
| |rapid convergence, but a smoother retrieved profile. |
|prior_decorrelation_length |(Real number) Only used if representation=direct or cubic_spline; this entry specifies the |
| |decorrelation length in metres to be used to define the prior error covariance matrix; this|
| |is denoted z0 in section 8.3. |
|max |(Real number) A term is added to the cost function that penalizes values above this number |
| |(if provided) with a strength specified by max_factor. |
|max_factor |(Real number) Factor scaling the strength of penalty for a retrieved variable x exceeding |
| |max. Thus for x>max, the term added to the cost function is max_factor × (x – max)2. |
|smoothness_factor |(Real number) Only used if representation=direct or cubic_spline; this entry specifies the |
| |strength of the smoothness constraint in the cost function, if required. |
|flatness_factor |(Real number) Only used if representation=direct or cubic_spline; this entry specifies the |
| |strength of the flatness constraint, if required. |
|gradient_factor |(Real number) Only used for liquid water content, this entry quantifies the strength of the|
| |gradient constraint in the cost function. |
|first_guess |(Real number) If provided, this will be used as the first guess of this state variable, |
| |rather than the prior. |
|gradient_factor |(Real number, liquid water content only) A one-sided gradient constraint should be used |
| |with number scaling the magnitude of the constraint in the cost function. |
|temporal_smoothing_error |(Real number) A Kalman smoother is to be employed for this variable to smooth between rays,|
| |by using the previous ray as a constraint on the current one, taking the retrieval error of|
| |the previous ray but inflating it in quadrature with this number (to represent spatial |
| |decorrelation between rays). |
|adiabatic_prior_peak |(Real number, liquid water content only) The prior liquid water content (LWC) will be |
| |adiabatic with this as the LWC value at the top of each independent layer of liquid cloud. |
|adiabatic_prior_slope_scaling |(Real number, liquid water content only) The prior liquid water content will be scaled down|
| |from adiabatic by this factor (usually less than 1). |
The sections of the configuration file describing each observation contain keys listed in Table 4. Within the observation sections is a variables key, which contains a list of the variables from that instrument to be assimilated (e.g. reflectivity factor and Doppler velocity). These in turn are described in a further set of subsections, as described in Table 5. The variable subsections contain, for example, the name of the variable in the input file, the name of the corresponding variable that contains its error, and in the case of active instruments it also contains a detection status key that contains information on which range gates in a profile contain a valid signal.
|Table 4. Valid entries in a section of the configuration describing an observation type. |
|Key name |Value |
|type |One of radar, lidar, ir_radiometer or sw_radiometer |
|wavelength |(Real number or numbers) The wavelength in metres; for imagers/radiometers multiple wavelengths might be|
| |simulated simultaneously in which case this would be a list of numbers. |
|beam_divergence |(Real number) The 1/e half-width of the transmitter beam in radians |
|receiver_fov |(Real number) The receiver field-of-view half-width in radians. If omitted this will be assumed to be |
| |the same as the beam divergence (which is the case for monostatic radar). |
|reference_dielectric_factor |(Real number) The reference dielectric factor |K0|2 used to define the calibration of a radar (e.g. 0.75|
| |for CloudSat). |
|wide_angle_multiple_scattering |(0 or 1) Option to switch on or off the forward modelling of multiple scattering. |
|variables |A list of the names of the individual observed variables that will be used, referring to a named |
| |subsection of the file describing that variable (see Table 4). |
|Table 5. Valid entries in a subsection of an observation section of the configuration describing an input variable. |
|Key name |Value |
|type |One of reflectivity_factor (assumed units dBZ), |
| |apparent_Doppler_velocity (assumed units m s-1), |
| |apparent_backscatter, apparent_backscatter_mie, |
| |apparent_backscatter_rayleigh, |
| |apparent_backscatter_cross (all with assumed units m-1 sr-1), |
| |apparent_surface_backscatter (assumed units dB), |
| |radiance (assumed units W m-2 sr-1 m-1), |
| |brightness_temperature (assumed units K) |
|name |The name of the variable in the input data file. |
|scaling |A real number that the variable should be scaled by to obtain the units expected according to its type |
| |(optional). |
|error_name |The name of the variable containing the corresponding 1-sigma error in the data file. |
|error_scaling |A real number that the error variable should be scaled by (optional). |
|detection_name |The name of the integer variable describing whether a valid signal is present (optional: if missing then|
| |a valid signal is always assumed present) |
|detection_values |A list of integers that if contained in the detection variable indicate a valid signal. |
|molecular_values |A list of integers that if contained in the detection variable indicate only molecular returns are |
| |detected (lidar only) |
|ground_values |A list of integers that if contained in the detection variable indicate the ground. |
|forward_model _error |(Real number) The error to add in quadrature to the observational error in the input file, representing |
| |the forward modelling error. |
|forward_model _log_error |(Real number) The fractional error to add in quadrature to the observational error in the input file, |
| |representing the forward modelling error. |
3 Output variables
This section describes the full product produced by the Cloud, Aerosol and Precipitation Best Estimate algorithm. It should be noted that the algorithm can also serve as a single-constituent algorithm (e.g. an ice-cloud retrieval), or as a single-instrument algorithm; in this case only the variables appropriate to the constituent retrieved or the instrument used would be reported in the output.
1 Product description
This product contains the following at the same resolution as the input merged observations, i.e. on the Joint Standard Grid (100 m vertically and around 1 km along-track):
• 1D coordinates variables: These are the variables that contain information on geolocation, time etc. The other variables in the product are in general functions of these coordinates. The dimensions of the coordinates are explicit in the column “dimension” (see below).
• Target classification information largely copied from the product described in section 5.1.2, and including details such as the location of supercooled liquid clouds
• Ice cloud properties: Vertical profiles of ice water content, ice effective radius and ice extinction coefficient; measures of the mean terminal velocity, equivalent snowfall rate and a measure of ice particle density may also be estimated in some situations
• Liquid cloud properties: Vertical profiles of liquid water content, liquid effective radius and liquid extinction coefficient
• Rain properties: Vertical profiles of rain rate and mean raindrop size, in both rain originating from melting ice and the light rain or drizzle falling from stratocumulus clouds
• Melting layer properties: Water path and radar two-way attenuation
• Aerosol properties: Some measure of aerosol type (copied from the classification product), aerosol extinction coefficient, total number concentration and median volumetric diameter
• Forward modelled observations at final iteration. These variables are forward modelled from the state vector at the final iteration of the algorithm, and hence should closely match the observations, except that the observational noise should not be matched. Agreement with the observations provides a useful confirmation that the retrieval has converged.
• Measures of convergence: We include variables that indicate how well the solution has converged, such as normalized “chi-squared”, which is the sum of [the squared differences between the observations and the forward model at the final observation divided by the observational error variance] divided by the total number of measurements.
• Status flags: These flags inform the user about the framework in which the retrieval was performed, e.g. the existence or not of targets (retrieval_flag), or the type of configuration used by the program (instrument_flag).
• Retrieval errors: These are the 1-sigma errors in the variables directly retrieved (be aware that often the errors are in the natural logarithm of the quantity being retrieved).
• Errors in derived variables: These are the corresponding 1-sigma errors in the derived quantities such as ice water content, liquid water content etc.
A future version of this product could contain variables derived from the Doppler radar, such as measurement of convective motions. This would probably be taken from the radar-only product “Vertical Motion”.
2 Variable List
The reader is referred to the [PDD] for the full list of output variables; the list is very long so it seems unnecessary to repeat it here and maintain two versions.
4 Algorithm overview and flow chart
A top-level description of the proposed algorithm is provided in Figure 2. The algorithm first loads the global configuration file then the particle scattering properties required for each of the instrument wavelengths and each of the particle types to be retrieved. Details of this step are given in section 5.5.2
The algorithm then analyzes each ray (i.e. each vertical profile) of the input data in turn (Box 2 in Figure 2). From the Target Classification product, it determines which atmospheric constituents (ice, liquid etc) are to be retrieved at each gate. Each constituent is described by a set of state variables, for example the liquid water content at each range gate (i.e. each grid-point in the vertical) containing liquid cloud. The variables used for each different constituent are described in section 5.5.8, and are concatenated into a state vector, x, with a total of n elements. Initially each state variable is assigned its first-guess value. From the Merged Observations product the algorithm extracts the observed values that will be used in the retrieval and concatenates them into an observation vector, yo, with a total of m elements. In the equations that follow, all vectors are column vectors.
|[pic] |
|Figure 2: Top-level flowchart of the operations carried out by the algorithm. |
The variational approach described in section 3.4 is used to find the optimum x for that ray: the L-BFGS method (indicated by Box 3) is used to minimize a cost function defined by (1), making use of a forward model that predicts the observations yf=H(x) that would correspond to the current contents of x. This is where most of the mathematical work lies and is described in section 5.5.5.
When the solution has been found, the associated error covariance matrix is calculated (Box 4 of Figure 1). This step is described in section 5.5.6. The retrieved properties of each atmospheric constituent are stored in the output file, along with their error. The algorithm then proceeds to the next ray of data.
5 Algorithm definition
The algorithm is designed to be extremely flexible, allowing arbitrary combinations of atmospheric constituents to be retrieved from arbitrary combinations of instruments. The combination required in a particular retrieval is specified by the user in a configuration file so that there is no need to recompile for each different instrument or constituent combination. Even though the algorithm is designed for EarthCARE, this flexibility will allow it to be tested on a very wide range of ground-based, airborne and spaceborne instrument combinations. An additional requirement is for the algorithm to be easily extensible, for example, to add new types of observations or retrieved atmospheric constituents. The only way to achieve this elegantly is to use an object-oriented language that supports inheritance and polymorphism, such as C++ or Fortran 2003. Section 5.5.1 describes the class hierarchy, and section 5.5.2 describes how the objects are created.
Readers interested purely in the mathematical steps of the core part of the algorithm should skip to sections 5.5.3-5.5.6.
1 Class hierarchy
We now describe the main objects that this algorithm needs. Some readers may feel that this is an implementation detail and detracts from description of the mathematical operations. However, we believe that the complexity of the data structures that need to be constructed, and the many ways in which they interact, means that a full description in terms of objects is necessary, otherwise this document would not serve its purpose as a guide to building the code to implement the algorithm.
Figure 3 depicts the main class hierarchy of the algorithm. The State class contains the interface to the top-level program. For example, after the object tree has been constructed and the observational data from a single ray have been read-in, a call to its member function “solve” will perform the retrieval using as many calls to the forward model as are required.
|[pic] |
|Figure 3: Main class hierarchy. The arrows with solid heads indicate that a class contains pointers to one or more of the |
|subclasses that the arrow points to. A class that points to another class with an open-headed arrow inherits from it. White |
|boxes denote abstract classes that define a generic interface but cannot be instantiated, while the grey boxes denote concrete|
|classes that can be instantiated. |
The State object contains a list of observations to assimilate and constituents to retrieve. These observations and constituents could be of any type, the details of which do not need to be known by State. All State requires is that they present a common interface, which we define in abstract base classes called Observation and Constituent. The minimizer makes repeated calls to a function that could be called State::calc_cost_function, which returns the value of the cost function for a particular trial state vector x. Looking at the definition of the cost function in (1), the first term on the right-hand-side is associated with the observations while the second and third terms are associated with the constituents (via the prior estimate xp, which is specific to each constituent). It can be seen from the summation form of (1) that each of the three terms on the right-hand-side can be thought of as the sum of the individual contribution to the cost function from each observation and each constituent. Therefore, State::calc_cost_function is implemented by calling Observation::calc_cost_function for each observation and Constituent::calc_cost_function for each constituent, and summing the results. This is where polymorphism comes in. Unknown to State, each Observation and Constituent object is actually an instance of a more specific object such as Radar or Aerosol. So when State calls Observation::calc_cost_function, it is the specific function for that observation type that is called, such as Radar::calc_cost_function.
2 Global initialization
This section describes the initialization of the algorithm carried out before any observational data are read in, depicted by Box 1 in Figure 1. The actions performed are:
1. Read the global configuration file and use this information to build the object lists (observations and constituents) in memory and configure each object: see section 5.5.2.1.
2. Create the necessary scattering look-up tables. This involves reading the appropriate particle scattering file for each observation-constituent pair, and computing look-up tables for full size distributions of particles. Section 5.5.2.2 describes the contents and behaviour of these look-up tables while section 5.5.2.3 describes how they are created.
1 Global configuration and building object lists
Section 5.2 described the contents of the global configuration file, and an example implementation is given in the SUM. These are used to create a list of Observation and Constituent objects corresponding to the observations and constituents entries in the configuration file. The configuration information for each observation and constituent is described in Table 2 and Table 4. These sections of the configuration file contain a type entry that is used by the program to decide which specific type of object to create. Thus constituent type entries ice, liquid, rain and aerosol will result in IceCloud, LiquidCloud, Rain and Aerosol objects being created. Likewise, observation type entries radar, lidar, ir_radiometer and sw_radiometer result in Radar, Lidar, InfraredRadiometer and ShortwaveRadiometer objects being created. This flexibility allows data from more than one instrument of the same type to be assimilated simultaneously.
Each Observation can be thought of as an instrument that may measure several variables (e.g. radar reflectivity factor and Doppler velocity), but to forward model these variables requires just one call to the relevant forward model. The variables required are listed in the variables entry for the observation, which are in turn configured by the relevant subsection of the configuration file outlined in Table 5.
2 Scattering look-up tables
After the configuration information has been read in and the corresponding objects created, the scattering look-up tables are created. These are required within the forward model to convert a meteorological description of the cloud (e.g. liquid water content and particle number concentration) into the variables that describe the interaction of radiation at a particular wavelength with the cloud particles. Since particle scattering calculations are slow, they are performed offline; section 8.2 describes how a library of particle scattering properties is created and stored. A further optimization is to create look-up tables once when the algorithm is first started and then use them for all the profiles in the file, rather than recomputing them for every profile. For each observation-constituent combination, it is necessary create a separate look-up table in order that the forward model for that observation can calculate the contribution from that constituent to the profile of scattering properties before performing the radiative transfer to calculate what the instrument would see.
Before describing how look-up tables are created, it is necessary to describe their interface and how they work internally. The interface to the look-up table is common for all constituent types: vectors of up to four variables are provided as input, corresponding to the properties of that constituent or the atmosphere at each vertical level where that constituent is present. The properties are:
• The abscissa of the look-up table, X. This is usually an extensive variable that is identical to or close to one of the state variables for that constituent. Permitted types are MASS_CONTENT, (e.g. for liquid clouds where the main state variable is the mass per unit volume) GEOMETRIC_EXTINCTION (equal to the extinction coefficient in the geometric optics regime, the main state variable for ice clouds), MASS_FLUX (e.g. for rain where the main state variable is rain rate) and MEAN_DIAMETER (one of the state variables for aerosol). It is important that for fixed number concentration, this variable increases monotonically with particle size. It is required that this variable does not vary with any of the others.
• The ordinate of the look-up table, Y. This provides an additional descriptor of the particle type and for most constituent types it is not used. An example of its use would be in ice clouds where the density of the particles is allowed to vary, in which case this can be RIMING_INDEX, a variable that scales the exponent in the mass diameter relationship between the default, corresponding to a density-size relationship appropriate for unrimed aggregates from the literature, e.g. Brown and Francis (1995), and a number corresponding to solid ice, such that the values in between correspond to different degrees of riming. Because of the requirement that each input variable is independent, if Y is used as a density multiplication factor then X cannot be related to the mass, so must be GEOMETRIC_EXTINCTION or MEAN_DIAMETER.
• The temperature, T. Liquid scattering properties at radar wavelengths vary with temperature, an effect that needs to be included. When there is no significant variation of the scattering properties with temperature, this input variable will be ignored.
• A measure of number concentration, N. When this is varied, the size of the particles does not change, simply their number concentration. Permitted types are (total) NUMBER_CONCENTRATION, suitable for liquid clouds where the droplet number concentration is a retrieved state variable assumed constant for each contiguous liquid layer, and NORMALIZED_NUMBER_CONCENTRATION, used by ice clouds and rain that are more usefully described by normalized size distributions (Illingworth and Blackman 2002, Delanoe et al. 2005).
The look-up table then outputs vectors of each of the scattering properties that are required to forward model the observation. The list of possible scattering properties that can be forward modelled is given in Table 6, followed by the scattering properties required by each observation type.
|Table 6. Scattering properties that are contained in a look-up table. Details of how these are computed are given in section 5.5.2.3 and also in |
|section 8.2. Note that since not all of these scattering properties are required for every observation-constituent pair, only those required will |
|be computed: see the following table. |
|Scattering property name |Symbol |Description |
|EXTINCTION_COEFFT |α |Extinction coefficient in m-1. |
|BACKSCATTER_COEFFT |β |Backscatter coefficient in m-1 sr-1. Note that this is the unattenuated |
| | |value, unlike the attenuated_backscatter observation variable. Moreover, |
| | |for radar forward modelling, the code works internally with this variable |
| | |and converts to radar reflectivity factor at the final step when comparing |
| | |to the observed values. Usually BACKSCATTER_EXTINCTION_RATIO is requested |
| | |instead. |
|SCATTERING_COEFFT |βs |Scattering coefficient in m-1, i.e. the extinction coefficient multiplied |
| | |by the single scatter albedo. Usually SINGLE_SCATTERING_ALBEDO is requested|
| | |instead. |
|BACKSCATTER_EXTINCTION_RATIO |s |The ratio of the backscatter and extinction coefficients in sr-1. |
|SINGLE_SCATTERING_ALBEDO |ω |The ratio of the scattering and extinction coefficients. This is needed for|
| | |radiance calculations, but is only needed for active instruments when |
| | |wide-angle multiple scattering is to be forward modelled. |
|ASYMMETRY_FACTOR |g |The mean of the cosine of the scattering angle. This is needed for radiance|
| | |calculations, but is only needed for active instruments when wide-angle |
| | |multiple scattering is to be forward modelled. |
|DOPPLER_VELOCITY |vd |The backscatter-weighted terminal velocity in m s-1. |
|EQUIVALENT_AREA_RADIUS |ra |The area-weighted mean radius (in m) of the particle distribution, required|
| | |to forward model small-angle multiple scattering for lidar by the Hogan |
| | |(2008) model, in order to predict the width of the forward scattered lobe. |
|MIE_DROPLET_FRACTION |f |The fraction of the backscatter coefficient that is attributable to |
| | |Mie-scattering liquid droplets. This plays a role in multiple scattering |
| | |calculations at lidar wavelengths, because the phase function for droplets |
| | |is peaked in the backscatter direction, while for ice particles it is flat.|
| | |Therefore, f=1 for droplets and f=0 for ice, and it is only when the total |
| | |scattering profile is computed, in which ice and liquid may coexist, that |
| | |an intermediate value of f can be produced. |
|Table 7. Scattering properties that are needed by each type of observation in the case of EarthCARE. |
|Observation |Required scattering properties |
|Radar |Particulates only: α, s, vd; if radar multiple scattering is to be simulated then also ω and g; the extinction|
| |coefficient of air is computed but kept separate: αair |
|Lidar |Particulates only: α, s, ω, g, ra and f; the extinction coefficient of air is computed but kept separate: αair|
|Infrared imager channels |Particulates only, separate values for each wavelength: α, ω, g |
|Shortwave imager channels |Particulates only, separate values for each wavelength: α, ω, g, delta-M-scaled Legendre expansion of the |
| |phase function |
Naturally, four-dimensional table look-up and interpolation can be slow so two optimizations are employed. Firstly, nearest-neighbour interpolation is used in temperature space; the look-up table holds the data at a number of different temperatures and the nearest one is used. Secondly, we know that variations in N simply scale extensive variables E (such as extinction coefficient) up and down while leaving intensive variables I (such as asymmetry factor) unchanged. Therefore, all extensive variables in the look-up table are stored divided by N to convert them into intensive variables. For example, for liquid clouds, the abscissa A is liquid water content and the number concentration variable N is total droplet number concentration. Internally the extinction coefficient at a particular wavelength, an extensive variable E, would be stored as a one-dimensional look-up table of E/N versus A/N, while the asymmetry factor, an intensive variable I would be stored as I versus A/N. So to calculate E for given A and N, the look-up table performs the look up for A/N and then multiplies the resulting E/N by N. If the abscissa were MEAN_DIAMETER, an intensive variable, then the look-up tables would be versus A rather than A/N. Thus if no ordinate variable is present, a three-dimensional look-up (A, T and N) is achieved with only one linear interpolation. If an ordinate variable is present then bi-linear interpolation is performed versus A/N and Y (or A and Y if A is intensive).
3 Creating scattering look-up tables
We now describe how the look-up tables are created. If we have m observations and n constituents then we require m×n different look-up tables. There are two exceptions. Firstly, if a particular constituent is essentially invisible to a particular observation (e.g. aerosols to radar) then it can be omitted. Secondly, the BackgroundConstituent type of constituent describes the properties of atmospheric gases that may scatter and absorb at a particular wavelength, but does not use look-up tables to calculate these properties. Since the look-up tables are required by the observations in their forward models, they should be stored with the data for the observation to which they refer. The information necessary to initialize such an object is as follows; note that all of the following except the wavelength are specific to the constituent:
1. The wavelength of the observation.
2. A reference density ρr used in the calculation of effective radius (1000.0 kg m-3 for LiquidCloud and Rain constituents, and 917.0 kg m-3 for IceCloud).
3. A vector of the mean particle sizes of the distributions used in the look-up table.
4. The type of size distribution normalization: either NUMBER_CONCENTRATION, or NORMALIZED_NUMBER_CONCENTRATION.
5. The type of the size distribution, from amongst EXPONENTIAL, GAMMA, LOGNORMAL, FIELD or DELANOE.
6. The shape parameter describing the form of some of the size distributions.
7. The type of abscissa: MASS_CONTENT, GEOMETRIC_EXTINCTION, MASS_FLUX or MEAN_DIAMETER.
8. A boolean flag indicating whether the abscissa will be stored in natural-logarithmic form (and will accept input abscissas in logarithmic form). This is useful in the common situation that the relevant state variable is actually the natural logarithm of a quantity such as extinction coefficient or water content and can be passed directly to the look-up table.
9. A boolean flag indicating whether the number concentration variable will be stored in natural-logarithmic form (as it usually is in the state vector).
10. The type of the ordinate vector: currently it can be NONE or FRACTION_SCALING (where the mass/density of the particle is scaled by this ordinate). In future, additional ordinate types may be added, for example to represent various aerosol properties.
11. A vector of ordinate values, if used.
Most of these variables are prescribed by the specific constituents, but some can be specified in the relevant constituent section of the configuration file, as described in Table 2. Each constituent accepts the keys minimum_size, maximum_size (both in m) and num_sizes, to indicate the range of mean sizes to be used in the abscissa of the look-up table. Note that the definition of “size” (e.g. diameter or radius) depends on the particle type. These three variables are used to create a logarithmically spaced vector of mean sizes that is used to create the vector given in item 3 above. Additionally the distribution_shape key can take the value exponential, gamma, lognormal, field or delanoe, while the numerical factor describing the width of that particular type of size distribution is given by shape_factor.
Before the algorithm is run (i.e. outside of the unified retrieval program), particle scattering calculations have to be performed for all particle types and wavelengths to be simulated, with the data for each combination stored as a NetCDF file. The data are stored as a function of up to 4 coordinate variables: wavelength, temperature, diameter and mass fraction. The last of these is specific to ice and represents the treatment of ice particles as homogeneous mixtures of ice and air, in which case the mass fraction variable is the ratio of the density of the particle to the density of solid ice. Often these coordinate variables are singletons, for example, for liquid water droplet scattering by a visible lidar, the file would contain scattering properties at just one wavelength, one temperature (since the refractive index of water can be treated as temperature-independent in the visible) and one mass fraction (we don’t consider liquid-air mixtures). The particle properties stored in the file are then cross-sectional area, mass, terminal fall speed, refractive index, extinction cross-section, scattering cross-section, backscatter cross-section (in m2) and asymmetry factor, all in SI units. A program to do these calculations using Mie theory is described in section 8.2. The terminal fall-speed of ice particles is calculated following Heymsfield and Westbrook (2010), while for liquid drops we use Beard (1976). Calculations are performed for a reference air density of 1 kg m-3, and scaled depending on the actual air density in the forward model. A configuration file, which is in ASCII and called scattering_database.cfg in the current implementation, contains a list of the available particle scattering files, the medium of the particles in that file (which may be ice, liquid-water or aerosol), and the wavelength range over which the scattering calculations in the file are valid. Since freely available code exists for both Mie scattering calculations and terminal fall-speed calculations, we treat these components as external libraries and no description of their inner workings is needed here.
The first task in creating the look-up table is to locate and load the relevant NetCDF file, which is done by searching scattering_database.cfg for files with a matching wavelength and medium. These particle scattering properties are then integrated across size distributions with the range of mean sizes D0 requested. If the NetCDF file contains n particle sizes where the ith size has diameter Di at which point the spacing in diameter space is ΔDi, then we may write the size distributions in discrete form in terms of the number concentration Ni (in m-3) of particles in size bin i, as follows:
• Exponential distribution: [pic], where here D0 is interpreted as the mean diameter;
• Gamma distribution: [pic], where here D0 is interpreted as the median volumetric diameter, and μ is the shape parameter;
The two measures of number concentration are given by:
• Total number concentration: [pic]
• Normalized number concentration: [pic]
Extensive variables such as extinction coefficient, α, are defined simply as [pic], where σi is the extinction cross-section of particles in size bin i. Other extensive variables in Table 6 are the backscatter coefficient, β, which is defined analogously as the sum over particle backscatter cross-section, and the scattering coefficient which is defined as the sum over particle scattering cross-section. The variable N0 in the definition of the number distribution provides a scaling for the distribution, but it is never required since we store the extensive variables divided by NT or N0*, which are themselves proportional to N0 and so the N0 cancels. For example, for ice clouds we would store a look-up table of α/N0* versus particle size, and when a forward model wants α as a function of size and N0*, then the 1D look-up table would provide α/N0* for the given size, which would then be multiplied by N0* to get α.
Intensive variables are also not dependent on N0. For example, we define the following measures of radius:
• Effective radius: [pic];
• Equivalent-area radius: [pic],
where ρr is the reference density for the medium, m is the particle mass and A is the particle cross-sectional area. Regarding the other intensive variables: single scattering albedo is defined as scattering coefficient divided by extinction coefficient, backscatter-to-extinction ratio is the backscatter coefficient divided by the extinction coefficient, while asymmetry factor for an entire distribution is defined as the scattering-coefficient-weighted average over the asymmetry factors of the individual particles.
|Table 8. List of state variables used to describe each of the five atmospheric species. As most state variables vary over several orders of |
|magnitude, the natural logarithm of each one is used in the state vector x, to improve the convergence behaviour and to prevent negative values. |
|The strings in brackets in the first column indicate the name used to configure this variable in the configuration file (see section 5.2). The |
|justification of this choice of state variables is given in section 5.5.7. |
|State variable |Representation with height |A-priori and error of natural |Unit |
| | |logarithm | |
|Ice clouds and snow |
|Visible extinction coefficient, αv |One variable per pixel but with smoothness |ln(10-6) ±5.0 |m-1 |
|(extinction) |constraint κ = 100 | | |
|N0’ = N0*/α0.6 (primed_number |Cubic spline basis functions with vertical |ln(N0’) = A+BT[°C]; A=22.46, |m-3.4 |
|_concentration) |decorrelation length 1 km |B=–.08932; error 1.0 | |
|Lidar backscatter-to-extinction ratio, S|One variable per pixel but with smoothness |ln(20)±0.5 |sr |
|(lidar_backscatter _extinction_ratio) |constraint κ = 200 | | |
|Riming factor (fraction_scaling) |Likely a single value per profile |Unity |None |
|Liquid clouds | | | |
|Liquid water content, LWC |One variable per pixel but with gradient constraint | |kg m-3 |
|(water_content) | | | |
|Total number concentration, NT |One value per liquid layer |Marine: ln(74×106); Continental: |m-3 |
|(number_concentration) | |ln(288×106); error 0.6 | |
|Rain | | | |
|Rain rate, R (precipitation_rate) |Cubic spline basis functions with strong vertical |None |kg m-2 s-1 |
| |correlation | | |
|Rain normalized number concentration, Nw|One value per profile |ln(8×106) for rain from melting |m-4 |
|(normalized_number _concentration) | |ice; ln(1.5×108) for warm rain; | |
| | |error 1.0 | |
|Melting layer |
|Melting-layer thickness scaling factor |One value per profile |Unity |None |
|(thickness_scaling) | | | |
|Aerosols | | | |
|Aerosol median volumetric diameter |One value per profile; prescribed if shortwave |None |m |
|(median_volume_diameter) |radiances not available | | |
|Aerosol number concentration |Cubic spline basis functions |ln(5×107); error 2.0 |m-3 |
|(number_concentration) | | | |
|Lidar backscatter-to-extinction ratio, S|Constant with height or cubic spline basis function |ln(0.05); error 1.0 |sr-1 |
|(lidar_backscatter _extinction_ratio) | | | |
3 Profile initialization
Once the look-up tables have been created, the program performs a retrieval on each profile in the file in turn. Before the retrieval of a particular profile can begin, the profile data need to be read in to memory and the relevant vectors and matrices constructed, indicated by Box 2 in Figure 2. This section describes the contents of these vectors and matrices, and the actions needed to populate them.
1 Contents of vectors and matrices describing retrieved quantities
In a particular profile, the target classification product is used to define which atmospheric constituents to retrieve and at what altitudes. In the most general case when ice, liquid, rain and aerosol are all present in the profile, the state vector is formed of a concatenation of sub-vectors for each constituent:
[pic].
The state vector for ice cloud is in turn a concatenation of the sub-vectors for each variable needed to describe the ice size distribution at each height:
[pic], where [pic], [pic] and [pic].
This matches the list of state variables given in Table 8 (except for riming factor for simplicity): [pic] is a vector of p variables describing the profile of the natural logarithm of ice extinction coefficient in the geometric optics approximation, [pic] is a vector of q variables describing the profile of the natural logarithm of a measure of number concentration N0’ = N0*/α0.6, where N0* is the normalized number concentration defined by Delanoe et al. (2005) defined in section 5.5.7.1.1. [pic] is a vector of r variables describing the natural logarithm of the ratio of lidar backscatter to extinction coefficient. The presence of the circumflexes above each symbol indicates that these vectors do not necessarily contain simply the values of these variables at each height; they may instead contain a reduced set of descriptors such as the control points of a set of cubic spline basis functions, or more simply just a single value that is constant for the entire profile. This is in order to improve efficiency by reducing the size of the state vector when a large number of retrieved degrees of freedom is not justified; further details are given in section 5.5.5.1.
The state vector for liquid cloud likewise is a concatenation of sub-vectors describing the profile of liquid water content L and total number concentration N:
[pic], where [pic] and [pic].
Similarly, the state vector for rain is a concatenation of sub-vectors describing the profile of rain rate R and normalized number concentration Nw.
[pic], where [pic] and [pic].
A similar vector may be defined for melting ice and for aerosols.
Each state variable has a corresponding prior estimate, written as a vector [pic], and these prior estimates have an associated error covariance matrix B. Appropriate values for the prior and its error are given in Table 8 and read in from the configuration file. Smoothness and flatness constraints may be applied via a Twomey-Tikhonov matrix T, defined in section 3.4.1. If x is of length n then both B and T are of size n×n. In practice the inverse of B appears in the equations, so this is stored (see section 8.3 for an explanation of why this is sparse). In fact, sub-matrices are stored for each constituent, i.e. for ice we would store B-1ice and Tice.
2 Contents of vectors and matrices describing observed quantities
The observations are read in to the observation vector, which is a concatenation of sub-vectors corresponding to the CPR, ATLID and the MSI:
[pic],
where the radar observation vector consists of a column of values of the logarithm of radar reflectivity factor at particular range gates, followed by a column of Doppler velocities at a subset of those gates where the signal was large enough for a reliable Doppler measurement. The lidar observation vector consists of a column of values of the logarithm of the Mie-channel apparent backscatter followed by a column of values of the logarithm of the Rayleigh-channel apparent backscatter. Finally, the imager observation vector is simply a list of radiances at different wavelengths:
[pic], [pic], and [pic].
The observations are assigned an error covariance matrix R, which is diagonal and each diagonal element is the square of the 1-sigma error reported in the input files. Since R-1 is used in the equations, this is what is actually stored.
Finally, we define yf as the forward modelled observations, but this is not filled until the forward model is run.
3 Populating the vectors and matrices
The actions to be carried out to populate the vectors and matrices defined in the two previous subsections are as follows:
1. Clear any data from the various objects that refers to the previous profile. The exception is that the retrievals in adjacent profiles will be correlated, so to speed up convergence it makes sense to use the solution from one profile as the first guess (but not the prior) at the next. This approach is now used in the Delanoe and Hogan (2010) ice-only retrieval in its application to A-Train data. Therefore some constituents may copy the previous solution to another variable at this point. Note that this does not imply that the cost function has an extra term favouring a solution close to the previous profile (as in, for example, Hogan 2007); to get that behaviour requires the Kalman smoother to be invoked specifically using the temporal_smoothing_error option in Table 3.
2. Read the properties describing the thermodynamic state of the atmosphere and other properties into the BackgroundConstituent object. Specifically, the satellite altitude is read along with the profile information: height, pressure, temperature, specific humidity and ozone mixing ratio. As indicated in the listing in the SUM, the names of these variables in the input file is specified in the configuration file.
3. Loop over each observation and read the current profile of observed variables into memory, together with its error and its detection variable (see Table 5). In the case of active instruments (radar or lidar), use the detection variable to determine which range-gates should be used. Firstly, all those gates where the detection variable is one of the values in detection_values (indicating that particulates are present) are added to the list. If the molecular_values list is present (for lidar variables able to detect Rayleigh scattering) then those gates where molecular scattering is present are added, but only if they lie beyond a gate where particulates are present. This way we only assimilate observations of molecular scattering where they can be used to provide information on the optical depth of a layer closer to the instrument. Thus, the following need to be stored for each variable of an active instrument:
• The observed data, yo, for that variable at the gates to be assimilated. In the case of backscatter variables, the natural logarithm is stored.
• The reciprocal of the error variance, R-1, for the corresponding observations (note that this is diagonal). So if the error in the ith element of yo is σi then the ith diagonal element of R-1 is 1/σi2. The forward model error is added at the same time; see Delanoe and Hogan (2010) for a discussion of how this is done mathematically for A-train observations.
• The indices to the gates of the full profile where the observations come from.
• The corresponding forward modelled data, yf, although at this stage in the algorithm it has not yet been used.
In the case of passive instruments, the process is simpler since the variables are not a function of range. In the case of radiances (or brightness temperatures) from the MSI, they will instead be functions of wavelength. This leads to the following vectors:
• The observed radiances, yo, at each wavelength.
• The reciprocal of the error variance, R-1, for these observations. This is diagonal if we assume the radiance errors to be uncorrelated. However, in the infrared, if there is an error in the atmospheric temperature estimated from ECMWF, then there will be an additional forward-model error that is correlated between channels, in which case R-1 would not be diagonal.
• The corresponding forward modelled data, yf, although at this stage in the algorithm it has not yet been used.
4. Loop over each constituent and allocate the necessary variables to represent it. The relevant classification variable and classification codes (see Table 2) are used to decide at which gates to retrieve that constituent. As shown in section 5.5.3.1 and Table 8, each constituent is described by several different variables (e.g. liquid water content and total number concentration in the case of liquid cloud), and for maximum flexibility, there are several different ways that the vertical distribution of each variable can be described. This description is specified in the configuration file as shown in Table 3. Thus for each variable [pic] of each constituent we need to store:
• The vector of state variables to retrieve [pic]. This can be initialized to the prior, or to the nearest solution value from the previous ray.
• The vector of prior estimates [pic]. This is set based on the values in the configuration file (see Table 3).
• The inverse of the error covariance for this prior estimate, Ba-1 .This is also set based on the values in the configuration file; Section 8.3 shows how to compute it for a given error variance and decorrelation length.
• A Twomey-Tikhonov matrix Ta, if a smoothness or flatness constraint is required. The structure of these matrices is given by (2) and (3), with the strength of the constraint governed by smoothness_factor and flatness_factor in the configuration file.
• A matrix W that will be used to convert [pic] to values at the range gates where the constituent is present via [pic]. Section 5.5.4.5 describes how it is formulated.
5. Allocate memory to store scattering properties. For each observation-constituent pair we will need to calculate the scattering properties needed to forward model that observation (given in Table 7) at each gate where the constituent is present. In addition, for each observation we will need to calculate a combined profile consisting of the same scattering properties but at all range gates in the atmosphere.
4 Minimizing the cost function
Now that the infrastructure for the retrieval has been set up, the cost function given by (1) can be minimized, which constitutes the core part of the algorithm depicted by Box 3 of Figure 1. The state vector that minimizes the cost function is denoted the solution. This is achieved using a packaged routine implementing the L-BFGS quasi-Newton algorithm (e.g. from the GNU Scientific Library), which was outlined in section 3.4.2 and Figure 1. This routine requires two functions to be provided, one which just computes the cost function for a given trial state vector, i.e. J = calc_cost_function(x) and the other which also computes the gradient of the cost function with respect to the trial state vector, i.e. [J, (xJ] = calc_cost_function_and_gradient(x). The steps that make up first of these functions are described in section 5.5.5. The second requires the “adjoint” of all these steps to be coded, which is a time consuming and error prone task if done by hand. Fortunately, a computational advance outlined in section 8.1 means that if we can write all the parts of the first function in C++ then it is possible for (xJ to be computed automatically. Therefore we do not spend any time here describing how to analytically differentiate an algorithm but assume that if code is available to compute simply the cost function, then its gradient can be computed automatically.
The minimizer is instructed to continue searching for the minimum unless one of the following conditions is met:
1. Convergence is achieved. This is deemed to have occurred occur if the L2 norm of the gradient of the cost function is less than a tolerance τ, i.e. [pic], where the summation is over each of the n state variables xi. The tolerance τ is set to the value of converged_gradient_norm in the configuration file, which would usually be set to 1.
2. The maximum number of iterations is reached. The maximum number of iterations is specified in the configuration file.
3. Algorithm fails to converge. This typically occurs quite close to the minimum of the cost function where the gradient of the cost function with respect to the state vector is quite small, and round-off error means that the gradient reported by the adjoint code is not quite accurate, such that if the minimizer tries to move in a downhill direction as determined by the gradient, it finds that the cost function actually increases.
When the minimizer stops searching, the convergence criterion is stored, errors are computed and the results are saved to the output file.
5 Calculation of the cost function
We now describe the steps needed to implement the function J = calc_cost_function(x) that is passed to the minimizer.
1. Compute the background scattering/absorption profile for each observational wavelength
o For each observation j, use the atmospheric thermodynamic profile compute the full profile of gas extinction coefficient αair, and interpolate it from the native model grid on to the Joint Standard Grid in the vertical. In the case of a lidar with a wavelength less than 1 micron, the density of the air is used to calculate the well-known Rayleigh profile of extinction coefficient. In the case of a radar with a wavelength larger than 1.3 mm (frequency less than 230 GHz), the Liebe (1985) model is used to calculate the gaseous extinction profile given the frequency, temperature, pressure and vapour pressure. In the case of wavelengths between 1 micron and 1.3 mm, it is assumed that the instrument is an imager/radiometer in which case the radiative transfer forward model needs to treat atmospheric absorption in a more complex way due to the spectral variation within the receiver frequency response window; therefore no action is taken. Since temperature and humidity are not retrieved, the gaseous scattering and absorption profile does not vary between iterations of the minimization, so after calculating the profile once for a particular profile, it is not recalculated in subsequent iterations.
2. Convert each constituent variable stored in the state vector to high resolution
o For each constituent i, convert the vectors describing each variable to high resolution via the operations [pic], [pic] and similarly for any other variables needed to describe the constituent. Here Wa,i is a weighting matrix described in section 5.5.5.3 that is computed at the profile initialization stage, and implements a particular representation such as cubic spline basis functions. In the case of ice, aice then represents the profile of the natural logarithm of extinction coefficient at the range gates at which ice is present.
o In the case of ice clouds, we convert the vector containing N0’ to one containing N0*; so if [pic] represents a vector containing entries ln(N0*) then we apply [pic].
3. Compute the profile of scattering properties for each constituent-observation pair
o For each observation j and each constituent i, use the look-up table specific to that pair to compute the scattering properties needed by observation j (Table 7) at the range gates where constituent i is present. Interpolation into the look-up table is as described in section 5.5.2.2.
o In the case of ice clouds observed by lidar, override the value of backscatter-to-extinction ratio returned from the look-up table with [pic] from the state vector. This is because the backscatter-to-extinction ratio is intrinsically variable in a way not captured in the look-up tables, and so is best to be retrieved.
4. Merge the scattering property profiles
o For each observation j, merge the scattering properties associated with each constituent to yield a single combined profile of scattering properties. The mathematical rules for how to merge each scattering variable are given in section 5.5.5.3.
5. Run the radiative transfer forward models
o For each observation j, run the appropriate radiative transfer forward model using the combined profile of scattering properties as input. The status of these models is described in section 5.5.8. In the case of radar and lidar we use the Multiscatter package (Hogan 2008, Hogan and Battaglia 2008). The internal details of this need not be described, but for our purposes we treat it as a function that is called as follows:
i. For radar: βradar = Mradar(α, s,ω, g, αair), and radar reflectivity factor is then computed from apparent backscatter coefficient using [pic] (Hogan and Battaglia 2008). Note that currently the Multiscatter package does not compute Doppler velocity and the effect of multiple scattering on it, so we simply report the Doppler velocity from the combined scattering profile.
ii. For lidar: [βmie, βray]= Mlidar(α, s,ω, g, ra, f, αair).
o In the case of radar and lidar modelling, it is inefficient to run a multiple scattering model on clear-sky gates between the instrument and the first particulates, since agreement or otherwise between the forward modelled and observed backscatters cannot contribute to the retrieval. It is valid to assume that the multiple scattering contribution from clear air is negligible, in which case we apply an optimization in which the optical depth δtop of the atmosphere from the highest gate n to the first particulate gate np is computed from [pic], where Δr is the range gate spacing. The forward model is run only on gates between the first particulate gate and the last gate where a signal is detected by the instrument, and then the backscatter profile is multiplied by exp(2τtop) to account for two-way transmission between the instrument and the first particulate gate.
o In the case of imager channels, we use the Delanoe and Hogan (2008) model in the longwave and the LIDORT model in the shortwave to compute the radiances I as follows:
i. In the infrared separately in each channel: IIR=MIR(α, ω, g, thermodynamic profile, ozone profile, CO2 profile, surface temperature & emissivity, satellite zenith angle).
ii. In the shortwave separately in each channel: ISW=MSW(α, ω, g, thermodynamic profile, ozone profile, surface albedo, satellite and solar zenith angle, azimuth between sun and satellite).
6. Compute the observational contribution to the cost function
o For each observation j, compute the contribution to the cost function [pic], where δyj = yjo – yjf. Sum the result to obtain the total contribution to the cost function from the observations.
7. Compute the constituent contribution to the cost function
o For each constituent i, compute the contribution to the cost function [pic], where δxi = xip – xi. Since B-1 is at most tridiagonal (see section 8.3) and T is pentadiagonal, this calculation is very efficient. Add any additional terms associated with prior constraints, specifically the gradient constraint applied to liquid water content (section 5.5.7.2.2) and the continuity constraint applied between ice flux just above the melting layer and rain rate just below (section 5.5.7.4). Sum the result to obtain the total contribution to the cost function from the constituents.
The sum of the cost function from the observations and the constituents then provides the total cost function, which is returned to the minimizer.
Steps 3 to 5 of this sequence constitute the forward model, depicted schematically in Figure 4. The coloured boxes on the top line indicate the state variables, which describe (in this example) the separate constituents of ice cloud (including snow), liquid cloud, rain and aerosol. The first yellow box indicates the use of look-up tables to calculate profiles of scattering properties (extinction coefficient, asymmetry factor etc) for each of the observations (radar, lidar and imager in this example). This yields separate profiles of scattering properties for each observation-constituent combination. Note that the radar-aerosol combinations is omitted since the radar sensitivity to small particles is much too low to detect aerosols. The second yellow box represents the act of combining these variables into a single combined profile of scattering properties for each observation. The final step in the forward model is the radiative transfer part, where the various radiative transfer schemes are run on the combined scattering property profiles to yield forward modelled estimates of the measurements, yf.
1 Different approaches to representing profiles of variables
Most retrieved atmospheric constituents are described by two variables (e.g. a measure of mass per unit volume and a measure of number concentration), implying that two variables need to be retrieved per range gate at which the constituent is present. Frequently the observations do not provide high-resolution information on the vertical structure of both variables, and it makes sense for efficiency and robustness to describe the profile with fewer variables, such as the coefficients of a set of cubic-spline basis functions, or to describe a profile with a simple linear slope, its mean and its gradient. For flexibility, it is also useful to be able to switch easily between different representations of the profile. A generic method for representing profiles in which the m-point high-resolution profile of a variable [pic] is derived from the n state variables describing it [pic] (note that this is only the part of the full state vector describing the profile of one variable) is via a matrix-vector multiplication of the form [pic], where W is an m(n matrix describing the representation.
The representation for a particular state variable is defined in the configuration file, as outlined in Table 3. Each representation is implemented in terms of the matrix W as follows:
o Predefined: If a is not retrieved then it is simply set to the prior value and W is not used.
o Constant: In the simple case when a single value is retrieved that is constant with height then [pic] is simply one value and W is an m(1 vector of ones. An example is the representation of total droplet number concentration, a variable that is assumed constant with height in liquid clouds.
o Linear: The next level of complexity is to allow the variable to vary linearly with height, in which case element i of a would be given by [pic], where zi is the height of the ith range-gate where the constituent is present, z0 is the average height where the constituent is present, and [pic] and [pic] are the two state variables representing the mean and the gradient of the variable. In this case W is an m(2 matrix whose first column consists of ones and whose value at the ith row second column contains the value zi – z0. An example is one of the options considered by Delanoe and Hogan (2010) for representing the variation of the natural logarithm of ice lidar ratio with height when no HSRL was available.
o Direct: If we want the state variable to correspond to the values at adjacent range gates, then W is simply an n(n identity matrix, although in practice we would simply apply [pic].
o Cubic spline: In this case W describes the cubic-spline basis functions, exactly as described by Hogan (2007).
2 Merging the scattering profiles
To merge the information on scattering properties from each constituent into the combined profile (step 4 in section 5.5.5), the observation object loops through each of the required properties and calculates the direct sum or the weighted sum of the contribution to that property from each constituent. The properties listed in Table 3 are summed as follows:
• Extinction coefficient, backscatter coefficient and scattering coefficient: direct summation over N constituents, so the combined extinction coefficient is [pic], and similarly for backscatter and scattering coefficient.
• Backscatter-to-extinction ratio and single-scattering albedo: summation weighted by extinction coefficient, which must therefore also be present as a requested property; so the combined backscatter-to-extinction ratio is given by [pic].
• Asymmetry factor: summation weighted by the scattering coefficient; therefore either the scattering coefficient must be present, or both the extinction coefficient and the single-scattering albedo ((: [pic].
• Mie droplet fraction f: The value of this variable is unity for liquid cloud and rain constituents observed at a wavelength of 2 microns or less, and zero in all other situations. The summation is weighted by backscatter coefficient; therefore, either the backscatter coefficient must be present, or both the extinction coefficient and the backscatter-to-extinction ratio. [pic]
• Doppler velocity vd: summation weighted by the backscatter coefficient; therefore, either the backscatter coefficient must be present, or both the extinction coefficient and the backscatter-to-extinction ratio. Note that if vertical wind w is non-zero then the contribution to the Doppler velocity must be treated differently: if the backscatter of the air at that wavelength is non-zero then it will contribute to the weighted summation in the usual way except with a zero Doppler velocity. The weighted summation will then contain the Doppler velocity relative to the air, i.e. the backscatter-weighted terminal fall speed vt. Then the air motion is simply added to this. [pic].
• Equivalent-area radius: the extinction-weighted average of the square of the equivalent-area radius from each constituent is computed, and the square-root of the result is taken to yield the combined value: [pic]
|[pic] |
| |
|Figure 4: Schematic of the operation of the forward model, which takes as input the state vector x and outputs the |
|forward modelled observations yf=H(x). |
6 Calculation of the retrieval error
1 Calculation of the solution error covariance matrix
We have now completed our description of the minimization process depicted by Box 3 in Figure 1. The packaged routine returns the state vector at the minimum of the cost function, or the “solution”. The final step is to estimate the error in the retrieval. This may be achieved by computing the Hessian matrix [pic] and inverting it to yield the error covariance matrix of the solution [pic]. Often we are interested only in the error standard deviation of each retrieved value, which is the square root of the diagonal of the error covariance matrix. In the quasi-Newton minimization method, the Hessian matrix is not normally calculated, but the automatic differentiation approach described in section 8.1 allows us to compute the full Jacobian matrix H quite efficiently, and we may then use (6) to obtain the Hessian.
2 Transformation of error covariance matrix
Assuming we can compute the error covariance matrix, some further calculations are required to compute error statistics that are useful for users. This is because the errors are in terms of state variables which may not make much sense on their own, such as our N0’ quantity for ice clouds (see section 5.5.7.1) which is used because it has a good temperature-dependent prior, but is otherwise not directly useful. The other factor is that often a reduced description of the profile of state variables is used internally, such as the coefficients for a set of cubic spline basis functions, in which case the errors would be in terms of these coefficients. A further difficulty with the raw error covariance matrix is that it is large (the length of the state vector squared) and its size varies from ray to ray which would necessitate a rather complex format if it were to be stored in its entirety. Therefore 1D error descriptors are derived from the 2D matrix.
Let us define the vector m as a vector of quantities describing a particular retrieved constituent for which we want error estimates. For example, for ice clouds this could be the logarithm of the ice water content and the effective radius at each range gate at which ice clouds are retrieved:
[pic].
These variables are either extensive or intensive, and the same look-up table code used to compute the scattering properties (described in section 5.5.2.2) can be extended a little to also compute these microphysical variables for the same inputs. Therefore, the necessary parts of the code to implement the function m = F(x) can be easily identified and used to compute m. The automatic differentiation is flexible enough that for any parts of the code for which we can define inputs (x) and outputs (m), we can compute the corresponding Jacobian matrix M=∂m/∂x. From this, the error covariance matrix of m, written as Sm, is derived from the error covariance matrix of x as follows:
[pic].
Appendix 2 of Delanoe and Hogan (2008) describes how to do this the “hard way”, i.e. manually differentiate all the steps between x and m.
3 One-dimensional error descriptors derived from the error covariance matrix
The main error descriptors reported are the RMS errors, which are simply the square-roots of the diagonals of the transformed error covariance matrices described in section 5.5.6.2. Two other descriptors are provided. The first is the error correlation between two variables describing a particular constituent at a particular range gate, such as the correlation of the error in lnIWC and in re. This is calculated as follows. Suppose element Sij contains the error covariance of lnIWC and in re at a particular range gate, then Sij/(SiiSjj)1/2 is the corresponding error correlation, having a value between –1 and 1. The second is the width of the diagonal band in the error covariance matrix for a particular variable, indicating the extent to which errors are correlated with adjacent values. For element i, this width is given by Equation 17 of Pounder et al. (2012):
[pic].
These are stored in the output file with the suffix _error_corr_scale, as shown in section 5.3.2, and have the units of metres.
4 Error descriptors derived from the averaging kernel
The averaging kernel is a matrix A associated with the retrieved variables x that describes their “information content” and is given by
[pic],
where Sx is the error covariance matrix of x, R is the error covariance matrix of the observations y and H is the Jacobian matrix.
The significance of the averaging kernel is that it can be thought of as the rate of change of the retrieved variables with respect to the “true” value of each state variable, and so gives an indication of whether the solution is contaminated with a priori information. Thus the identity matrix would indicate a retrieval for which the information derives entirely from the observations. We follow Pounder et al. (2012) and summarise the p×p averaging kernel for a particular p×1 vector of derived quantities (e.g. lnIWC) by two one-dimensional descriptors. This is done by recognising that the row i of an averaging kernel (which we write as aiT) is a smoothing function that describes how retrieved element i depends on the true values at each other point. The first descriptor is aiTu (where u is a column vector of unit elements) and can be considered as a rough measure of the fraction of the retrieval that comes from the observations rather than the prior or the additional constraints. These are stored in the output file with the suffix _avgkern_sum, as indicated in section 5.3.2. The second descriptor describes the width of the diagonal band of the averaging kernel and gives a measure of the degree to which the true information is smoothed by the retrieval. This is given the suffix _avgkern_corr_scale, and is calculated in the same way as the correlation scale of the error covariance matrix described in the previous subsection. It has yet to be seen how these two correlation scales differ from each other.
7 Retrieved atmospheric constituents
So far the description of the algorithm has been mostly in terms of the generic aspects that are common to all observations and all constituents. Naturally, reliable retrievals require that suitable assumptions and state variables are chosen for each constituent, and that each instrument is forward modelled appropriately. In this section we describe the representation of the four retrieved constituents. It should be noted that the flexibility offered by the generic part of the code means that alternative representations can quite easily be tested. Table 8 lists the state variables use to describe each retrieved component. These are discussed much more fully and justified in the following subsections.
1 Ice cloud
The treatment of ice clouds is largely the same as that described by Delanoë and Hogan (2008b, 2010) and in [CASPER-ATBD], and has now been applied successfully to three years of CloudSat and Calipso data. However, it is necessary to ensure that the algorithm represents all ice phase hydrometeors above the melting layer, including snow and graupel. It was suggested by Hogan et al. (2001) that from cloud radar observations there was no sharp distinction between ice clouds and snow. Therefore, rather than introducing snow as a new component entirely, it makes more sense to ensure that the ice cloud retrieval at warmer temperatures and large ice water contents behaves as a snow retrieval.
1 State variables
To describe the ice microphysics, three or four state variables are required at each gate, as shown in Table 8 (along with the state variables for the other species to be retrieved, as described in the following sections). It should be stressed that these are not necessarily the same as those used to describe the cloud in the output product, but rather are chosen because they may have a good a priori estimate. In the case of ice clouds, for example, the main output variables would be ice water content, effective radius, visible extinction coefficient and equivalent snowfall rate.
The primary variable describing ice clouds in the state vector is the ice visible extinction coefficient, αv, a key property describing the radiative properties of the cloud. This is the most logical choice also because it is what is most naturally retrieved by the HSRL lidar.
|[pic] |Figure 5. (a) The temperature dependence of N0* for|
| |each size distribution within the large in situ |
| |database of Delanoë et al. (2005) (dots), |
| |superimposed by the mean N0* in 5[pic]C temperature|
| |ranges and various ranges of ice water content IWC |
| |(lines and symbols). (b) The same but for the |
| |variable N0’ = N0*/αv0.6. From Delanoë and Hogan |
| |(2008). |
In order to relate αv to other moments of the size distribution such as radar reflectivity factor (Z) or ice water content (IWC), it needs to be supplemented in the state vector by another intensive or extensive variable, such as a measure of particle size or number concentration. This additional variable should ideally have two key properties. Firstly, a good a priori estimate of it should be available as a function of temperature. This ensures that in regions where only the radar or the lidar is available, the scheme will tend towards existing empirical relationships involving temperature, such as the formulae for IWC as a function of Z and temperature (e.g. Liu and Illingworth 2000, Hogan et al. 2006a, Protat et al. 2007). It was demonstrated by Hogan et al. (2006a) that the temperature dependence in these relationships must arise via the temperature dependence of the number concentration parameter of a size distribution, commonly referred to as N0. Secondly, it should be easy to combine this additional variable with αv to estimate any other property of the size distribution. A good candidate is the ice “normalized number concentration parameter”, N0*. This is defined in section 5.5.3.2, but for a full description of the properties of this variable (including to what it is normalized), the reader is referred to Delanoë et al. (2005). but for our purposes, the key property that we exploit is that for any intensive variable y and extensive variable Y there is a near-unique relationship between the ratio αv / N0* and both y and the ratio Y / N0*. However, rather than using N0* directly in the state vector, it turns out to be better to use the quantity N0’= N0*/αv0.6. In Figure 6 it is shown that this mathematical combination of variables has the useful property of being dependent on temperature but independent of ice water content; it therefore has a robust prior estimate that can be used (defined in Table 8).
We follow Delanoe and Hogan (2010) and represent N0’ by a reduced set of cubic-spline basis functions, resulting in the retrieved N0’ being continuous in itself and its first and second derivatives. This also ensures a shorter computation time. However the radiative transfer forward model works on the lidar range grid, so at the beginning of each iteration, the amplitudes of the basis functions within the state vector, xn, have to be converted to a vector n containing the values of N0’. This transformation is achieved via [pic], where the elements of the matrix containing the basis functions are set as described in the appendix of Hogan (2007). The equivalent adjoint statement to this is [pic]. Section 5.5.4.3 discusses different approaches to reducing the size of the state vector while still being able to retrieve the main structure in the profiles of cloud properties.
The third state variable for ice clouds is the lidar extinction-to-backscatter, S. Even though the HSRL Rayleigh channel provides a “direct” measurement of extinction, it is often quite noisy, and therefore it is beneficial to make use of the HSRL Mie channel, which while less directly dependent on extinction, will be less noisy in cloud. By retrieving S but with a constraint on its smoothness with height, it is possible to retrieve a more accurate extinction than would be possible with either HSRL channel on their own.
The final variable, the “riming factor”, is intended to represent the transition from low density ice aggregates in cirrus clouds to high density ice (e.g. hail and graupel) in convective and some intense stratiform precipitating systems. It is also be useful for modelling snow. This avoids the need to have independent species for graupel and snow that would coexist with smaller ice particles. The riming factor scales-up the density of all the ice particles in the profile for any cloud extending up continuously from the 0°C level (but of course not allowing the density of any individual particle to exceed that of solid ice). It essentially represents the fact that when supercooled liquid is present in sufficient concentrations, ice particles can grow by riming. The information for it would come from the Doppler velocity in stratiform clouds (high density particles fall faster for a given ice water content), the path-integrated attenuation (high density ice can attenuate the radar signal) and from the multiple scattering signal (the most extreme examples of multiple scattering in CloudSat are due to high density ice high above the melting layer). However, the difficulty in interpreting these measurements means that it is not justified to include a very complicated representation, and so initially only a single variable will be used to scale the density through the ice above the rain.
Figure 6 illustrates how Doppler velocity provides information on the riming factor. The black line shows the default un-rimed relationship between particle size and ice fraction (i.e. density), with the colours beneath indicating that generally these particles have a fall speed of 0.5-1 m s-1. The riming factor scales up the ice fraction, leading to larger fall velocities for a given size. Therefore, if the Doppler velocity is larger than would normally be expected given the inferred size of the particles, then the algorithm will retrieve a riming factor greater than 1.
2 Microphysical and scattering assumptions
In order to forward model the various observations from these state variables it is necessary to make several assumptions for ice clouds:
• The shape of the size distribution has a particular form; we use the unified size distribution of Delanoë et al. (2005), which has been found to fit in-situ aircraft observations, but the Field et al. (2005) distribution is also possible. This is listed along with the size distribution assumption of the other retrieved constituents in Table 9.
• The mass-size (or density relationship must be specified. We use the Brown and Francis (1995) relationship, which was found by Hogan et al. (2006a) to produce excellent agreement between radar reflectivities calculated from aircraft probes and simultaneously measured by radar. When the riming factor is included as a state variable, then it is envisaged that it would scale the density such that 0 corresponds to the Brown and Francis (1995) value at a given size, and 1 corresponds to the value for solid ice.
|[pic] |
|Figure 6. The colours represent the fall speed of ice particles in m s-1, as a function of their diameter (maximum |
|dimension) and the fraction of the volume of the smallest sphere that encloses the particle that is ice, according to |
|the model of Heymsfield and Westbrook (2010). The black line depicts the relationship between ice fraction and diameter|
|of Brown and Francis (1995). |
• The backscattering of an ice particle at 94-GHz is described by the self-similar Rayleigh Gans method of Hogan and Westbrook (2013), applicable to ice aggregates both smaller and larger than the wavelength. This is coupled to the evidence presented by Hogan et al. (2012) that aggregates (which dominate ice clouds and snow) have an aspect ratio of around 0.6 and fall with their longest dimension in the horizontal (although the Hogan and Westbrook scattering model is much better than that of Hogan et al. for particles larger than the wavelength. The scattering model along with those for other retrieved constituents is listed in Table 10.
• The extinction efficiency of the ice particles at lidar and imager wavelengths may be derived from Baran’s (2004) ice aggregate database. In the case of lidar, this would be an improvement over [CASPER-ATBD], in which it was assumed that the geometric optics approximation holds and the extinction efficiency of all ice particles is exactly two, meaning that the extinction cross-section of a particle is twice its geometric cross-sectional area. In practice this is a good assumption because the particles are much larger than the wavelength. We also use Baran’s (2004) phase functions for the scattering asymmetry factor, which is required by the wide-angle multiple scattering model. Baran’s phase functions are preferred over those of Yang et al. (2000), since the former are supported by aircraft observations of scattered sunlight while the latter are dominated by pristine crystal types when large in-situ aircraft databases have shown that irregular particles (assumed by Baran) dominate.
• We assume that the extinction-to-backscatter ratio is not something well constrained by theoretically calculated phase function, but rather that it is something that needs to be retrieved (and indeed its retrieval is straightforward with an HSRL). Indeed, it was shown by Hogan (2008) that the phase functions of the Baran (2004) library (principally ice aggregates) and the widely used Yang et al. (2000) library (principally idealized particle shapes such as columns and dendrites) are very different. The main difference is in the backscatter direction, where Baran’s phase functions are flat and Yang’s upturn significantly, leading to a range in extinction-to-backscatter ratio of a factor of 8.
In the extension of this approach to include all ice hydrometeors including graupel and snow, the validity of some of these assumptions, found to be good for ice clouds, may need to be re-examined.
|Table 9. List of assumptions on the shape of the size distribution for each particle type, together with the reference. |
|Constituent |Size distribution type |Reference |
|Ice cloud |Delanoe |Delanoe et al. (2005) |
|Liquid cloud |Lognormal with shape factor 0.38; |Miles et al (2000); |
| |or Weibull |Liu and Hallett (1997) |
|Rain |Gamma with shape factor 5.0 |Illingworth & Blackman (2002) |
|Aerosol |Lognormal with shape factor 2.0 |Hesse et al. (1998) |
|Table 10. Summary of the particle scattering models used for each different particle type (rows) and each part of the spectrum |
|(columns), together with the model used to estimate particle terminal fall speed for the Doppler radar. Note that the thermal |
|infrared, solar radiances and UV lidar all use the same scattering models but with the calculation performed for the required |
|wavelength. |
|Particle type |Radar (3.2 mm) |Radar Doppler |Thermal IR, Solar, UV |
|Aerosol |[Aerosol not detected by radar]|[Aerosol not detected by radar] |Mie theory, Highwood refractive index |
|Liquid droplet |Mie theory |Beard (1976) |Mie theory |
|Rain drop |T-matrix: Brandes et al. (2002)|Beard (1976) |Mie theory |
| |shapes | | |
|Ice cloud particle |Self-similar Rayleigh Gans |Heymsfield & Westbrook (2010) |Baran (2004) |
| |(Hogan and Westbrook 2013) | | |
|Graupel and hail |Mie theory |Heymsfield & Westbrook (2010) |Mie theory |
|Melting ice |Wu & Wang (1991) |TBD |Mie theory [Or ignore since these |
| | | |observations virtually insensitive to |
| | | |melting ice given ice cloud also in the|
| | | |profile] |
3 Liquid cloud
The key variables reported from the liquid-cloud component of the retrieval are the liquid water content, the effective radius, the visible extinction coefficient and the total droplet number concentration, although only two of these four variables are independent. These are listed as liquid-cloud output variables in the [PDD] along with visible optical depth, which is simply the vertical integral of extinction coefficient. As discussed in the [PARD], even with the many instruments available on EarthCARE, the retrieval of liquid cloud microphysical properties represents a significant challenge, because the profile of cloud variables can be inaccessible to the measurements, particularly towards cloud base where the lidar has lost signal. Therefore, extra constraints on the profile are required for those profiles where the active instruments provide less information.
1 State variables
The state vector contains the liquid water content (LWC) at each range-gate containing liquid, and a single value for the total number concentration for each contiguous liquid layer (see Table 8). LWC is chosen as a state variable because (1) the radar path-integrated attenuation (PIA) is proportional to it, and (2), its rate of increase with height rarely exceeds the known adiabatic value, potentially enabling a gradient constraint (described below) to be applied and hence lead to more physically realistic profiles towards cloud base where there is little information from the observations. An optional alternative to the gradient constraint is to use a very low prior LWC value together with a smoothness constraint in height; this ensures that where the radar and/or lidar detect the cloud well, the liquid water content is well constrained by the observations, but in other regions (such as where the lidar signal is largely extinguished), the liquid water content tends smoothly (but not necessarily adiabatically) to the low prior value. Either of these approaches allows path constraints, such as from the radar path-integrated attenuation and the shortwave radiances to add or remove liquid from the parts of the profile most weakly constrained by the active sensors.
The other state variable describing liquid clouds is the cloud-base droplet number concentration. The primary source of information for this variable is the MSI shortwave radiances, which provide information on cloud-top droplet size, which is related to the number concentration via the liquid water content. Since there is no real information on the vertical profile of size or number concentration, it is only valid to retrieve one value per continuous cloud layer. In perfectly adiabatic clouds, number concentration would be expected to be constant with height. The a-priori value for number concentration in boundary-layer clouds depends on whether the cloud is a marine or continental airmass; since this is difficult to diagnose in real time for any particular cloud, we take the simple approximation of using different values over land and sea, taken from Miles et al. (2000) and shown in Table 9.
2 Microphysical and scattering assumptions, and additional constraints
The following assumptions must be made in order that the observations can be forward modelled:
• The size spectra of liquid water droplets are assumed to follow a lognormal distribution with a shape parameter of 0.38, found to be a good fit to both marine and continental clouds (Miles et al. 2000). An alternative to consider is the Weibull distribution; this form has a sound theoretical basis (Liu and Hallett 1997), and fits the observational data presented by Martin et al. (1994).
• The lidar extinction-to-backscatter ratio for liquid droplets is taken to be constant at 18.5 sr (O’Connor et al. 2004).
|[pic] |Figure 7. Aircraft measurements of the |
| |vertical profile of liquid water content and |
| |stratocumulus clouds (solid lines), together |
| |with the adiabatic profile calculated from the|
| |cloud-base temperature and pressure (Slingo et|
| |al. 1982). |
A potentially problematic area is that often it is not possible to locate cloud base unambiguously, since the lidar will be rapidly extinguished and the radar return may be dominated by drizzle falling from cloud base. The target classification product provides an estimate of the location of cloud base, but this may not be accurate. Our approach is to retrieve liquid cloud properties not only in the gates indicated by the target classification, but also in an additional 300 m of gates below the lowest liquid cloud base. A one-sided gradient constraint is then used for each liquid layer in the profile to ensure that in the absence of information on the shape of the profile from the observations, the liquid water content (LWC) will fall smoothly to zero towards cloud base at a rate expected of an adiabatic cloud (supporting measurements shown in Figure 8). This is implemented by adding an additional term to the cost function of the form:
[pic],
where
[pic].
Note that the summation in the definition of Jgrad only includes a gate i if liquid cloud is present in both gates i and i+1. In this equation, ci+1/2 is the adiabatic LWC gradient between gates i and i+1. The term λ weights this constraint against the other terms in the cost function. Since we wish to prevent gradients that are steeper than adiabatic, we can use a large value, and satisfactory results have been found with λ =109.
4 Rain
The retrieval of liquid precipitation is potentially the most difficult component of the “Best Estimate” algorithm, because most of the information will be from the radar, and its high wavelength is optimized for clouds rather than precipitation. Moreover, the microphysics of precipitating systems can be very complex, and dependence on prior information will be necessary when this information cannot be determined from observations.
1 State variables
Two components of the state vector are used to represent rain, as listed in Table 9. The primary component is the rain rate. This variable has the advantage that the fast fall-speed of raindrops means that it varies only slowly with height. Therefore it may be represented adequately with a reduced set of cubic-spline basis functions in those regions of the profile identified as rain or drizzle by the Target Classification, plus a constraint on its “flatness” with height, as described in section 3.4.1. For “warm rain” (e.g. generated within stratocumulus clouds), the flatness or smoothness constraint will be weaker due to the greater rate of evaporation beneath cloud base by the smaller drops. For rain originating from melting ice, variables describing the melting layer are introduced as outlined in section 5.5.7.4.
The second variable is the normalized number concentration parameter, Nw (Bringi and Chandrasekhar 2001, Testud et al. 2001, Illingworth and Blackman 2002). This plays exactly the same role as N0* in ice clouds (Delanoë et al. 2005, Delanoë and Hogan 2008), and moreover it is reasonable to assume that it is constant with height in rain. Therefore a single value would be retrieved for every profile containing rain. The information on this variable would come in from the Doppler velocity, and from the consistency of the gradient in reflectivity in rain and the path-integrated radar attenuation. The “Marshall-Palmer” value of 8000 m-3 mm-1 may be used as the prior constraint and is consistent with recent work (Illingworth and Blackman 2002). In principle this should be correlated with the N0* of the ice above the melting layer, but there have been suggestions that aggregation and break-up means that particle number flux density is not conserved across the melting layer, although other studies suggest that this is a weak effect (Fabry and Zawadski 1995). Moreover, the N0* of ice clouds originates primarily from the radar-lidar combination but the lidar rarely penetrates down to the melting layer and so the N0* here will be almost entirely from the prior.
2 Microphysical and scattering assumptions, and additional constraints
The following assumptions must be made in the formulation of the rain component of the algorithm:
• A shape for the size distribution of raindrops. Initially we propose to use a gamma distribution with a shape factor of 5 (Illingworth and Blackman 2002), but with the flexibility to try other values.
• A relationship between the aspect ratio of oblate raindrops and their size; we will use the recent formula of Brandes et al. (2002).
6 Melting layer
The final retrieved constituent is the melting layer. There are two approaches that could be taken:
• The simplest approach is to treat the melting layer only for radar attenuation, neglecting its contribution to radar backscatter, lidar backscatter/attenuation and the optical properties that are required in radiance calculations in the shortwave and infrared. This can be justified in the case of spaceborne observations since (1) at 94-GHz there is no high reflectivity “bright band” at the melting layer, (2) the lidar has almost always been completely extinguished before penetrating down to the melting layer, (3) infrared and solar radiances are dominated by ice and liquid clouds in the profile, not by a thin layer of melting ice. In this approximation, an additional attenuation is added to the highest range gate containing rain.
• The complete approach employs a combination of a microphysical model for the characteristics of the size distribution of particles and their liquid-to-ice fraction, both as a function of depth into the melting layer. This is then coupled to scattering models to describe the internal structure of melting particles in a way that can be used in scattering calculations, such as the stratified sphere model of Wu and Wang (1991).
Here we will consider the first case. It was reported by Matrosov (2008) that at 94 GHz, there is not a great deal of sensitivity to assumptions about particle morphology, and from his Figure 11 (covering rain rates up to 14 mm h-1), we see that to within 5% the two-way attenuation through the melting layer (in dB) is 2.2 times the rain rate (in mm h-1). Nonetheless, Matrosov (2008) acknowledged that if the temperature gradient of the atmosphere in the vicinity of 0°C is such that the melting layer is thicker or thinner than average, then an error will result. Since these subtle properties of the atmospheric temperature structure are unlikely to be captured in the ECMWF analysis, we instead choose to retrieve the natural logarithm of a multiplying factor on the thickness of the melting layer, which acts to scale the attenuation up or down. The natural logarithm ensures that the scaling factor does not go negative. Information on the scaling factor comes from the profile of reflectivity factor coupled to the path-integrated attenuation from the radar surface return over the ocean.
To implement the constraint that the precipitation rate should be conserved across the melting layer, we add the additional term to the cost function that is only used when a melting layer is present:
[pic],
where I is the ice flux just above the melting layer and R is the rain rate beneath. The λ factor scales the importance of this term relative to others in the cost function. Since we have not yet tested this term, the optimum value of λ has yet to be determined. The use of natural logarithms ensures that it is the fractional difference in precipitation rate that is penalized, and is also the most natural since the natural logarithms are the quantities that are actually held in the cost function.
7 Aerosol
Compared to liquid clouds and rain, the retrieval of aerosols is relatively straightforward in principle: the HSRL will provide detailed information on the profile, with less ambiguity than an ordinary backscatter lidar due to the more direct extinction information. Two additional constraints on optical depth are available to be used: (1) over the ocean the lidar surface backscatter provides is stable enough to infer the path-integrated attenuation, which in this case is the two-way optical depth at the lidar wavelength (Josset et al. 2008); (2) in the daytime the shortwave MSI radiances will provide an estimate of aerosol optical depth, being more accurate over darker surfaces such as the ocean. The approach is somewhat similar to the way that optically thin ice clouds are retrieved when only the HSRL lidar is available. However, in practice aerosol information derived purely from the HSRL on a profile-by-profile basis is likely to be rather noisy due to the fact that aerosol is quite weakly attenuating compared to cloud, and the aerosol will provide only a weak target leading to quite a noisy Mie channel. The ATLID-only aerosol extinction algorithm will use an adaptive horizontal smoothing approach to overcome this problem, but such an approach is not directly compatible with the profile-by-profile approach needed to retrieve clouds, which have a much stronger degree of horizontal variability. Therefore, three possible approaches can be conceived at this stage for the aerosol:
1. Aerosol is retrieved on a profile-by-profile basis as part of the Best Estimate algorithm, but incorporating a Kalman smoother to enact horizontal smoothing. This should yield similar results to the adaptive smoothing approach of the ATLID-only algorithm, but is an approach compatible with the profile-by-profile variational retrieval. It is described in Chapter 7 of Rodgers (2000) and essentially uses the retrieval of properties in one profile as a constraint on the next (using a term in the cost function that looks like an a-priori constraint). In order that the smoothing is symmetric in time, a reverse pass of the algorithm is required in which each profile is re-retrieved but this time constrained by both the profile after and the profile before.
2. The smooth ATLID-only aerosol retrieval is imported into the Best Estimate algorithm and output unchanged. It might be used in the shortwave and infrared forward models used in the retrievals of clouds in the same profile.
3. The smooth ATLID-only aerosol retrieval is imported into the Best Estimate algorithm and used as a prior constraint on the retrieval performed within the Best Estimate algorithm itself. A variant on this approach would be to apply the Kalman smoother as well, but this would be applying smoothing in two different ways and so is conceptually not very appealing.
So far, only (1) above has been implemented, but simplified in that the Kalman smoother only operates in the forward sense; the reverse pass has yet to be implemented. The aerosol retrieval works with both simple backscatter or the two HSRL channels, but solar radiances have yet to be implemented.
– a forward-only Kalman smoother In terms of development, the priority is to add an aerosol retrieval (option 1 above) and test its performance before making the additional (smaller) modifications to also enable options 2 or 4. We propose three levels of complexity that would be implemented in turn:
1. The simplest aerosol retrieval would rely purely on the HSRL information. Extinction and lidar ratio at the wavelength of the lidar would be retrieved. This approach requires no microphysical information at all on the nature of the aerosol particles other than the a-priori value for the lidar ratio of the aerosol particles.
2. The next step would be to incorporate the aerosol into the forward modelling of MSI solar radiances. This would require an assumed aerosol size distribution and type, in order to provide the scattering properties at each MSI wavelength in order that the MSI radiances can be forward modelled. However, the wavelength-dependent information would not be used to retrieve aerosol type or size.
3. The most sophisticated retrieval would add an extra term to the state vector, most likely mean aerosol size, in order that the wavelength-dependence of the aerosol scattering inferred from the various MSI solar radiances can be exploited. This would be the physical equivalent of the Angstroem exponent.
These will only be implemented in VARSY if time permits.
1 State variables
The main state variable is aerosol total number concentration, which is retrieved independently at each lidar range gate (see Table 8). A weak smoothness constraint is used to ensure lidar instrument noise does not feed through into the retrievals. The reason that number concentration is used rather than, for example, the geometric-optics extinction coefficient (i.e. twice the sum of the geometric area of the particles in a unit volume) is that if mean size is held constant then the extensive scattering properties such as extinction and backscatter coefficients vary linearly with number concentration, providing a unique route from the backscatter and extinction observations to the state variables. However, if number concentration is held constant while geometric extinction is varied, then this implies a mean size variation, and since the wavelength is of order the size of the particles, the relationship between geometric extinction and extinction/backscatter at the wavelength of the lidar is not necessarily monotonic. This leads to multiple minima in the cost function and the retrieval does not converge correctly.
In terms of exploiting the HSRL information in aerosol conditions that are not too optically thin, the other key state variable is the lidar backscatter-to-extinction ratio. This will most likely be taken to be constant with height in each layer, although at night when there is no solar background in the HSRL, it may be possible to estimate its vertical profile in the optically thicker aerosol layers. In this case the extinction-to-backscatter ratio would be represented by a few cubic-spline basis functions. The reason for forcing backscatter-to-extinction to be smoothly varying is that it enables the high precision Mie channel to be used to infer the vertical structure, but with the Rayleigh channel ensuring the correct optical depth. If an arbitrary variation of backscatter-to-extinction is allowed, then the vertical structure comes entirely from the Rayleigh channel, which is much noisier. The a-priori for the backscatter-to-extinction ratio would ideally be based on a climatology of aerosol types in different locations around the globe (such a database is currently used by CALIPSO). It should be noted that so far the ACM-CAP algorithm has been tested only on HSRL data simulated from a scene retrieved from A-Train data in which no variation of backscatter-to-extinction ratio could be observed. Therefore, the retrieval of backscatter-to-extinction ratio has not yet been tested.
An additional state variable would be the aerosol median volumetric diameter; in the daytime it is possible that this can be estimated from the available shortwave radiances, provided that aerosol composition is assumed (e.g. from the ATLID aerosol classification). At night there are no observations that have any significant sensitivity to particle size, and therefore a constant value would need to be assumed (the retrieval would tend towards the a priori value naturally). To date, the aerosol component of the algorithm has only been tested on A-Train data at night, so a fixed aerosol median volumetric diameter of 0.5 microns has been assumed.
2 Microphysical and scattering assumptions
The following assumptions must be made in the aerosol retrieval:
• The shape of the size distribution; most commonly in aerosol work a lognormal distribution is assumed.
• The median volume diameter of the distribution (unless this information can be retrieved from shortwave radiances).
• The particle type, which informs the wavelength-dependent refractive index. Operationally this will be taken from an ATLID-only aerosol classification.
• The particle shape; most commonly spheres are assumed meaning that Mie scattering can be applied.
8 Radiative transfer forward models
In this section we outline the radiative transfer forward models that are used within each iteration of the algorithm to predict the observations from a profile of scattering properties that have been derived from the current values in the state vector, as described in section 5.5.5.4. Note that this just the last of the three yellow boxes in Figure 5. The scattering models used to convert physical properties of the particles into scattering properties were outlined in the previous section and summarised in Table 10. The radiative models used in this section are listed in Table 11.
|Table 11. Summary of radiative transfer forward models to be used for each observation type, together with their computational cost |
|(expressed in powers of the number of pixels in the profiles N). Multiscatter refers to the general-purpose active instrument |
|simulator code (but with particular capability for multiple scattering, MS), available at |
| (released under the GNU Lesser General Public License). A status of “OK” implies |
|that an adjoint code is also available. |
|Observation |Application |Speed |Status |
|Radar reflectivity |Multiscatter: single scattering option |N |OK |
|Radar reflectivity in deep |Multiscatter: single scattering plus TDTS MS model (Hogan and |N2 |OK |
|convection |Battaglia 2008) | | |
|Radar Doppler velocity |Single scattering easy; fast MS model with Doppler does not |N2 |Not available for MS |
| |exist | | |
|HSRL lidar in ice and aerosol |Multiscatter: PVC model (Hogan 2008) |N |OK |
|HSRL lidar in liquid cloud |Multiscatter: PVC plus TDTS models |N2 |OK |
|Lidar depolarization |Multiscatter: Developed and a publication is planned |N2 |Beta version only |
|Infrared radiances |Delanoe and Hogan (2008) two-stream source function method |N |Adjoint will be |
| | | |provided by Adept |
| | | |library |
|Infrared radiances |RTTOV (EUMETSAT license) |N |Tested (performance |
| | | |unsatisfactory) |
|Solar radiances |LIDORT (permissive license) |N |Testing |
1 CPR Reflectivity factor
Radar multiple scattering can be neglected in non-precipitating liquid and ice clouds, but in anything more intense than moderate rain from frontal rainclouds, and especially in deep convective clouds, multiple scattering becomes important (Battaglia et al. 2010). Therefore, a multiple scattering forward model is needed in any profiles containing rain. The time-dependent two-stream method of Hogan and Battaglia (2008) is a suitable method for calculating radar multiple scattering very rapidly. It has been demonstrated to agree well in comparison to Monte Carlo calculations despite being around seven orders of magnitude faster. An exact adjoint has been written for this code, enabling it to be used in quasi-Newton type solvers; this approach has been demonstrated with this forward model by Pounder et al. (2012). The code runs in milliseconds, but the execution time scales as the number of points in the column squared. This is fast enough to use in a retrieval algorithm, but to provide even more speed, an adaptive component has been added: this triggers the more expensive multiple-scattering part of the calculation only if the transport mean free path is less than 20 times the receiver footprint at the altitude of the cloud. The surface return is not yet included in this model.
2 CPR Doppler velocity
Forward modelling of the reflectivity-weighted terminal fall speeds is straightforward since the terminal fall speed for individual particles is quite well known given their mass, maximum dimension and projected area, and is done as part of the first two yellow boxes in Figure 5. In the absence of multiple scattering, no forward model needs to be run; this profile of Doppler velocity is the best estimate of the measured Doppler velocity. In the presence of multiple scattering, the results of Battaglia et al. (2011) suggest that the sideways spacecraft motions from lateral scattering somehow cancel out, with the result that the measured Doppler velocity is the mean of the vertical velocity from the places in the cloud where the radiation has been scattered from. Note that often CloudSat profiles in deep convection have the appearance of extending down to the surface, when in fact the return is dominated by multiple scattering near cloud top. In this case the apparent Doppler velocity from low in the profile would actually be close to the mean value over a range of heights much higher in the cloud. This effect could be added to the Hogan and Battaglia (2008) code relatively straightforwardly, but this has yet to be done.
It is anticipated that the greatest difficulties in practice will be in characterizing the errors in the measurements (which need to be included in a variational retrieval), something that is difficult to achieve completely before launch. The sources of error are:
• Antenna mispointing. If the radar is not pointing completely at nadir then a component of the spacecraft motion will be measured. In principle this can be removed by subtracting the Doppler of the ground return.
• Non-uniform beam filling (NUBF). When there is a large horizontal gradient in reflectivity from one side of the volume to the other (such as in the vicinity of convective clouds), the radar will receive a Doppler return preferentially from one side of the volume and therefore measure a component of the spacecraft motion. The C-CD product will attempt to use the resolved along-track gradient in reflectivity to correct for this effect, but such a procedure would not be perfect. We will therefore require the Doppler velocity error to be provided with this product in order that the best estimate algorithm can weight the Doppler information less in regions of NUBF.
• Misrepresentation of multiple scattering. When Doppler is added to the multiple scattering forward model, it will be in an approximate form. It remains to be seen whether this modelling is accurate enough to use directly, or if we should use the multiple scattering flag in the C-CD product to disregard Doppler measurements likely to be subject to multiple scattering.
• Mistreatment of vertical wind. In stratiform clouds, particularly cirrus, the vertical wind is much smaller than the terminal fall velocity (particularly when averaged over the 1-km radar footprint) and it is reasonable to neglect it in the retrieval except in estimating its standard deviation to use as the error in the forward modelled Doppler velocity. Unfortunately, vertical wind is correlated in height, which means that to treat it correctly, the observation error covariance matrix R ought to have off-diagonal components.
The last three of the four bullets above are all particularly problematic in the vicinity of deep convection: here horizontal gradients of reflectivity are largest, multiple scattering is the most significant and the vertical winds are strongest and most turbulent. Therefore it may not be possible to use the Doppler signal for microphysical retrievals in these circumstances, in which case the target classification information would be used to identify these regions and ensure that the Doppler signal in them is not added to the observation vector. This implies that it may not be essential to include the Doppler return in the multiple scattering model. In testing so far we have simply calculated what CloudSat would see if it had Doppler capability but in the absence of either multiple scattering or vertical wind.
3 CPR Surface echo
When viewed from nadir, the surface echo return can be thought of as specular backscatter modified by the surface ocean waves induced by wind stress (Kudryavtsev 2003). Accurate modelling of the ocean surface cross section is required for the use of this information in both liquid cloud and precipitation retrieval techniques that use the path integrated attenuation (PIA) method. The stability of the ocean as a target at 94-GHz was investigated using airborne radar by Li et al. (2005). The PIA is simple to forward model from the optical depth δ without using a full radiative transfer model, via PIA[dB] = 10log10(e2δ) = 8.686δ. Estimating PIA observationally requires modelling not only its dependence on wind and temperature, but also interpolating spatially between cloud-free regions. Since PIA is of use to other algorithms, it is not intended that it forms part of the Best Estimate algorithm, rather that PIA is a variable in the merged observations or target classification products.
4 ATLID HSRL Backscatter
For the lidar forward modelling, it is important to take account of multiple scattering. We will use the fast but accurate model “Photon Variance-Covariance” (PVC) model of Hogan (2006, 2008) for ice clouds and aerosols, in which only small-angle multiple scatter is significant. An exact adjoint model has been coded and debugged. Moreover, this model naturally predicts the cloud/aerosol and molecular returns separately, and therefore can be used to forward model the separate Mie and Rayleigh channels of an HSRL lidar. This was used in [CASPER-ATBD] to adapt the Delanoë and Hogan (2008) algorithm to EarthCARE.
In optically thick liquid clouds there can be significant returns from wide-angle multiple scattering, which results in apparent “pulse stretching”. Hogan and Battaglia (2008) provided a fast algorithm based on the Time-Dependent Two-Stream (TDTS) equations, which can be used in conjunction with the PVC algorithm. Indeed, the multiscatter software package (see Table 11) merges these two components automatically. Pounder et al. (2012) demonstrated that this model can be used in a retrieval scheme. The model assumes that multiple scattering by cloud and aerosol particles does not widen the spectral return, and therefore that these signals are still measured only by the Mie channel of the HSRL. It is possible that the lateral transport in multiple scattering leads to a random Doppler frequency shift that has the effect of widening the Doppler spectrum and thereby causing a leakage into the Rayleigh channel, but the exact characteristics of the HSRL filters (as well as the spacecraft motion) would be required to model this effect.
5 ATLID Depolarization
Depolarization in ice clouds is very difficult to interpret quantitatively due to the unpredictable dependence on particle shape. In liquid clouds subject to single scattering the depolarization is essentially zero, but depolarization can occur when multiply scattered photons are detected. Recently we used the framework provided by the TDTS methodology to estimate the distribution of scattering orders at each apparent range gate, and then parametrically predict the typical depolarization associated with each order of scattering. The idea is that the TDTS model predicts the average number of wide-angle scattering events that have been experienced by radiation returning from a particular gate at a particular time. The depolarization of that return may then be parameterized as a function of the number of scattering events, since each one will increase the degree of depolarization. The PVC part of the multiple scattering package was also modified to simulate depolarization. The result has been shown to work for both idealized case of isotropic scattering, and for Mie phase functions in clouds with a simple geometry (a “top-hat” extinction profile and a constant effective radius). An earlier stage of the work in this area is described in [RATEC-FR], but a little further development would be needed before the code could be used in a retrieval algorithm. In any case, the information on extinction coefficient near cloud top duplicates the information available more directly from the HSRL capability of EarthCARE so it is of higher priority to work on other forward modelling aspects than to add a depolarization forward modelling capability.
6 MSI Thermal infrared radiances
We use the Delanoë and Hogan (2008) infrared radiance model, which includes scattering by use of the two-stream source function technique (Toon et al. 1989), as well as an accurate treatment of gas absorption and emission via correlated-k-distribution look-up tables (taken from ECSIM). An approximate Jacobian matrix may be computed rapidly using the absorption approximation, and the Adept library will be able to provide the adjoint. The method performs well in comparison to both DISORT and the ECSIM Monte Carlo model in both cloud and clear sky, as is shown in [RATEC-FR].
The standard package for forward modelling infrared and microwave radiances is “RTTOV” which is widely used in data assimilation, particularly at ECMWF and the UK Met Office. The most recent version, RTTOV-9 (Saunders et al. 2007) has detailed treatment of clouds and includes both adjoint and Jacobian calculations, so would appear to be an alternative to use in the framework of a combined active-passive retrieval. This package is widely used and the work was funded by EUMETSAT, so no licensing problems are anticipated with its use in an ESA retrieval algorithm. However, as shown in [RATEC-FR], agreement with the Delanoe and Hogan, DISORT and ECSIM Monte Carlo codes is poor, so we do not have plans currently to make it the primary infrared radiance model for the unified retrieval.
7 MSI Solar radiances
It was demonstrated by Polonsky et al. (2008) how solar radiances could be modelled within the optimal estimation framework using the “RADIANT” radiance forward model, and this is now used as part of the official CloudSat-MODIS retrieval of cloud properties. Little further development is being performed with RADIANT, but the LIDORT model provides a similar capability (for example, a fast analytic Jacobian model, which is equivalent to an adjoint model for passive instruments with only one output radiance) and has been recommended as a replacement. These models provide a multi-stream simulation of the radiation field, but need fewer vertical levels than an infrared model due to there being no need to resolve the variation of the Planck function with temperature and therefore height. We have compared LIDORT to the widely-used SBDART model (the former is designed for retrievals while the latter is designed for off-line research calculations). [RATEC-FR] presented preliminary results that appeared to show large errors for LIDORT, but an issue in the comparison has been resolved and the performance is much better. Therefore we are believe LIDORT will be a viable model for operational use.
There are three main source of inaccuracy in modelling solar radiances:
• Uncertainty in the surface albedo. Except for the most optically thick clouds, the surface albedo has a significant role in determining a solar radiance, and this is uncertain over many land surfaces.
• Uncertainty in the detailed shape of the scattering phase function, which is most important for optically thin clouds. The Baran (2004) database of ice particle phase functions extends from the UV to the far infrared and so a consistent model could be used across this range of the spectrum.
• Three-dimensional scattering effects. These would be virtually impossible to incorporate into a forward model, which for efficiency must make the plane-parallel approximation. It is conceivable that the texture of the MSI image could be used to estimate the likely error due to this effect.
Algorithm performance, sensitivity studies, limitations
1 Computational cost
Since introducing automatic differentiation to replace the use of numerical gradients (perturbing each state variable in turn), the algorithm has sped up considerably. Figure 8 shows a breakdown of the mean time spent per profile on different activities, for a dataset of test EarthCARE data generated from an A-Train scene. In this case the mean total time per profile was 0.3 s. We see that negligible time is spent in the minimiser computing the next search direction, or reading or writing data. The error computation, which involves the most intense matrix operations, is a relatively small fraction of the total cost. The two largest bars are:
o Active instrument radiative transfer. It is not surprising that this occupies a significant fraction of the overall time, since the computational time of the Hogan and Battaglia (2008) model for wide-angle multiple scattering goes as the number of range-gates squared. We have already optimized this code by only applying it when it is needed. A possible further opportunity for optimization lies in the possibility to run it at half the resolution of the raw observations; since multiple scattering tends to smooth out sharp features in the cloud field, this is not expected to introduce significant error, but could gain a factor of four in speed.
o Automatic adjoint. Hogan (2013) showed that the Adept library is significantly faster than other comparable libraries and computes adjoints not much slower than hand-coded versions, so it is difficult to see how much further reduction could be obtained in this bar.
An important point to note is the fifth to the eleventh bars scale with the number of iterations required, so there is much to be gained by using a more efficient minimizer that can find the minimum of the cost function in fewer iterations. Typically at least a factor of 5 reduction in the number of iterations can be achieved by using the Levenberg-Marquardt method, but at the cost of having to compute the Jacobian matrix each iteration. If this Jacobian computation can be parallelized then Levenberg-Marquardt algorithm suddenly becomes much more efficient than the L-BFGS quasi-Newton method shown here. We have recently found very good parallelization of Jacobian computations using an OpenMP on multi-core machines, so there is reason to believe that significant speed-up of the ACM-CAP algorithm than shown in Figure 8 is possible.
|[pic] |
|Figure 8. Summary of the mean time per profile spent in each part of the algorithm, averaged over the processing of 1400 rays of simulated |
|EarthCARE data (“test2” in the Algorithm Test Plan). In this scene the mean total time spent on each profile was 0.3 s. |
2 Limitations
There are numerous limitations of the Best Estimate retrieval, but virtually all are limitations that are shared by other retrieval schemes, or are related to optional capabilities of the Best Estimate algorithm that other algorithms do not even have. These include the following:
• In ice clouds, radar and lidar retrievals of Delanoe and Hogan (2008), which this algorithm will closely replicate, produce some of the most rigorously derived remote retrievals of ice cloud properties. It had previously been found that if infrared radiances were incorporated, they always tended to increase the extinction towards the cloud top, implying that there was a bias in the assumed ice particle scattering properties at cold temperatures. However, a bug has since been found in the implementation of the infrared model in the Delanoe and Hogan (2008) algorithm, not replicated in ACM-CAP, so further work is required to explore how well the radar-lidar retrieved clouds are consistent with the information from infrared radiances.
• The use of longwave radiances requires a good estimate of atmospheric temperature in order for them to improve the radar-lidar retrievals, and if the cloud is optically thin then an accurate surface temperature and emissivity is also required. If these are inaccurate then there is a danger that the extra observational information can actually degrade the retrieval. There is a similar need for accurate surface albedo to be provided when shortwave radiances are to be assimilated. There are regions of the globe where these variables are likely to be less accurate, particularly areas with seasonal snow or ice cover. At the very least, realistic errors in these quantities will need to be provided in order that they can be incorporated as forward model errors (see Delanoe and Hogan 2010).
• As discussed in section 5.5.7.2, liquid clouds are rather troublesome to retrieve from space due to the frequent lack of information on the location of cloud base or the presence and properties of multiple liquid layers beneath which liquid or ice precipitation is falling. This limitation is common to all retrieval algorithms, but the methods proposed in section 5.5.7.2 are hopefully able to extract what information there exists in the observations, making use of physical constraints, better than any other in current use.
• The retrieval of microphysical information in deep convection is very difficult because the lidar will lose signal rapidly, the Doppler signal will be dominated by convective motions, the radar is affected by multiple scattering and there are a wide range of particle types present in the mixed-phase regions of such clouds.
Validation status
1 Maturity
Development of the Best Estimate algorithm was started in an earlier ESA project (reported in [RATEC-ATBD]) with a working code produced. There is still considerable development needed both in terms of coding and refining the assumptions and constraints used for each retrieved constituent. The ice-cloud part of the algorithm is the most mature, being based on the Delanoe and Hogan (2008) scheme, which has been applied successfully to ground-based, airborne and over 3 years of A-Train data.
To illustrate the workings of the Best Estimate algorithm, Figure 9 shows the results of running the retrieval algorithm on a north-Atlantic cloud system observed by CloudSat and Calipso. Comparing the forward modelled radar and lidar signals to those measured, we can see that the algorithm has successfully converged, matching the observations in each constituent type. The fact that the lidar forward model matches the molecular return below the ice cloud (but without reproducing the speckle noise) indicates that it has successfully used this information to provide an optical depth constraint on the ice cloud retrieval. The bottom-right image shows the forward modelled Doppler velocity in the absence of either vertical wind or multiple scattering. Encouragingly, this field has a very similar appearance to ground-based Doppler radar observations, with the ice falling at a rate that increases down towards the melting layer, followed by a rapid increase in the rain.
Some of the three main retrieved consistuents are depicted in Figure 10 for ice, rain and liquid cloud. Naturally we have no truth to compare against, but some comments can be made on the realism of the retrievals. In the case of liquid clouds, the information in this case is essentially coming only from the lidar, since the radar return is dominated by drizzle drops (which are also retrieved). Since the lidar backscatter is rapidly extinguished, it is very difficult to infer the structure of optically thick clouds. It was shown by Pounder et al. (2012) that multiple-scattering ought to enable optical depth up to around 35 to be retrievable from Calipso within a variational scheme that models the multiple scattering, like this one. However, in this case it appears that there is a problem with the target classification that makes this difficult to achieve; specifically, the radar sees cloud top higher than the lidar due to the long radar pulse length. The target classification has diagnosed liquid cloud in the region above the cloud top where the radar receives an erroneous signal. This needs to be fixed before we can really investigate in detail the liquid cloud retrieval. Note the final EarthCARE algorithm we will also incorporate the HSRL information for the extinction coefficient at cloud top, and the solar radiances and path-integrated attenuation for constraints on the optical depth and liquid water path of the cloud, respectively. The performance of the gradient constraint for liquid cloud retrievals is explored in the next section in a simulated example.
|[pic] |
|[pic] |
|Figure 9. Demonstration of the retrieval for a deep cloud observed by CloudSat and Calipso. (Top panel) Forward-modelled attenuated lidar |
|backscatter coefficient at the final iteration of the algorithm, below which are the corresponding observed values. (Third panel) Forward |
|modelled radar reflectivity factor, below which are the corresponding observed values. (Bottom panel) The forward modelled Doppler |
|velocity (positive downwards) that would be measured in the absence of vertical wind and multiple scattering. The retrieved variables are |
|shown in Figure 10. |
The ice cloud retrieval reveals a generally smooth field but with some horizontal discontinuities due to the fact that each ray is retrieved independently. This retrieval is quite similar to that of Delanoe and Hogan (2010), although the actual retrieved values have not been compared for the same case.
The rain retrieval has values of several mm h-1 in the deep cloud, and less than 0.1 mm h-1 in the drizzling stratocumulus at the beginning of the period. A flatness constraint has ensured a slowly changing profile in the vertical, as well as to extend the rain retrieval (essentially to extrapolate) down to the ground where the radar suffers from ground clutter. The attenuation by the melting layer has been included in a simple way: by assuming that the two-way attenuation is uniquely and predictably related to the rain rate.
These results are encouraging and demonstrate that the use of a Best Estimate algorithm can retrieve realistic values. Naturally, validation is needed, although it is likely that any differences can be remedied simply by changing assumptions made in the retrieval, rather than being due to fundamental problems in the way the variational retrieval has been formulated.
A further example is shown in Figure 11, this time with a mid-level mixed-phase cloud and an ice cloud above. The high lidar backscatter at the top of the mid-level cloud is interpreted as a layer of liquid water, with ice falling from it. An interesting aspect of this case is that there is probably liquid cloud atop the entirety of this cloud, but the optical depth of the cloud above varies so the lidar backscatter from the liquid is often strongly attenuated. The algorithm has difficulty in accounting for this and tends to retrieve lower liquid water content in clouds more obscured by ice clouds above. It is hoped that with visible radiances providing a better constraint on optical depth, some of this could be remedied.
|[pic] |
|[pic] |
|[pic] |
|Figure 10. Main retrieved atmospheric constituents for the case shown in Figure 9. |
|[pic] | |
|[pic] |[pic] |
|Figure 11. (Left first two panels) Calipso lidar backscatter and CloudSat radar reflectivity factor for a mixed-phase cloud with an ice cloud above. |
|(Bottom-left panel) Forward modelled Doppler velocity. (Right panels) Retrieved variables. |
2 Liquid cloud example
Figure 12 illustrates the retrieval of liquid cloud using the Best Estimate algorithm, but applying it to simulated profiles using only an HSRL lidar with the EarthCARE specification, including instrument noise.
|[pic] |[pic] |[pic] |
|Figure 12. Demonstration of retrieval of liquid water content profile from ATLID using simulated profiles and the Best Estimate algorithm. The top row |
|shows the Mie and Rayleigh channels while the bottom row shows the true and retrieved profiles in each case. The left profile is of an adiabatic cloud |
|retrieved using a smoothness constraint, the middle profile considers the same profile retrieved using a one-sided gradient constraint, while the right |
|panel considers a cloud that is adiabatic only in the lower half again retrieved using the one-sided gradient constraint. |
The left column illustrates the retrieval of an adiabatic profile using a smoothness constraint on the logarithm of liquid water content. The result is a retrieved profile that has approximately the correct optical depth and liquid water path, but the vertical distribution of liquid water is not correct (since the constraint makes it tend towards a straight line in logarithmic space). The centre collumn illustrates the retrieval of the same profile but using a one-sided gradient constraint described above. Here the retrieval is almost perfect. The right column illustrates the retrieval of a cloud that is adiabatic at the base but “flat” in the top half. The retrieved shape is not quite so good. It should be noted that when the wider Calipso field-of-view is used (not shown), the retrieval is very much better, because the extra multiply scattered radiation from deeper in the cloud provides information on the extinction profile there.
Annex A: Implementation details and supporting software
The annex describes some of the implementation details in the particular implementation that has been coded so far.
1 Automatic differentiation
1 The problem
Coding up both the exact adjoint needed by the quasi-Newton method (section 3.4.2), and the exact Jacobian needed for the error covariance matrix of the solution, is laborious and error-prone. This can significantly slow down algorithm development since every change to the algorithm that changes the computation of the cost function (which is almost everything) also requires a working adjoint in order to be tested and a working Jacobian to be used operationally. Before presenting a solution to this problem, we first illustrate how an adjoint code would formally be written for this algorithm.
Often the minimizer requests not just the cost function J for a given trial state vector x, but also its gradient (xJ. The process of computing the gradient is to split the computation into the contribution to the gradient from each observation and each constituent. As shown in (1), the contribution to the cost function from constituents (specifically from any constraints associated with these constituents) comes from the last two terms, which can be differentiated to get an expression that can be computed directly:
[pic]
The contribution to the gradient from the observations is then added. It is inefficient to use the first term on the right-hand-side of (3) explicitly (also shown by Eq. 7). Instead we would code up the adjoint of all the steps in the forward model. This can be thought of as a standard exercise in mathematics and programming in which all the mathematical statements in the program are converted into adjoint form, then executed in reverse order. Thus for every function in the program we must write its equivalent adjoint function. Suppose a function is of the form b = f(a), where the inputs and outputs are represented as vectors, then the adjoint function is of the form dJ/da = f*(a, dJ/db). A relatively accessible description of how to do adjoint coding was provided by Giering and Kaminski (1989), which we will not repeat here. Nonetheless, it is useful to see how the order of operations in section 5.5.4 are reversed. We may summarize the calculation of the contribution to the cost function from a particular observation by the following:
1. Calculate scattering properties for that observation:
a. Loop through each constituent and use look-up table to calculate scattering properties of that constituent at the wavelength of the observation, e.g. to yield extinction coefficient αi for constituent i
b. Merge properties from different constituents, e.g. to get total extinction coefficient α
2. Run radiative transfer forward model for that observation to get yf
3. For each variable in that observation, calculate cost function contribution J
The equivalent adjoint code would do the same, but then follow it by:
3* For each variable, calculate gradient dJ/dyf = R-1[yf – yo]
2* Run adjoint of radiative transfer code for that observation, e.g. to get dJ/dα
1* Calculate gradient of scattering properties for that observation
b* Extract gradient with respect to the scattering properties of each constituent from total, e.g. to get dJ/dαi for constituent i
a* For each constituent: update gradient dJ/dx
Since the order of every operation is reversed and lots of intermediate variables need to be stored in the forward pass in order to be used in the reverse pass, adjoint coding is error prone and difficult to debug. For this reason, during early development and testing of the algorithm, we implemented a function that calculates the gradient numerically simply by individually perturbing each element of the state vector by 0.001 and comparing the resulting cost function to the unperturbed value. Naturally this was many times slower than a full adjoint code, especially for large state vectors, but allowed the forward model to be developed and tested easily. Naturally this is too slow for an operational algorithm, so was not a long-term solution: some way of computing the adjoint exactly and cheaply (both in computer and programmer time) is required.
2 A solution
An solution is to generate the adjoint or Jacobian automatically, for which two possibilities exist. The first is a source-to-source compiler, which analyses the original code and generates the source for the equivalent adjoint code, which can then be compiled into the executable. The problem with this approach is that existing programs for doing it are very limited in the range of language features that they recognise, preventing user-defined classes, operator overloading and templates. The second possibility is the operator overloading technique, where variables in the code for which a gradient is required are declared as a special “active” type, and then their use in mathematical expressions leads to then necessary information being automatically stored such that the adjoint can be computed automatically. This is compatible with all language features, but a concern is that existing libraries for automatic differentiation by operator overloading, such as CppAD and ADOL-C, are typically 5-10 times slower than hand-coded adjoints.
However, a breakthrough was recently made by Hogan (2013) to show that if an operator overloading library is coded using “expression templates” in C++, then the adjoint can be computed only 5-10% slower than hand-coded adjoints, i.e. 5-9 times faster than CppAD and ADOL-C. The resulting library, Adept (Automatic Differentiation using Expression Templates) is now being used to compute adjoints and Jacobian matrices automatically, resulting in a very significant speed-up compared to the use of numerically computed gradients.
2 Scattering and terminal fall-speed calculations
The scattering look-up tables generated internally by the algorithm (described in section 5.5.2.2) provide scattering properties for full distributions of particles in a volume. The function that calculates these look-up tables takes as input the precalculated scattering properties of individual particles stored in a scattering library file. This section is concerned with how these particle scattering library files are generated.
In the case of ice particles observed by instruments with a wavelength shorter than around 15 microns (e.g. lidar and imager), we use the Baran (2004) database of ice-particle scattering properties. For all other particles, and for ice particles observed at radar wavelengths, scattering libraries for individual particles may be created using the scatter program that is part of the unified algorithm distribution. This program currently uses either Mie theory (using the freely available BHMIE Fortran back end), T-matrix for oblate spheroids or Self-Similar Rayleigh Gans for low-density aggregates (Hogan and Westbrook 2013). Scattering properties may be calculated as a function of wavelength, temperature, size and mixing fraction (for aggregated ice particles that may be treated as a mixture of ice and air). Thus the output variables may be up to four dimensional, although frequently not all these dimensions are required. In addition to electromagnetic scattering, the scatter program also calculates particle terminal fall speeds using the Heymsfield and Westbrook (2010) approach for ice particles and the Beard (1976) approach for liquid droplets and raindrops. The ASCII input format is as in section 5.5.2.1; to create the particle scattering properties for ice at 94 GHz, the input file ice_94-GHz.cfg contains the following:
# Wavelength corresponding to 94-GHz, in m
wavelength 0.0031893
# Ice spheres
medium ice
# Use Mie theory
engine bhmie
# Temperatures from -80 degC to 0 degC in 10 degC intervals
temperature { 193.15 203.15 213.15 223.15 233.15
243.15 253.15 263.15 273.15 }
# Define size in terms of radius (m), with logarithmically-spaced
# sizes from 0.5 microns to 1 cm
size_variable radius
size_start 0.5e-6
size_end 0.01
size_values 100
size_scale logarithmic
# Include ice-air mixtures with a full range of ice fractions, using
# the Maxwell-Garnet mixing rule
fraction_start 0.02
fraction_end 1.0
fraction_values 50
fraction_scale linear
# Output format
format netcdf
# Output filename
output ice_94-GHz_scat.nc
Note that the refractive index is computed from the medium, temperature and fraction. The output NetCDF file, ice_94-GHz_scat.nc, then contains the following coordinate variables [with units]:
• wavelength [m]
• temperature [K]
• fraction
• mean_dimension [m]
It contains the following particle physical properties (functions of the dimensions specified in brackets):
• area (size) [m2]: the geometric projected area of the particle when viewed from zenith or nadir – this is required, for example, by the multiple scattering forward model.
• mass (fraction, size) [kg]
• fall_speed (fraction, size) [m s-1]: the terminal fall speed at a reference air density of 1.0 kg m-3 – this is required by the Doppler forward model and to calculate the snow rate.
• n_r (wavelength, temperature, fraction): the real part of the refractive index.
• n_i (wavelength, temperature, fraction): the imaginary part of the refractive index.
It also contains the following particle scattering properties:
• ext (wavelength, temperature, fraction, size) [m2]: extinction cross-section.
• scat (wavelength, temperature, fraction, size) [m2]: scattering cross-section.
• bscat (wavelength, temperature, fraction, size) [m2]: backscattering cross-section.
• g (wavelength, temperature, fraction, size): asymmetry factor (the mean of the cosine of the scattering angle.
At solar wavelengths for radiance forward models, it may be necessary to include more details of the scattering phase function; these will be added when solar radiances are incorporated into the algorithm.
3 Efficient storage of inverse of background error covariance matrix
Commonly the background (or prior) error covariance matrix, B, has a simple form, resulting in its inverse being sparse. If this can be exploited then matrix-vector multiplications involving B-1 can be much more efficient. In the trivial case where B is diagonal and the error variances on the diagonal are equal at σ2, then B-1 is also diagonal with the diagonal elements all equal to σ -2. Commonly we use off-diagonal elements in B to represent error correlations, which has the effect of smoothing the retrieval. If the correlation coefficient coefficient between errors in adjacent elements of the state vector is ρ, then we would write the error covariance matrix as
[pic].
Equivalently we may write each element in terms of a decorrelation length z0:
[pic]
The inverse of this is tridiagonal, and therefore efficient to store and use:
[pic],
Where a = 1/(1 – ρ2), b = (1 + ρ2)/(1 – ρ2) and c = – ρ/(1 – ρ2). This is trivially extended to larger matrices.
-----------------------
[1] 湧牯⽧潳瑦慷敲术汳ᔯഠ 瑨灴⼺眯睷挮潨歫湡漮杲猯景睴牡⽥楬汢晢獧യ –奈䕐䱒义⁋栢瑴㩰⼯睷潲煣椮牮慩昮⽲松汩敢瑲洯摯汵灯⽴灯楴業慺楴湯爭畯楴敮⽳ㅭ湱⼳ㅭ湱⸳瑨汭•ᐁ瑨灴⼺眯睷爭捯湩楲牦縯楧扬牥⽴潭畤潬瑰漯瑰浩穩瑡潩潲瑵湩獥洯焱㍮洯焱㍮栮浴ᕬഠ̍ഄ̍ഄ䄍䵃䌭偁䄠䉔ൄ䅖卒⁙牐橯捥ݴ潃敤䰇戲䄭䵃䌭偁䄭䉔݄܇獉畳ݥ㌰܇䐇瑡ݥ〳〯⼹〲㌱܇倇条ݥ–䅐䕇ᐠ〱―景ጠ丠䵕䅐䕇⁓尠䕍䝒䙅剏䅍w.software/gsl/
[2]
[3]
-----------------------
Not converged
Converged
3. Line search
Repeated calls to forward model to find minimum J along this direction
First guess of state vector
2. Chose search direction
Either L-BFGS or conjugate gradient method (first iteration will be in direction of steepest descent)
4. Test for convergence
1. Forward model
Use state vector x to calculate cost function J and, if requested, gradient (xJ
1. Forward model
Use state vector x to calculate cost function J and its gradient (xJ
Converged
Not
converged
1. Load scattering libraries
2. Load new ray of data
Use classification to define state vector and its first guess
3. Minimize cost function
Iterative method to calculate best estimate of the state vector
4. Calculate retrieval error
Store retrieved variables from ray
Not converged
Converged
Aerosol
Constituent
Observation
State
Contains lists of Observations and Constituents
ActiveInstrument
Radar
Lidar
Radiometer
BackgroundConstituents
LiquidCloud
IceCloud
ConcreteClass
AbstractBaseClass
InfraredRadiometer
ShortwaveRadiometer
x
H(x)
................
................
In order to avoid copyright disputes, this page is only a partial summary.
To fulfill the demand for quickly locating and searching documents.
It is intelligent file search solution for home and business.
Related searches
- analytics and reporting best practices
- finance and accounting best practices
- poets and quants best online mba
- planning and budgeting best practices
- diversity and inclusion best practices
- 90 day precipitation forecast
- best cloud management software
- best free cloud app
- planning and scheduling best practices
- melatonin what it does and the best ways to take it
- noaa 90 day precipitation forecast
- best cloud based phone systems