VARSY Project / Code / L2b-ACM-CAP-FR
Issue / 01
Date / 12/09/2013
Page / 1 of 24
EarthCARE Level 2 Documentation
ACM-CAP Cloud, Aerosol and Precipitation “Best Estimate”
Final Report
VARSY Project
Code: / L2b-ACM-CAP-FR
Issue: / 01
Date: / 12/09/2013
Reference: / University of Reading
Name / Function / Signature
Prepared by / Robin Hogan / Project Scientists
Reviewed by / Pavlos Kollias / Project Scientist
Approved by / Pavlos Kollias / Project Manager
Signatures and approvals on original
This page intentionally left blank
Document Information
Contract DataContract Number: / 4000104528/11/NL/CT
Contract Issuer: / ESA-ESTEC
Internal Distribution
Name / Unit / Copies
Robin Hogan / University of Reading / 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
Pavlos Kollias / McGill University / 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-READING-FR-001
Document Status Log
Issue / Change description / Date / Approved01 / First version delivered to ESTEC / 12/09/2013
Table of Contents
1. Purpose and Scope______
2. Documents______
2.1. Applicable Documents______
2.2. Reference Documents______
2.3. List of Abbreviations______
3. INTRODUCTION______
4. Achievements within the project______
5. Current Status______
List of Tables
Table 1: Applicable Documents......
Table 2: Reference Documents......
Table 3: List of abbreviations......
1.Purpose and Scope
This document briefly summaries the work performed at Reading within VARSY project. The specific work that has been performed revolved around the development of the software for the ACM-CAP product.
2.Documents
2.1.Applicable Documents
The following table specifies the applicable documents that shall be complied during the project development.
Table 1: Applicable Documents
Reference / Code / Title / Issue[SOW] / EC-SW-ESA-SY-0310 / Statement of Work: VARSY - 1-Dimensional VARiational Retrieval of SYnergistic EarthCARE Products / 1.0
[CC] / Appendix 2 to
AO/1-6823/11/NL/CT / Draft Contract (attachment to SOW) / 1.0
[AD 1] / EC-SW-ESA-SY-0152 / EarthCARE Level 2 Processor Development
General Requirements Baseline / 1.0
[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:
Common Product Definitions / 1.0
[AD 4] / EC.ICD.ASD.ATL.00021 / EarthCARE Product Definitions. Vol. 2b:
ATLID level 1 / 1.0
[AD 5] / EC.ICD.ASD.BBR.00022 / EarthCARE Product Definitions. Vol. 3b:
BBR level 1 / 1.0
[AD 6] / EC.ICD.ASD.MSI.00023 / EarthCARE Product Definitions. Vol. 4b:
MSI level 1 / 1.0
[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 EarthCARE Ground Segment / 2.0
2.2.Reference Documents
The following table specifies the reference documents that shall be taken into account during the project development.
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
[ATLAS-FR] / EC-FR-KNMI-ATL-027 / ATLAS Final report / 1.0
[ATLAS-ACM-TC] / EC-TN-KNMI-ATL-ACM-TC-024 / L2b Classification ATBD / 1.2
13/03/08
[ATLAS-EBD] / EC-TN-KNMI-ATL-ATBD-A-EBD-021 / L2a ATLID Extinction, Backscatter and Depolarization algorithm ATBD / 1.1 27/04/09
[ATLAS-FM] / EC-TN-KNMI-ATL-ATBD-A-FM-010 / L2a ATLID Feature mask ATBD / 2.2
[RATEC-FR] / RATEC-FR-READING-1 / RATEC Final Report / 1.0, April 2011
2.3.List of Abbreviations
Table 3: List of abbreviations
Abbreviation / Name1D-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)
3.INTRODUCTION
This document summarizes the work carried out at the University of Reading during the VARSY activity in the development of the “Cloud, Aerosol and Precipitation Best Estimate” product. It should be read in combination with the Product and Algorithm Requirements Document [PARD], which contains an extensive review of the previous relevant synergy retrieval work, the Product Definition Document [PDD], describing the product, and the Algorithm Theoretical Basis Document [ATBD], which describes in detail how the algorithm works. Given the high level of detail in these other (longer) documents, section 4 contains a briefer overview of the work and several short accounts of contributing studies that assisted in the ongoing development of the Best Estimate algorithm. Section 5 then summarises our progress in the development of the algorithm since the end of the previous RATEC project, along with the remaining work to be done.
4.Achievements within the project
4.1.One-sided gradient constraint for liquid clouds
An important physical constraint for liquid clouds is that the gradient of liquid water content with height should not exceed the adiabatic value, which is a function of temperature and pressure. In the [ATBD] it is explained how a “one-sided gradient constraint” has been added, which adds a term to the cost function if this gradient is exceeded. Here we evaluate the liquid-cloud part of the retrieval by comparing it to independent data.
Figure 1. (from top) Calipso lidar backscatter observed for a stratocumulus cloud observed over ocean; the forward-modelled values at the final iteration of a retrieval using ACM-CAP making use only of the lidar; the corresponding retrieved liquid water content; the corresponding atmospheric optical depth at the lidar wavelength.Figure 1 depicts a retrieval of the liquid water content using purely the Calipso lidar return. A one-sided gradient constraint has been used. It can be seen that by forward modelling the multiple scattering signal of the lidar, optical depths to 60 are retrieved; if it were not for multiple scattering, retrievals above an optical depth of around 5 would not be possible.
Independent verification of these retrievals is provided in Figure 2a, which shows a very good agreement between observed and simulated radar path-integrated attenuation; this variable is proportional to liquid water path. The one location where the agreement is not so good is in the region of highest PIA; presumably this is where the optical depth is so high that even with multiple scattering, there is not enough information in the multiply scattered returned light to infer the full optical depth of the cloud. Figure 2b demonstrates that this may be overcome by assimilating also the PIA in the retrieval.
Figure 2. Comparison of the observed and forward modelled radar-path integrated attenuation (which is proportional to liquid water path) for the retrieval shown in Figure 1. The top panel shows a comparison for the retrieval using only the lidar, while the bottom panel shows a retrieval where the path-integrated attenuation was assimilated too.These results are for the Calipso lidar, which has a field of view of around 90 m. The EarthCARE lidar has a much narrower field of view, which reduces the strength of the multiple scattering effect that is providing most of the information on optical depths above around 5. We have performed simulations and find that there is a critical field of view of around 50 m, below which it is no longer possible to use multiple scattering to infer the optical depth of thick clouds. Unfortunately, EarthCARE has a smaller value than this. Therefore, the retrieval will need to rely more on the other information available, such as the solar radiances (in the day only), the path integrated attenuation (useful over oceans only), the HSRL capability of the lidar and of course the gradient constraint.
4.2.Automatic differentiation library “Adept”
An outstanding question at the end of the RATEC project and the start of the VARSY project was the problem of coding up by hand all the adjoints of the various components of the retrieval algorithm, given how time consuming and error prone this activity is. As solution has been found in that we have developed a new automatic-differentiation software library “Adept”, which is easy to incorporate into an existing C++ project and can compute both adjoints and Jacobian matrices only slightly slower than it would take a hand-coded function (although without the long and difficult process of doing the hand-coding). Note that other automatic-differentiation libraries are available that build on the same operator-overloading idea (e.g. ADOL-C, CppAD and Sacado), but they are all considerably slower than Adept. This work is in the process of being published (Hogan 2013, submitted to ACM Trans. Math. Softw.), and the software has been released at
Very recently we have been working on a parallelized version of the Jacobian computation in Adept, using OpenMP; this exploits the multi-core architecture of many modern workstations. This will speed up the retrieval substantially. It will also make it more robust, as it would be coupled with a move to the Levenberg-Marquardt minimization scheme, which uses the curvature of the cost function in addition to its gradient which makes it easier to find the minimum of the cost function.
4.3.Storage of error descriptors
In previous projects there was considerable discussion about how to compress the information on errors, which for N retrieved variables is in the form of an NxN error covariance matrix and an NxN averaging kernel matrix. We now routinely report
(a)1-sigma errors; these are the square-root of the diagonal of the error covariance matrix
(b)Error correlations between variables at the same height; these are derived by normalizing specific off-diagonal elements of the error covariance matrix
(c)The area of the averaging kernel for a particular variable; this is found by integrating a row of the averaging kernel and indicates what fraction of the information for a particular retrieved variable is from the observations, as opposed to the a-priori constraints
(d)The spatial width of the peak in the averaging kernel; this indicates the extent to which the nature of the observations have smoothed the retrieval at each height.
An example for a simulated ice cloud profile retrieved by radar and lidar is shown in Figure 4. In this example, the radar and lidar both sampled the cloud down to 4 km, below which it was sampled by radar alone. It can be seen that the error rises rapidly in the radar-only region, particularly for the number concentration parameter N0’ which transitions back to the a-priori value when only one instrument is available. The errors in extinction and ice water content are found to be well correlated, but less so between extinction and either particle size or number concentration. The right panel shows that down to 4 km, the averaging kernel area is close to 1, indicating that most of the information comes from the observations, but below this height, where the lidar no longer detects the cloud, the retrieval of number concentration parameter N0’ in particular is forced to rely more on the a-priori estimate so the area of the averaging kernel falls.
Figure 3. Example of error descriptors for a retrieved ice cloud profile: (left) vertical profile of percentage errors in various variables; (centre) error correlation between various pairs of variables; (right) area of averaging kernel for two variables.4.4.New model for radar scattering by snowflakes
At the start of the VARSY project, the radar scattering model by ice particles was oblate spheroids, based on the findings of Hogan et al. (2012, J. Appl. Meteorol. Soc.).However, this model works only up to sizes equal to the wavelength. For larger particles, the internal structure of the ice particle, which is ignored in the oblate spheroid model, becomes important. This presents a problem because all ice particles are different and it is very time consuming to perform “discrete dipole approximation” calculations on a large number of representative ice particle shapes.
We have made a breakthrough in characterizing ice particle shape that allows a formula to be derived for the mean backscatter cross section of an ensemble of particles of the same size but different internal morphologies. This makes use of the Rayleigh-Gans scattering approximation, which recent work shows to be valid for low-density snowflakes at millimetre wavelengths, where the shape of the particle observed by a vertically pointing radar only needs to be described by a function A(z), which is the area of ice in a horizontal plane through the particle at height z.
The parameters of the model have been derived by examining 50 ice aggregates generated by the Westbrook aggregation model. One such aggregate is depicted in Figure 3a. The solid black line in Figure 3b shows the mean A(z) function, and a simple parametric fit by the dashed line. The smaller-scale features in the A(z) function have been found to follow a power law that may be characterized by two coefficients. This enables an analytic expression to be found for the backscatter cross section.
Figure 4. (left) Example ice aggregate from the Westbrook aggregation model, viewed with its longest axis horizontal (representing the ordinary fall mode of ice particles), where the shading is proportional to how much mass is present integrated along the third dimension. (right) The area of ice A(z) intersected by a horizontal plane through ice aggregates as a function of vertical distance z, where the mean area of 50 aggregates is shown by the solid black line, and a parametric fit is shown by the dashed line. The dotted line shows the specific case of the particle shown in the left panel.Figure 4 depicts the mean backscatter cross section from both the ensemble of aggregates (using the Rayleigh Gans approximation) and the new equation. It can be seen that the new equation generally fits very well, except at the far right of the plot corresponding to very high frequencies. At this point, the wavelength is smaller than the individual crystals that make up the aggregate, so the power law breaks down. This situation does not occur for ice aggregates at 94 GHz, so the new formula is suitable to use for both CloudSat and EarthCARE. To illustrate the importance of the effect of internal structure in the particles, consider a 1cm particle at 94 GHz, which corresponds to kDz=20 in Figure 4. It can be seen that here the backscatter of the aggregates is around 2 orders of magnitude larger than the backscatter of a homogeneous spheroid with the same mass and size.
This new scattering model has been implemented into the retrieval algorithm and it has been found that it approximately halves the ice water content and extinction coefficient retrieved from the radar alone in regions of high radar reflectivity factor, such as in the first few kilometres above the melting layer.
Figure 5. Mean backscatter cross section versus normalized wavenumber (where k=2/ and is the wavelength) for the 50 aggregates (solid black line) and the new equation (dashed line).4.5.Retrieval of riming factor
In principle, the EarthCARE Doppler velocity provides information on the density of ice particles, since particles with the same ice mass (and hence radar reflectivity) will fall faster if that mass is concentrated in a smaller volume. The density of ice particles is largely controlled by riming, i.e. accretion of supercooled water droplets in mixed-phase clouds. With reasonable observational support, the algorithm currently assumes the Brown & Francis mass-size relationship for ice clouds, applicable when riming is not present. Figure 5 depicts fall speed as a function of mass and size, and it can be seen that even for large particles, the Brown & Francis relationship corresponds to fall speeds no higher than around 1 m/s.
We have added the capability to retrieve a “riming factor”, which scales the exponent of the mass-size relationship such that a riming factor of 0 corresponds to Brown & Francis while a value of 1 corresponds to solid ice. It can be seen that for larger particles the fall speed, and hence Doppler velocity, increases monotonically with riming factor, indicating that riming factor can be retrieved if Doppler velocity is available.
Figure 6 shows a simulated retrieval of riming factor, where the top row shows the “truth”, the bottom row shows the simulated observations and the red dashed lines compare the retrieval to the truth and the simulated observations. In this case the riming factor of up to 0.5 is successfully retrieved. Work is still required to evaluate this approach, and in particular its applicability in that often riming occurs when there are significant vertical motions present, in which case the Doppler velocity could not be interpreted unambiguously in terms of terminal fall speed.
Figure 6. The colours show the terminal fall speed of ice particles according to the Heymsfield and Westbrook model, as a function of their diameter and their ice fraction (the fraction of a circumscribing sphere if this diameter that contains ice). The black line shows the default Brown & Francis mass-size assumption. The blue lines show a parameterization for how ice fraction is allowed to vary in a retrieval representing riming, where the number shown is the “riming factor”.Figure 7. Simulated retrieval of riming factor. The blue lines in the top row show the “truth”. These have been used to forward model the simulated EarthCARE observations on the bottom row, with noise. The retrieval has been run, with the dashed red lines on the bottom panels showing the forward modelled observations at the final iteration, and the dashed red lines on the top panel showing the retrieved cloud/rain variables. The grey lines depict the individual iterations of the algorithm.
4.6.Retrieval of aerosol properties using a Kalman smoother
The HSRL offers the potential to unambiguously retrieve the vertical profile of extinction coefficient, without any assumptions on the nature of the particles. However, in the case of aerosol which has a very low optical depth, it is challenging to extract the signal from the noise. The approach taken in the lidar-only L2a algorithms is to average the ATLID data first, but for the ACM-CAP retrieval, this is not really compatible with the retrieval of the other atmospheric constituents which can be retrieved independently on a profile-by-profile basis. Therefore, the capability has been added to use a “Kalman smoother” for any retrieved variable; currently this is forward-only, meaning that each profile uses the retrieved values at the previous profile as an additional constraint, which has a smoothing forward-in-time only. Soon a forward-reverse Kalman smoother will be added, in which the retrieval is repeated in the reverse direction, at which time the retrieval is constrained by both the previous and subsequent profiles, leading to smoothing that is symmetric in time.