Cloud and precipitation microphysical retrievals from the EarthCARE Cloud Profiling Radar: the C-CLD product

The Earth Clouds, Aerosols and Radiation Explorer (EarthCARE) satellite mission is a joint endeavour developed by the European Space Agency (ESA) and the Japan Aerospace Exploration Agency (JAXA) and features a 94 GHz Doppler Cloud Profiling Radar. This paper presents the theoretical basis of the cloud and precipitation microphysics (C-CLD) EarthCARE Level 2 (L2) algorithm. The C-CLD algorithm provides the best estimates of the vertical profiles of water mass content and hydrometeor characteristic size, obtained from radar reflectivity, path-integrated signal attenuation and hydrometeor sedimentation Doppler velocity estimates using optimal estimation (OE) theory. To obtain the forward model relations and the associated uncertainty, an ensemble-based method is used. This ensemble consists of a collection of in situ measured drop size distributions that cover natural microphysical variability. The ensemble mean and standard deviation represent the forward model relations and their microphysics-based uncertainty. The output variables are provided on the joint standard grid horizontal and EarthCARE Level 1b (L1b) vertical grid (1 km along track and 100 m vertically). The OE framework is not applied to liquid-only clouds in drizzle-free and lightly drizzling conditions, where a more statistical approach is preferred.

1 Introduction

Clouds and precipitation systems play a critical role in Earth's energy and hydrological cycle (Stephens et al., 2010, 2012). The accurate representation of cloud and precipitation systems in numerical models is essential for improving the predictability of weather and climate models. While surface-based observatories (Illingworth et al., 2007; Mather and Voyles, 2013; Kollias et al., 2020) can provide high-resolution observations suitable for process studies, satellite-based active remote sensors have the potential to obtain global estimates of cloud and precipitation microphysics and dynamics (Battaglia et al., 2020b). The National Aeronautics and Space Administration (NASA) A-Train satellite constellation (Stephens et al., 2002, 2018) first demonstrated the potential of active remote sensing from space. The Earth Clouds, Aerosols and Radiation Explorer (EarthCARE) mission (Wehr et al., 2023) scheduled for launch in 2024 features the first spaceborne Cloud Profiling Radar (CPR) with Doppler capability (Illingworth et al., 2015; Kollias et al., 2022a, b). The EarthCARE CPR observations will offer a unique opportunity for the collection of a global dataset of vertical motions and microphysics in clouds and precipitation.

Compared to CloudSat, the EarthCARE CPR has higher sensitivity (5–6 dB more sensitive), better vertical sampling (100 versus 240 m), higher along-track resolution (500 versus 1100 m) and a smaller instantaneous field of view (IFOV, 800 versus 1400 m) and includes Doppler velocity measurements and improved detection in the lowest atmosphere levels (Illingworth et al., 2015; Battaglia et al., 2020b; Burns et al., 2016; Lamer et al., 2019; Kollias et al., 2014). Based on these characteristics, the EarthCARE CPR is expected to provide an improved set of CPR observables, i.e. radar reflectivity, path-integrated attenuation (PIA) and Doppler velocity, which after their post-processing and quality control by the CPR processor (C-PRO) algorithms (Kollias et al., 2022b) will be used for the development of the CPR-only cloud and precipitation microphysics (C-CLD) product.

The long record of CloudSat observations and the parallel development and validation of the CloudSat data products provide a strong heritage for the C-CLD algorithm development. In particular, the use of the CloudSat CPR path-integrated attenuation (PIA in dB) for estimating the total liquid water path (LWP) in the atmospheric column and for constraining surface and profile estimates of the rainfall rate (Haynes et al., 2009; Lebsock and L'Ecuyer, 2011) is applied in a similar manner in C-CLD. Another factor that influenced the C-CLD algorithm development is the development of sophisticated ground-based networks such as the U.S. Department of Energy Atmospheric Radiation Measurement (ARM) observatories and the pan-European Aerosol, Clouds and Trace Gases Research Infrastructure (ACTRIS) (Illingworth et al., 2007; Mather and Voyles, 2013; Kollias et al., 2020). The measurements from these surface-based networks have stimulated the development of several algorithms that utilize the combination of radar reflectivity and mean Doppler velocity (Delanoë et al., 2007; Heymsfield et al., 2008; Mason et al., 2018; Oue et al., 2019). These efforts highlighted the information content of the Doppler velocity that is a new EarthCARE CPR observable from space compared to CloudSat. In addition, the C-CLD algorithm utilizes our latest understanding of solid hydrometeors scattering at 94 GHz (Hogan and Westbrook, 2014; Kneifel et al., 2020) and the availability of extensive ground-based observations of particle size distributions (Williams, 2012; von Lerber et al., 2017).

2 Description of the algorithm

The C-CLD cloud and precipitation retrieval algorithm is based on a profile-by-profile approach. At each profile, it uses information available from the radar-only measurements provided in the form of the following products, as described by Kollias et al. (2022b): the CPR feature mask and radar reflectivity (C-FMR), CPR cloud Doppler (C-CD) parameters, and CPR target classification (C-TC). The C-CLD algorithm derives the best estimates of cloud and precipitation microphysics that feed into the composite cloud and aerosol profiles product (ACM-COM, ATmospheric LIDar–Cloud Profiling Radar–MultiSpectral Instrument composite; Cole et al., 2022) as explained in Eisinger et al. (2023). The main retrieved quantities consist of the water mass content and particle characteristic size. First, the output of the C-TC hydrometeor classification is used to determine the occurrence of the specific hydrometeor type (ice cloud, snow, rimed snow, melting snow, cold rain, warm rain, non-drizzling liquid cloud, drizzling liquid cloud). This information is used to determine which branch of the C-CLD retrieval will be employed.

As illustrated in the flowchart shown in Fig. 1, the C-CLD processor contains specific algorithms designed to retrieve distinct cloud system types:

• a.

  liquid cloud retrieval with the separation between non-drizzling and drizzling liquid clouds;

• b.

  ice cloud retrieval and precipitation retrieval with specific algorithms designed to retrieve ice cloud, snow, riming snow, cold rain and warm rain.

The optimal estimation (OE) variational approach is applied as described in Rodgers (2000). It is based on a Gauss–Newton minimization algorithm that allows for a quantitative evaluation of the uncertainty in the retrieved quantities. The forward model within the OE approach maps two moments of the particle size distribution (PSD), the particle characteristic size (D_m; i.e. mean-mass-weighted equivalent melted diameter) and mass water content (MC) to the CPR reflectivity and Doppler velocity (previously corrected in C-CD for vertical air motion).

Figure 1 Cloud and precipitation retrieval scheme flow chart. OE theory is applied in the retrieval of solid and liquid precipitating clouds. The retrievals in the liquid clouds and lightly drizzling clouds (yellow box) are not performed using the OE method. HB: Hitschfeld and Bordan (1954).

The OE method is not applied for retrieval of drizzle-free clouds (i.e. non-precipitating liquid clouds) and for lightly drizzling clouds (group a). This is justified by the following facts.

• In the case of drizzle-free clouds, observed Doppler velocity does not provide any relevant information (fall velocity of the cloud droplet is negligible), so the reflectivity is the only measurement available.

• In the case of lightly drizzling clouds, the observed Doppler velocity could be heavily dominated by vertical air motion, leading to a large uncertainty in the reflectivity-weighted velocity. Moreover, the observed reflectivity is, in general, dominated by drizzle.

Therefore, a retrieval approach involving the optimal estimation method and much simpler methods based on the use of power-law relationships will result in similar uncertainties in microphysical parameters. The main retrieved variables for liquid clouds and drizzle are liquid water content and particle effective diameter, i.e. the ratio of the third and the second PSD moments. Note that the case of heavy drizzle is included in the rain retrieval, as a subcategory of warm rain.

The individual algorithms will now be described in detail in the next sections.

2.1 Liquid cloud and light-drizzle retrieval

Two distinct situations are analysed: non-precipitating liquid clouds and clouds that generate light drizzle. These two scenarios do not overlap spatially with each other or with pixels classified as precipitation, where the OE algorithm is used. When a radar volume is classified as a liquid cloud or light drizzle, an appropriate power-law formula is utilized to estimate the combined liquid water content of the cloud and any drizzle that may be present. In the case of the OE retrieval, the liquid water path is one of the retrieved state vector unknowns. Once retrieved, it is distributed adiabatically in the column to provide an estimate of the liquid water content. The algorithm clearly distinguishes between the cloud water and the precipitation mass content, with the latter being equal to 0 for liquid clouds and light drizzle.

2.1.1 Drizzle-free clouds

The cloud liquid water content (LWC) vertical structure is determined from the reflectivity values using the relationship LWC−Z_e, derived in a power-law form:

$$\begin{array}{}\text{(1)}& \text{LWC}\left(z\right)=\langle A\rangle Z_e(z{)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where 〈A〉 is an average value assumed constant across the whole height (Frisch et al., 1998). This relationship assumes that both the cloud droplets' number concentration and the PSD spectral width are constant with height. While this is reasonable in marine clouds (e.g. Miles et al., 2000), Löhnert and Crewell (2003) concluded that this assumption is the dominant error factor in continental clouds. If a measurable PIA signal is available, we can estimate the liquid water path (LWP) and subsequently estimate the parameter 〈A〉. Otherwise, two constant values of 4.7 and 2.4 g mm⁻³ m${}^{-\mathrm{3}/\mathrm{2}}$ are assumed over land and ocean, for 〈A〉, respectively.

The corresponding cloud LWP obtained from such an approximation is evaluated against the estimated cloud mass content adiabatic profile that represents an average in-cloud profile, independent of the reflectivity vertical variability with LWC increasing with the distance from the cloud base, as

$$\begin{array}{}\text{(2)}& {\text{LWC}}_{\mathrm{ad}}\left(z\right)=f_{\mathrm{wN}}{\mathrm{\Gamma }}_{\mathrm{w}}\left[{\displaystyle \frac{z_{\mathrm{0}\mathrm{w}}}{z_{\mathrm{0}\mathrm{w}}+(z-h_{\mathrm{cb}})}}\right](z-h_{\mathrm{cb}}),\end{array}$$

where the term in the square brackets accounts for a decrease in adiabaticity in the thicker clouds as proposed in Wood et al. (2009), z_0w is a scaling parameter set to 500 m, Γ_w is the average vertical gradient of the change in adiabatic LWC (parameterized in Eq. B1), h_cb is the height of the cloud base and f_wN is a normalization factor that is set to 1. In the case of unreliable LWP estimates (i.e. when there is more than one cloud layer and the PIA corresponding to the cloud is smaller than 2 dB), a minimum and a maximum limit of the adiabatic profile (0.3LWC_ad and 0.9LWC_ad, respectively) is enforced on the estimate from the power-law formula.

The cloud effective radius is computed as a mean between the relationship proposed by Fox and Illingworth (1997),

$$\begin{array}{}\text{(3)}& r_{\mathrm{eff}}=\mathrm{23.3}{Z}^{\mathrm{0.177}},\end{array}$$

and a relation derived for a lognormal PSD with a spectra width of 0.38 µm reported by Miles et al. (2000) as an average value,

$$\begin{array}{}\text{(4)}& r_{\mathrm{eff}}=\mathrm{46.5}\sqrt[\mathrm{6}]{Z/N_{\mathrm{cl}}}.\end{array}$$

In Eqs. (3) and (4), r_eff is in micrometres; Z is in mm⁶ m⁻³; and the number concentration N_cl is per cubic centimetres and is assumed to be equal to 288 and 74 over land and ocean, respectively (Miles et al., 2000).

2.1.2 Lightly drizzling clouds

The retrieval of LWC for lightly drizzling clouds combines two estimates:

1. the LWC derived from reflectivity based on power laws derived by Sauvageot and Omar (1987) and Baedi et al. (2000),

   $$\begin{array}{}\text{(5)}& {\displaystyle }\begin{array}{rl}& {\text{LWC}}_{\mathrm{SO}\mathrm{87}}\left[\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}\right]\\ & =\mathrm{12.25}{\left(Z\left[{\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}\right]\right)}^{\mathrm{0.763}}\text{for}Z<-\mathrm{22}\mathrm{dBZ},\end{array}\text{(6)}& {\displaystyle }\begin{array}{rl}& {\text{LWC}}_{\mathrm{B}\mathrm{00}}\left[\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}\right]\\ & =\mathrm{0.457}{\left(Z\left[{\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}\right]\right)}^{\mathrm{0.193}}\text{for}Z>-\mathrm{15}\mathrm{dBZ},\end{array}\end{array}$$

   and a linear interpolation between these two for the intermediary reflectivity regime (−22 to −15 dBZ);

2. the LWC profile derived from the adiabatic model (Eq. 2) with f_wN set to fit LWP estimates, if present, or otherwise set to 1.

The final LWCs are computed by combining the LWC_ad derived from the adiabatic model and LWC_Z from the reflectivity profile. The used LWC–Z relation for drizzling conditions represents an average relation and can introduce a large bias, mainly close to the cloud base and cloud top; in such regions, large differences between LWC_ad and LWC_Z are expected. Therefore, in the calculation of the LWC the weight attached to LWC_Z is progressively reduced where the absolute relative difference between LWC_ad and LWC_Z becomes large.

Note that the estimate of the liquid water content reported here includes both cloud and drizzle water content.

2.2 Optimal estimation framework components

For solid/liquid precipitation and ice clouds, the C-CLD algorithm applies a variational approach (Rodgers, 2000). It assimilates radar measurements and aims at balancing these data with the prior information to provide an optimal estimate of the state vector. Gauss–Newton iterations are used to find the best solution, which allows for a quantitative evaluation of the uncertainty in the retrieved quantities. This approach has been applied to similar radar-based microphysical retrievals in the past years (Lebsock and L'Ecuyer, 2011; Szyrmer et al., 2012; Battaglia et al., 2016, 2020c; Tridon et al., 2019a; Mason et al., 2017) and in the EarthCARE synergistic microphysical retrieval product (ACM-CAP, ATmospheric LIDar–Cloud Profiling Radar–MultiSpectral Instrument CAPTIVATE; Mason et al., 2022).

Radar measurements depend on a number of microphysical properties of hydrometeors in the sampled volume, including the particle size distribution (PSD) and, for solid-phase particles, the shape and mass distribution. Assumptions on any of the parameters listed above lead to a variety of microphysical relations reported in the literature between radar observables and microphysical properties (e.g. Protat et al., 2007; Matrosov and Turner, 2018). To assess the effect of the uncertainty associated with the microphysical description, the ensemble-based method is used to obtain the forward model relations and the associated simulation uncertainty. The ensemble consists of a number of particle size distributions collected at the ground for the Global Precipitation Measurement (GPM) mission ground validation programme (Dolan et al., 2018). Scattering models are applied to these data to map the microphysical quantities to the radar observables. The ensemble mean relations and its spread, defined as 1 standard deviation, represent the forward model relations and their microphysics associated uncertainty, respectively.

2.2.1 State vector

The PSD is parameterized using the concept of double-moment normalization. Following Delanoë et al. (2005), the normalizing moments are defined as

$$\begin{array}{}\text{(7)}& M_p=\underset{\mathrm{0}}{\overset{D_{max}}{\int }}{D}^{p}N\left(D\right)\mathrm{d}D,\end{array}$$

where p is the moment, D is the equivalent liquid sphere diameter and D_max is the diameter of the largest particle. The ratio of the fourth and the third moment represents the mean-mass-weighted melted diameter $D_{\mathrm{m}}=M_{\mathrm{4}}/M_{\mathrm{3}}$ and is used as the size scaling parameter, while M₃ is proportional to the water-equivalent mass content $\text{MC}=(\mathit{\pi }/\mathrm{6}){\mathit{\rho }}_{\mathrm{w}}{M}_{\mathrm{3}}$ that controls the magnitude of the PSD. D_max is set to be equal to 5 mm or 2.5D_m, whichever value is smaller. By selecting M₃ and M₄ moments, the PSD can be expressed as

$$\begin{array}{}\text{(8)}& N(D;D_{\mathrm{m}},\text{MC})=\text{MC}\times f(D;D_{\mathrm{m}}),\end{array}$$

where f represents functional forms that are reported in the literature, derived from large datasets for each given hydrometeor category.

The goal of the C-CLD algorithm is to retrieve two moments of the PSD, i.e. the mass content, MC, and D_m, from radar reflectivity and Doppler velocity measurements. Therefore, the vector of the retrieval unknowns has the following form:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathit{x}=[& {\log }_{\mathrm{10}}{D}_{\mathrm{m}}^{\mathrm{1}},{\log }_{\mathrm{10}}{D}_{\mathrm{m}}^{\mathrm{2}},\mathrm{\dots }{\log }_{\mathrm{10}}{D}_{\mathrm{m}}^N,{\log }_{\mathrm{10}}{\text{MC}}^{\mathrm{1}},\\ & {\log }_{\mathrm{10}}{\text{MC}}^{\mathrm{2}},\mathrm{\dots }{\log }_{\mathrm{10}}{\text{MC}}^N{]}^T,\end{array}\end{array}$$

where N is the total number of the CPR range gates, regardless of whether they include ice/snow or rain. It is important to note that we do not use a separate notation for the mass content of ice and rain in this study, although it is commonly referred to as “ice water content” (IWC) and “rain water content” (RWC) in the literature. In case of warm- and cold-rain retrieval, the vector x also includes the cloud liquid water path with liquid water content distributed according to Eq. (2). In cold rain, the attenuation of the melting layer is an additional unknown. Note that the errors in the variables in the logarithmic units can be converted to fractional errors in the variable in the linear scale by the following error propagation formula:

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{\Delta }z}{z}}\approx {\displaystyle \frac{z^{\prime }\mathrm{\Delta }x}{z}}={\displaystyle \frac{\ln \mathrm{10}z\mathrm{\Delta }x}{z}}=\ln \mathrm{10}\mathrm{\Delta }x,\end{array}$$

where z=10^x and x is either log ₁₀D_m or log ₁₀MC, e.g. the root mean square error in log ₁₀z of 0.3 corresponds to the fractional error of 69 % in z.

2.2.2 Vector of measurements

The forward model maps the retrieved microphysical parameters to the space of radar measurements (attenuated radar reflectivity Z_m; Doppler velocity corrected for air motion U_D; and, in some cases, path-integrated attenuation).

The equivalent reflectivity factor for a radar operating at the wavelength λ is given by

$$\begin{array}{}\text{(11)}& Z_e={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{5}}|K_{\mathrm{w}}{|}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathit{\sigma }}_{\mathrm{b}}(D,\mathit{\lambda })N\left(D\right)\mathrm{d}D,\end{array}$$

where σ_b is the backscattering cross-section of a particle and K_w is the dielectric factor of liquid water at a reference temperature and frequency. For this study, K_w is assumed to be equal to 3.195+1.667i, which represents its value at 10 ^∘C according to the model of Turner et al. (2016). The reflectivity is usually expressed in mm⁶ m⁻³ or, due to its high variability, in logarithmic units of dBZ $=\mathrm{10}{\log }_{\mathrm{10}}\left({\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}\right)$.

When analysing millimetre wavelength radar data, the attenuation due to gases (mainly water vapour and oxygen) and the one caused by hydrometeors cannot be neglected (Battaglia et al., 2020b; Tridon et al., 2020; Lamer et al., 2019). The measured reflectivity at distance r from the radar is given by $Z_{\mathrm{m}}\left(r\right)=Z_e\left(r\right)\exp \left[-\mathrm{0.2}\ln \mathrm{10}{\int }_{\mathrm{0}}^{r}k\left(s\right)\mathrm{d}s\right]$ or, in more commonly used logarithmic units, $Z_{\mathrm{m}}\left(r\right)\left[\mathrm{dBZ}\right]=Z_e\left(r\right)\left[\mathrm{dBZ}\right]-\mathrm{2}{\int }_{\mathrm{0}}^{r}k\left(s\right)\mathrm{d}s$, where k is the so-called specific attenuation given in decibels per unit length; its component associated with the hydrometeors can be computed as the extinction cross-section-weighted (σ_e) integral of the PSD:

$$\begin{array}{}\text{(12)}& k={\displaystyle \frac{\mathrm{10}}{\ln \mathrm{10}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathit{\sigma }}_{\mathrm{e}}\left(D\right)N\left(D\right)\mathrm{d}D.\end{array}$$

Over waterbodies, the total path-integrated attenuation (PIA $\equiv \mathrm{2}\int k\left(s\right)\mathrm{d}s$) can be estimated from the surface return, and then it is used as an additional observational constraint.

The mean Doppler velocity is the backscattering-weighted line-of-sight velocity (v_LOS) of targets relative to the radar:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{U}_D& =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}N\left(D\right){\mathit{\sigma }}_{\mathrm{b}}(D,\mathit{\lambda })v_{\mathrm{LOS}}\left(D\right)\mathrm{d}D\\ & \times {\left(\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}N\left(D\right){\mathit{\sigma }}_{\mathrm{b}}(D,\mathit{\lambda })\mathrm{d}D\right)}^{-\mathrm{1}}.\end{array}\end{array}$$

Here, positive velocities correspond to downward motions (away from the CPR).

The CPR processor (C-PRO; Kollias et al., 2022b) derives an estimate of the CPR measurements with their associated uncertainties. This includes the attenuated radar reflectivity, the PIA provided by the C-FMR product and the sedimentation Doppler velocity (C-CD product). The estimation of the sedimentation velocity from raw EarthCARE CPR Doppler velocity measurements is a multistep, complex process consisting of non-uniform beam-filling correction; velocity unfolding; spatial averaging; and finally the sedimentation velocity estimate, where the contribution of the vertical air motion has been removed (based on the methodology of Kalesse and Kollias, 2013). The C-CLD algorithm takes in radar reflectivity, sedimentation velocity and PIA measurements as inputs. The measurement vector is composed of 2N+1 entries, and it is given by

$$\begin{array}{}\text{(14)}& \mathit{y}={\left[Z_{\mathrm{m}}^{\mathrm{1}},Z_{\mathrm{m}}^{\mathrm{2}},\mathrm{\dots }{Z}_{\mathrm{m}}^N,U_D^{\mathrm{1}},U_D^{\mathrm{2}},\mathrm{\dots }{U}_D^N,\text{PIA}\right]}^T,\end{array}$$

where N is the number of retrieval layers. However, because the normalized radar cross-section of the surface over land varies widely depending on factors such as vegetation, surface slope, soil moisture and snow cover, estimates of PIA are only provided over the ocean. The retrieval is still performed without PIA estimates, but results are significantly more uncertain and should be used with caution.

The vertical resolution of the retrieval matches the radar sampling, and it is equal to 100 m. Note that the actual vertical resolution of the radar is 500 m, which implies a factor of 5 oversampling. Thanks to a large antenna (2.5 m) and low aircraft altitude (400 km), the CPR is expected to achieve an unprecedented result in space sensitivity and collect measurements as low as −36 dBZ.

2.2.3 OE procedure

The aim of the OE is to provide the most probable value of the microphysical state vector x given the information provided by the measurements and prior knowledge about the state of the atmosphere. This is done by an iterative search that minimizes the cost function Φ:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}\mathrm{\Phi }& ={\left[\mathit{y}-F\left(\mathit{x}\right)\right]}^T{\mathbf{R}}_y^{-\mathrm{1}}\left[\mathit{y}-F\left(\mathit{x}\right)\right]\\ & +{\left[\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right]}^T{\mathbf{R}}_{\mathrm{a}}^{-\mathrm{1}}\left[\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right]+{\mathit{x}}^T{\mathbf{R}}_{\mathrm{TT}}\mathit{x}.\end{array}\end{array}$$

Here, F denotes the forward model (radar simulator), ${\mathbf{R}}_{\mathrm{a}}={\mathbf{R}}_{\mathrm{m}}+{\mathbf{R}}_{\mathrm{F}}$ represents the sum of the measurement error covariance matrix R_m and the forward model error covariance matrix R_F, and R_a represents the prior covariance matrix. The measurement errors are assumed to be uncorrelated, and so the matrix R_m is diagonal. The reflectivity and Doppler velocity errors depend mainly on the number of independent samples and on the signal-to-noise ratio (SNR); for a typical measurement, they are 1 dB and 0.2 m s⁻¹, respectively. An uncertainty estimation of the PIA is more complex (e.g. it depends on the surface characteristics within the radar field of view), but it is provided by the C-PRO (for more detail on the PIA estimator, see Kollias et al., 2022b). A smoothness constraint is introduced in the form of the Twomey–Tikhonov matrix R_TT, including a scaling coefficient, as described in Hogan (2006). The state vector at the ith iteration can be found using Gauss–Newton minimization steps, i.e.

$$\begin{array}{}\text{(16)}& {\displaystyle }\begin{array}{rl}{\mathit{x}}_i=& {\mathit{x}}_{i-\mathrm{1}}+{\mathbf{H}}^{-\mathrm{1}}[{\mathbf{J}}^T{\mathbf{R}}_{\mathrm{a}}^{-\mathrm{1}}\left(\mathit{y}-F\left({\mathit{x}}_{i-\mathrm{1}}\right)\right)\\ & -{\mathbf{R}}_{\mathrm{a}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{i-\mathrm{1}}\right)-{\mathbf{R}}_{\mathrm{TT}}{\mathit{x}}_{i-\mathrm{1}}{]}^T,\end{array}\text{(17)}& {\displaystyle }\mathbf{H}\equiv {\mathbf{J}}^T{\mathbf{R}}_y^{-\mathrm{1}}\mathbf{J}+{\mathbf{R}}_{\mathrm{a}}^{-\mathrm{1}}+{\mathbf{R}}_{\mathrm{TT}},\end{array}$$

where J denotes the Jacobian (gradient) of the forward model F. Usually, after a few iterations, the algorithm converges to the minimum that provides the final solution. If convergence criteria are not met within a set number of iterations, the state variables are set to the missing value.

The advantage of using the OE approach is that it provides a method for propagating errors in the measurements and uncertainties in the algorithm assumptions. The error covariance matrix R_x associated with the retrieval variables is given as

$$\begin{array}{}\text{(18)}& {\mathbf{R}}_x^{-\mathrm{1}}={\mathbf{J}}^T{\mathbf{R}}_y^{-\mathrm{1}}\mathbf{J}+{\mathbf{R}}_{\mathrm{a}}^{-\mathrm{1}}.\end{array}$$

The diagonal elements of R_x provide the estimates of the variance of x, i.e. decimal logarithm of the retrieved quantities (MC and D_m). The off-diagonal elements give the cross-correlations between errors. The errors for any related quantity, like precipitation rate, can be computed by propagating these errors.

2.3 Warm rain

2.3.1 Forward model of rain reflectivity, attenuation and Doppler velocities

For rain, a gamma model is used to analytically approximate the PSD shape, i.e. the function f in Eq. (8) is

$$\begin{array}{}\text{(19)}& \begin{array}{rl}f(D;D_{\mathrm{m}},\mathit{\mu })=& {\displaystyle \frac{\mathrm{6}(\mathit{\mu }+\mathrm{4}{)}^{\mathit{\mu }+\mathrm{4}}}{\mathit{\pi }{\mathit{\rho }}_{\mathrm{w}}{D}_{\mathrm{m}}^{\mathrm{4}}\mathrm{\Gamma }(\mathrm{4}+\mathit{\mu })}}{\left({\displaystyle \frac{D}{D_{\mathrm{m}}}}\right)}^{\mathit{\mu }}\\ & \exp \left(-{\displaystyle \frac{D(\mathit{\mu }+\mathrm{4})}{D_{\mathrm{m}}}}\right),\end{array}\end{array}$$

where Γ denotes the gamma function and μ is a shape-controlling parameter. Schulte et al. (2022) have demonstrated that, in warm-rain retrievals, single-moment PSD models can lead to large biases, of the order of 100 %, when retrieving rain rates. The selection of the shape parameter μ is based on the methodology presented by Williams et al. (2014), where the expected value of μ is found for a given D_m, based on the statistical analysis of in situ microphysical measurements. In this study, in situ PSD data collected during field campaigns and from the permanent sites of the ground validation programme of the Global Precipitation Measurement mission (GPM; Hou et al., 2014) are exploited (for more detail, see Mróz et al., 2019). The analysis is restricted to the measurements from the two-dimensional video disdrometer (2DVD; Kruger and Krajewski, 2002) with a series of quality checks performed beforehand. These checks include discarding frozen precipitation or insufficient PSD sampling that happens for small rainfall rates (≤0.1 mm h⁻¹) and large sizes (D_m≥4 mm) that disdrometers are not well suited to capture (Guyot et al., 2019). These filtering criteria are set to have statistically and physically meaningful PSDs. The final dataset includes almost 150 000 samples of rainy measurements over different latitudes, thus thoroughly covering natural variability.

Our analysis confirmed previous findings of Williams et al. (2014) about the microphysical properties of PSDs (i.e. the mass-weighted standard deviation of D); the so-called PSD width (σ_m) is highly correlated with D_m, and its expected value is given by

$$\begin{array}{}\text{(20)}& {\mathit{\sigma }}_{\mathrm{m}}^{\mathrm{expected}}={\displaystyle \frac{D_{\mathrm{m}}^{\mathrm{1.4}}}{\sqrt{\mathrm{10}}}},\end{array}$$

where σ_m and D_m are in millimetres. Although these statistics are based on binned PSD measurements with no underlying assumptions about the PSD shape, they can be translated into a gamma-model-specific relation via ${\mathit{\sigma }}_{\mathrm{m}}^{\mathrm{\Gamma }}=D_{\mathrm{m}}/\sqrt{\mathrm{4}+\mathit{\mu }}$ giving

$$\begin{array}{}\text{(21)}& {\mathit{\mu }}^{\mathrm{expected}}=\mathrm{10}{D}_{\mathrm{m}}^{-\mathrm{0.8}}-\mathrm{4}.\end{array}$$

For the forward model simulation, backscattering and extinction cross-sections are computed with the T-matrix approximation assuming the axial-ratio formula of Brandes et al. (2005). The Doppler velocity is computed using the raindrop terminal fall speed as determined by Gunn and Kintzer (1949). It is assumed that the PSD shape can be parameterized by Eq. (19) with μ given by Eq. (21) to reduce the number of free parameters in the retrieval. The uncertainty in such an approximation was estimated via analysis of the radar simulations for the binned PSDs collected at the ground. The forward model errors for the reflectivity, specific attenuation and mean Doppler velocity are 0.42 dB, 10 % and 0.12 m s⁻¹, respectively. Note that the specific attenuation uncertainty is given in terms of a fractional error, as it strongly varies with the absolute value. The simulated radar observables corresponding to the in situ PSD measurements and the forward model parameterization used in this study can be found in Appendix A1.

2.3.2 Cloud liquid water correction

A crucial component in the warm-rain algorithm for W-band radars is the cloud liquid water correction (Haynes et al., 2009; Battaglia et al., 2020a). The following strategy has been followed: first, the cloud boundaries are identified based on the lifting condensation level (cloud bottom) and by the highest altitude of the detectable reflectivity (cloud top); then the shape of the profile of cloud MC given by Eq. (2) is attributed to the measurement column. Once the shape is fixed, the magnitude of the liquid water content is controlled by the cloud liquid water path (CLWP) that, in the logarithmic units, is one of the retrieved unknowns. Because the radar reflectivity of cloud droplets is much smaller than the one of raindrops, the retrieval of log ₁₀(CLWP) is mainly driven by the PIA estimate.

2.3.3 A priori

One of the essential elements of the OE procedure is the initial estimation of the microphysical parameter values along with their uncertainties. This can be done by providing climatological statistics based on long-term observations. This approach usually involves very large uncertainties that correspond to the natural variability in the rain microphysics. Alternatively, a much more constrained a priori estimate can be obtained by statistical analysis of in situ PSD measurements in relation to their radar simulations, as was done by Tridon et al. (2019b). For example, an estimate of the mean value and standard deviation of log ₁₀D_m and log ₁₀MC in correspondence to a given reflectivity range (Z±SD(Z)) can be provided. This approach is adopted in this study.

The a priori information on log ₁₀D_m and log ₁₀MC is obtained from the rain microphysics statistics, and their corresponding reflectivity simulations collected in the PSD dataset are described in Sect. 2.3.1. Regression analysis reveals a moderate correlation with the Pearson correlation coefficient (CC) of 0.53 between the state vector parameters via the following linear formula:

$$\begin{array}{}\text{(22)}& {\log }_{\mathrm{10}}\text{MC}\left[\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}\right]=\mathrm{1.863}{\log }_{\mathrm{10}}{D}_{\mathrm{m}}\left[\mathrm{cm}\right]+\mathrm{0.757}.\end{array}$$

The root mean square error (RMSE) of this fit is estimated to be 0.33 B for $\mathrm{0.8}<D_{\mathrm{m}}<\mathrm{2}$ mm. Since the PSD dataset does not include small raindrop sizes, we use regressions (21) and (22) together with the related uncertainties to supplement the in situ data with a low precipitation rate and/or low reflectivity points. This leads to the following a priori relations in rain:

$$\begin{array}{}\text{(23)}& {\displaystyle }\begin{array}{rl}& {\log }_{\mathrm{10}}{D}_{\mathrm{m}}\left[\mathrm{cm}\right]\\ & =\left\{\begin{array}{l}\mathrm{0.020}{Z}_e\left[\mathrm{dBZ}\right]-\mathrm{1.446}\text{if}{Z}_e<\mathrm{6.75}\mathrm{dBZ},\\ \mathrm{0.036}{Z}_e\left[\mathrm{dBZ}\right]-\mathrm{1.554}\text{if}\mathrm{6.75}\le Z_e<\mathrm{17}\mathrm{dBZ},\\ \mathrm{0.012}{Z}_e\left[\mathrm{dBZ}\right]-\mathrm{1.147}\text{if}{Z}_e\ge \mathrm{17}\mathrm{dBZ};\end{array}\right.\end{array}\text{(24)}& {\displaystyle }\begin{array}{rl}& {\log }_{\mathrm{10}}\text{MC}\left[\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}\right]\\ & =\left\{\begin{array}{l}\mathrm{0.038}{Z}_e\left[\mathrm{dBZ}\right]-\mathrm{2.043}\text{if}{Z}_e\le \mathrm{12.5}\mathrm{dBZ},\\ \mathrm{0.109}{Z}_e\left[\mathrm{dBZ}\right]-\mathrm{2.932}\text{if}{Z}_e>\mathrm{12.5}\mathrm{dBZ}.\end{array}\right.\end{array}\end{array}$$

Uncertainties in these relations over the whole range of reflectivity values are estimated to be 0.15 and 0.2 B, respectively (i.e. a factor of 1.41 and 1.58), representing the maximum RMSE value for PSD simulations partitioned into 1 dB reflectivity bins from −15 to 32 dBZ. Note that, for large reflectivity, the slope of the Z–D_m relation is very small compared to the uncertainty estimate, which indicates a weak correlation between these parameters. In practice, it reduces the Z–D_m relationship to the climatological value of log ₁₀D_m provided by the in situ dataset.

The derived regressions require effective reflectivity estimates; therefore, the radar measurements are corrected for attenuation by using the Hitschfeld and Bordan (1954) methodology before a priori estimates are derived.

The expected value of the cloud liquid water path is estimated to be weakly related to the rain water path (RWP), i.e.

$$\begin{array}{}\text{(25)}& {\log }_{\mathrm{10}}{\displaystyle \frac{\text{CLWP}}{\text{RWP}}}=\mathrm{0.344}\pm \mathrm{0.26}.\end{array}$$

This formula is based on the statistical analysis of warm-rain simulations over the Cabo Verde islands described in Sect. 3.

2.4 Ice and snow

Large natural variability in ice microphysics results in a variety of solid-phase hydrometeor structure models. In this study, the mass of the snowflakes is modelled using the parameterization of Morrison et al. (2009), where riming is simulated by filling the gaps between the ice crystal branches with supercooled liquid droplets (Heymsfield, 1982). The mass of snowflakes is parameterized by the power-law formula of m[kg]=α(D[m])^β, with α and β varying for different size regimes. For the unrimed aggregate, it is assumed that α=0.01 and β=2, agreeing with the simulations (e.g. Leinonen and Szyrmer, 2015; Westbrook et al., 2004) and in situ measurements of aggregates (Brown and Francis, 1995; Erfani and Mitchell, 2017; Moisseev et al., 2017). For sizes where the power-law formula would exceed the mass of solid ice spheres, the latter is used. In riming conditions, the smallest aggregates are fully filled with rime and they grow by accretion, so their mass–size relation follows the one for graupel (α=86.6, β=3). During riming, large aggregates do not increase their size due to the collection of supercooled droplets, but they only increase their mass proportionally to their projected surface area and the amount of supercooled liquid water the snowflake passes through. This implies that the exponent in the mass–size formula for partially rimed snow remains the same as for unrimed aggregates (β=2) and it is only α that increases with the degree of riming. It is implicitly assumed that the mass of rimed aggregates is always larger than or equal to the mass of unrimed snow; therefore the maximum between the power-law formulas for rimed and unrimed aggregates is taken. For more detail on this conceptual model, see Mroz et al. (2021b) and their Fig. 1, which shows the transition points between different mass–size relationship regimes. With this parameterization, a degree of riming is fully represented by the value of α that is equal to 0.01 for unrimed aggregates and reaches 0.5 for heavily rimed large graupel particles. The OE retrieval for snow profiles is performed for five different values of α, and the one that provides the lowest cost function (see Eq. 15) is used as a final state estimate.

2.4.1 Forward model of ice reflectivity, attenuation and Doppler velocities

The scattering properties of snow particles are obtained by using discrete dipole approximation corresponding to realistic snowflake shapes (see Leinonen et al., 2016). These snowflakes are composed of dendrites of different sizes, and they are subject to various degrees of riming. In the computations, the radar is pointing vertically, the particles are aerodynamically aligned with the maximum dimension oriented horizontally and particles are discretized to a collection of 40 µm dipoles. The original dataset of Leinonen et al. (2016) is complemented by large aggregates generated by the authors using the same aggregation model (https://github.com/jleinonen/aggregation, last access: 1 October 2022). The terminal velocity of particles is simulated for standard atmospheric conditions (relative humidity of 50 %, $\stackrel{\mathrm{̃}}{T}=\mathrm{20}$ ^∘C, $\stackrel{\mathrm{̃}}{P}=\mathrm{1013}$ mbar) using the parameterization of Böhm (1992). The physical and scattering properties of individual snowflakes are freely available at https://doi.org/10.5281/zenodo.7510186 (Mroz and Leinonen, 2023). The velocity U_D(p,T) at any temperature T and pressure p is computed via an air density correction as suggested by Foote and du Toit (1969):

$$\begin{array}{}\text{(26)}& U_D(p,T)=U_D(\stackrel{\mathrm{̃}}{p},\stackrel{\mathrm{̃}}{T}){\left[{\displaystyle \frac{\stackrel{\mathrm{̃}}{p}T}{p\stackrel{\mathrm{̃}}{T}}}\right]}^{\mathrm{0.4}}.\end{array}$$

Consistency between the microphysical parameterization and radar simulations is achieved by assuming that the scattering properties of snowflakes are functions of their mass and size only. For this purpose, for a selected mass–size formula (i.e. a selected degree of riming, α), the scattering database is searched for aggregates in the proximity of that relation. More specifically, for a given size D only snowflakes that satisfy ${\mathrm{\Delta }}_m=\mathrm{10}{\log }_{\mathrm{10}}\left|m_{\mathrm{snowflake}}\right(D)/m_{\mathrm{expected}}(D\left)\right|<\mathrm{3}$ are considered; i.e. the mass is within a factor of 2 of the formula. Next, depending on the distance from the expected mass–size relation, the particle is assigned its weight, $w\left(m\right)=\exp (-{\mathrm{\Delta }}_m^{\mathrm{2}})$. The scattering properties for a given mass (and α) are computed by locally fitting a fifth-degree polynomial to the decimal logarithm of the cross-sections as a function of log ₁₀m. The fitting of logarithmic values is adopted because of the large variability in the cross-sections with respect to the mass. Moreover, it reduces the variability in the averaged variables. The terminal velocities are fitted without the logarithmic transformation.

Once the snowflake density model is chosen and the corresponding scattering and falling velocity simulators are obtained, it is necessary to characterize the particle size distribution so that the description of the forward model of snow is complete. Due to the complexity of snow crystal shapes, the wide range of their densities, the ambiguities in the size definition (von Lerber et al., 2017) and the related difficulties in the PSD measurements, we decided not to use in situ snow PSD measurements to derive their statistical properties. Instead, it is assumed that the rain that was captured by the disdrometers of the GPM ground validation programme formed from snow melting and thus, by taking into account the differences in raindrop and snowflake terminal velocities, can be used to fully describe the natural variability in PSDs in snow. Implicitly, we assume that melting is the only process that occurs while snowflakes melt; i.e. no collision coalescence, breakup, condensation or evaporation takes place. By doing so, the particle size distributions in rain (N_r) and snow (N_s) are linked via the following relation:

$$\begin{array}{}\text{(27)}& \begin{array}{rl}& N_{\mathrm{s}}(D_{\mathrm{eq}},\mathit{\alpha })V_{\mathrm{s}}(D_{\mathrm{eq}},\mathit{\alpha })\mathrm{d}{D}_{\mathrm{eq}}=N_{\mathrm{r}}\left(D_{\mathrm{eq}}\right)V_{\mathrm{r}}\left(D_{\mathrm{eq}}\right)\mathrm{d}{D}_{\mathrm{eq}}\\ & \Rightarrow N_{\mathrm{s}}(D_{\mathrm{eq}},\mathit{\alpha })=N_{\mathrm{r}}\left(D_{\mathrm{eq}}\right){\displaystyle \frac{V_{\mathrm{r}}\left(D_{\mathrm{eq}}\right)}{V_{\mathrm{s}}(D_{\mathrm{eq}},\mathit{\alpha })}},\end{array}\end{array}$$

where V_r and V_s denote the terminal velocity in rain and snow, respectively, and D_eq is the equivalent melted diameter. The statistics about the microphysical properties of rain derived in Sect. 2.3.1 translate naturally, through melting-only assumption formula (27), into characteristics of snow. In particular, the PSD of snow, after melting, converts into the gamma PSD (Eq. 19) with $\mathit{\mu }=\mathrm{10}{D}_{\mathrm{m}}^{-\mathrm{0.8}}-\mathrm{4}$. The radar forward model is obtained by combining the electromagnetic and microphysical properties of snow. The scattering properties for selected values of α are shown in Fig. A2.

2.4.2 A priori

The a priori profiles of MC and D_m are generated using the empirical relations that take Z and temperature into account. Estimates of the mass content and a priori D_m are based on the relationships provided by Matrosov and Heymsfield (2008, 2017):

$$\begin{array}{}\text{(28)}& {\displaystyle }\text{MC}\left[\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}\right]=\mathrm{0.086}{\left(Z_e\left[{\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}\right]\right)}^{\mathrm{0.92}};\text{(29)}& {\displaystyle }{D}_{\mathrm{m}}\left[\mathrm{cm}\right]=\left\{\begin{array}{l}\mathrm{0.052}{\left(Z_e\left[\mathrm{dBZ}\right]\right)}^{\mathrm{0.28}}\text{for snow},\\ \mathrm{0.047}{\left(Z_e\left[\mathrm{dBZ}\right]\right)}^{\mathrm{0.294}}\text{for cirrus}.\end{array}\right.\end{array}$$

The reflectivity profiles are corrected for attenuation before the above relationships are applied. First, the cloud liquid water correction is performed. In the presence of riming, a constant amount of supercooled LWC (SLWC) is present across the ice layer for all pixels flagged as riming snow in the C-TC product. Attenuation is computed according to the parameterization provided in Sect. A2, and the reflectivity profile is corrected for the SLWC attenuation. Then the ice profile is further corrected for ice attenuation using the Hitschfeld and Bordan (1954) approach with the two-way attenuation coefficient proposed by Protat et al. (2019):

$$\begin{array}{}\text{(30)}& k_{\mathrm{ice}}\left[\mathrm{dB}{\mathrm{km}}^{-\mathrm{1}}\right]=\mathrm{0.0325}Z\left[{\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}\right].\end{array}$$

In the presence of a PIA measurement, if the attenuation is overestimated, the LWC is reduced to match the PIA. On the other hand, if the correction underestimates the PIA, the coefficient in the Z–k relation is scaled to match the PIA.

2.5 Cold rain

The cold-rain retrieval capitalizes on the modelling for the liquid phase described in Sect. 2.3 and on the solid phase described in Sect. 2.4. In cold rain, in the layer where temperatures become warmer than 0 ^∘C, hydrometeors transition between the solid and liquid phases. This region is very well identified by the target classification (C-TC). The modelling for Doppler velocities and reflectivities for the solid and the liquid phase follows what is described in Sects. 2.3.1 and 2.4.1. The melting layer is not modelled, and observables within the melting layer are not fitted as in Tridon et al. (2019b). The melting layer attenuation coefficient is estimated to be proportional to the mean rain rate underneath:

$$\begin{array}{}\text{(31)}& k_{\mathrm{ML}}\left[\mathrm{dB}{\mathrm{km}}^{-\mathrm{1}}\right]={\mathit{\gamma }}_{\mathrm{ML}}\langle \text{RR}\left[\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}}\right]{\rangle }^{{\mathit{\delta }}_{\mathrm{ML}}},\end{array}$$

with γ_ML=2.6 and δ_ML=0.87 as proposed by Matrosov (2008). This estimate is used as a soft constraint only; i.e. the bright band extinction is added to the vector of the unknown variables. During the OE iterations, the difference between its expected and state vector value is minimized, assuming the uncertainty in the Matrosov (2008) formula to be a factor of 2. The liquid cloud content (in logarithmic units) is also retrieved in cold rain. It is assumed that the liquid cloud is distributed between the freezing level and the height of the lifting condensation level (LCL) according to Eq. (2). Due to the high uncertainty regarding the occurrence of the cloud and its possible water content, it is assumed that the a priori estimate of the cloud water path (CWP) is very small, i.e. 0.1 g m⁻², which has no effect on the radar measurements. The relative uncertainty in this estimate is set to be 100 dB, which reflects no prior knowledge of this parameter.

Unlike the retrieval of snow profiles, the cold-rain retrieval is performed for one value of α only. Selection of the best α value is based on the continuity of the mass flux between the solid and the liquid phase, and it follows these steps: first, utilizing Eqs. (A1) and (A2), the mass content and the characteristic size of rain below the melting layer are estimated from the mean Doppler velocity and radar reflectivity measurements corrected for attenuation using the PIA-constrained Hitschfeld and Bordan (1954) technique. Once the water content and the size of rain are known, the radar simulations in the ice part are performed for logarithmically sampled values of α ranging from 0.01 to 0.5, assuming that the equivalent melted diameter D_m and the precipitation rate in rain and ice are the same; i.e. melting is the only process within the melting zone (Mróz et al., 2021a). Then, the difference between the radar simulations and the measurements in the radar bin above the melting zone is computed for all considered values of α, taking into account corresponding measurement uncertainties. Finally, an α value that minimizes this distance is selected for the retrieval.

3 Validation of the algorithm

The validation of the algorithm was performed with the synthetic precipitation scenes generated by Global Environmental Multiscale (GEM) model (Côté et al., 1998; Girard et al., 2014).

3.1 Warm rain

The C-CLD algorithm has been tested with warm-rain simulations over the Cabo Verde islands. The cloud microphysical processes were represented by the Predicted Particle Properties (P3) two-moment bulk microphysics scheme (Morrison and Milbrandt, 2015; Milbrandt et al., 2016). In the P3 scheme, the ice-phase hydrometeors are represented by three ice categories whose physical properties evolve continuously and were proved sufficient to represent the co-existence of cloud ice particles of different sizes (Qu et al., 2022). In addition to the three distinct ice species, rain and cloud droplets are also simulated. The horizontal resolution of the simulation is 250 m, which allows for resolving fine-scale convective cells that are characteristic of warm rain. The readers are referred to Qu et al. (2022) for more detail.

To simulate the radar measurements, the effective reflectivity and the specific attenuation of rain are estimated using formulas (A1) and (A3) in each model bin. The cloud contribution is simulated with an exponential PSD and summed up with the rain components. Then, the attenuated reflectivity, at the native model resolution, is computed by integrating the attenuation along the vertical path. The resulting three-dimensional reflectivity field is averaged horizontally over 3×3 pixels to provide a resolution of 0.75×0.75 km² that is comparable with the one of the EarthCARE CPR. Similarly, the mean Doppler velocity is first simulated at the native resolution. Next, it is averaged over 3×3 pixels using the attenuated reflectivity (in the linear units) as the weights. This provides the Doppler measurements at the radar scale. An estimate of the PIA aims to reflect as closely as possible the values that would be observed with the surface reference technique. The normalized radar cross-section, σ₀ [dB], is assumed to be uniform in the field of view. Then, the apparent PIA is given by

$$\begin{array}{}\text{(32)}& \begin{array}{rl}{\text{PIA}}_{\mathrm{SRT}}& ={\mathit{\sigma }}_{\mathrm{0}}-\mathrm{10}{\log }_{\mathrm{10}}\left({\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\mathrm{10}}^{\mathrm{0.1}({\mathit{\sigma }}_{\mathrm{0}}-{\text{PIA}}_i)}\right)\\ & =-\mathrm{10}{\log }_{\mathrm{10}}\left({\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\mathrm{10}}^{-\mathrm{0.1}{\text{PIA}}_i}\right),\end{array}\end{array}$$

where PIA_i denotes the path-integrated attenuation in the ith column and n=9 is the number of the spatially averaged profiles of the simulations. The water mass content, similarly to the reflectivity field, is averaged over nine neighbouring pixels. The characteristic size at the radar resolution is the mean of the fine-scale D_m values weighted by the corresponding mass content. Both rain and cloud components are taken into account in the state and measurement vector computations. The ice and snow species are neglected in these simulations because only warm-rain columns are considered.

Contour frequency altitude diagrams (CFADs) of the radar observables simulated for the warm-rain profiles are shown in Fig. 2. The freezing level is located at about 5 km. Two distinct hydrometeor populations can be seen in the reflectivity and in the Doppler velocity data. In the dominant mode, the cloud top height is about 1 km above the freezing level, where the Doppler measurements do not exceed 1.5 m s⁻¹. This corresponds to raindrop diameters less than 0.3 mm (Fig. A1c) that are characteristic of drizzle and cloud droplets. The velocity tends to increase towards the ground, indicating an increase in the size of the raindrops caused either by collision-coalescence processes or growth by condensation. The reflectivity profiles reach their maximum at approximately 4 km, and then they tend to decrease toward the ground, which may be due to the signal attenuation, a decrease in the water mass content, non-Rayleigh scattering effects (Kollias et al., 2002) or a combination of some of these factors. The secondary mode of the radar observables corresponds to more shallow precipitation columns, with the cloud top height between 2 and 4 km above the ground. This suggests the presence of a liquid cloud at this altitude too. Although similar peak reflectivity values are observed, the Doppler velocity is reduced compared to the deeper profiles, indicating smaller rain drops with a higher concentration and thus completely different microphysics. The presented simulations cover precipitation rates up to 15 mm h⁻¹, with a mean value of 0.4 mm h⁻¹. This is reflected in the PIA values, normalized by the cloud top height, shown in Fig. 2c.

Figure 2 Contour frequency altitude diagrams of the (a) radar reflectivity (Z_W) and (b) mean Doppler velocity (MDV_W) of the warm-rain profiles used for the C-CLD validation. Panel (c) shows a histogram of the PIA_SRT estimates (one-way) normalized by the cloud top height (CTH).

Validation of the retrieval was performed using approximately 8000 warm-rain columns, and its performance is illustrated in Fig. 3. The algorithm accuracy and precision are quantified by the mean error (ME) and root mean square error (RMSE) in the retrieved variables. The correlation coefficient (r) and normalized RMSE $\left(\text{NRMSE}\right(x)=\text{RMSE}(x)/\text{SD}(x)$, where SD(x) denotes the standard deviation of x) are computed as additional quality metrics. Since the considered variables are given in the logarithmic units, i.e. log ₁₀CLWP, log ₁₀MC and log ₁₀D_m, the ME and RMSE are given in bels (B). On average, the algorithm is overestimating the liquid cloud water path by about 32 % ($\text{ME}=-\mathrm{0.12}$ B and 10^0.12≈1.32). For profiles with higher cloud water content, the overestimation is reduced but scattered more around the 1–1 line. An opposite behaviour is observed for rain; the algorithm underestimates the rain MC by approximately 26 % (ME=0.13 B and ${\mathrm{10}}^{-\mathrm{0.13}}\approx \mathrm{0.74}$) to compensate for the PIA overestimation due to the cloud droplets. The retrieval of D_m shows very good accuracy; for $\mathrm{0.1}\le D_{\mathrm{m}}<\mathrm{3}$ mm, the algorithm tends to underestimate the characteristic size only by 5 % ($\text{ME}=-\mathrm{0.02}$ B). Because the same forward model was used for the retrieval and the scene simulations, the systematic underestimation for large sizes is believed to be caused by non-uniform beam-filling (NUBF) effects; i.e. the antenna-pattern-averaged mean Doppler velocity is smaller than the Doppler velocity corresponding to the footprint-averaged D_m because of the shadowing effect due to attenuation in correspondence to the fraction of the footprint with larger reflectivities (see Fig. 9 in Mroz et al., 2018). The precision of the algorithm is greatly reduced when PIA estimates are not assimilated in the retrieval, which is reflected in a reduction in the correlation coefficient and an increase in the RMSE values, as can be seen in Fig. 3d and e. The RMSE value increases from 0.24 and 0.27 to 0.48 and 0.45 B, while the correlation drops from 0.81 and 0.94 to 0.31 and 0.76 for CWP and rain MC, respectively. The estimate of the characteristic size is not affected by the lack of PIA measurements because it is mainly retrieved from the Doppler velocity measurements.

Figure 3 The warm-rain algorithm performance histograms. The x axis represents the model values, while the y axis corresponds to the retrieval. Panels (a), (b) and (c) show the cloud liquid water path, rain water content and rain characteristic diameter, respectively. Panels (d), (e) and (f) show the cloud liquid water path, rain water content and forward path-integrated attenuation when assuming that the PIA is not available. The reported values of the ME, RMSE, normalized RMSE (NRMSE) and correlation coefficient (r) are calculated for unknowns in the logarithmic units, i.e. log ₁₀CLWP, log ₁₀MC and log ₁₀D_m.

When PIA measurements are available, the simulated PIA is practically the same as the one being assimilated, with small differences due to the assumed error in the PIA measurements (i.e. 1 dB), giving a correlation of 0.99 and RMSE of 0.07 dB. When the PIA measurements are not available, the algorithm estimates the PIA using the maximum value of the reflectivity profile and the value close to the surface. While this approximation is useful, the lack of an integral constraint makes the correlation between measured and retrieved PIA drop to 0.58 and the RMSE increase to 0.39 dB as shown in Fig. 3f. When raindrops are present in the CPR radar sampling volume, they dominate the CPR observables. In this case, the information provided by the radar reflectivity and mean Doppler velocity is not sufficient to predict the PIA values reported by the model well. This results in a tendency to overestimate the amount of liquid cloud water content and thus to overestimate the observed attenuation.

The quality of the mass content retrieval can be further improved when the PIA estimate based on the surface reference technique (Eq. 32) is corrected for NUBF. We quantify this by replacing the PIA_SRT estimate with the fine-scale antenna-pattern-averaged attenuation values. In that case, the bias and RMSE in the rain MC estimate are reduced by 14 percentage points for both metrics. This indicates the need for more research on the NUBF and the related forward model adjustments, even in the case of satellite systems with such small footprints as the EarthCARE CPR (Battaglia et al., 2020a).

3.2 Cold rain and snow

The cold-rain and snow retrieval was applied to all the simulated scenes; in Fig. 4 the “Halifax” scene over eastern Canada is presented (for more detail on the simulated scenes, see Donovan et al., 2023). The left-hand side panels show the model output, while the right-hand side panels depict the retrieval and the simulated radar observables. The first part of the scene is occupied by light and moderate snow, with the cloud tops below 5 km. The second part presents ice clouds reaching 8 km and the associated stratiform precipitation with the melting layer between 2 and 3 km a.s.l., clearly highlighted by a sharp change in the Doppler signal. The cold-rain part features a heavy precipitation band associated with convection where the rain rates exceed 10 mm h⁻¹.

Figure 4 Panels (a), (b), (c) and (d) show the model radar reflectivity, Doppler velocity, mass content and mean-mass-weighted characteristic diameter for the ice and cold-rain regions of the Halifax scene. Panels (e), (f), (g) and (h) show the forward radar reflectivity and Doppler velocity and the retrieved mass content and mean-mass-weighted characteristic diameter applied to all regions, where $Z>-\mathrm{21}$ dBZ. The grey band represents the melting layer where the retrieval is not applied.

Overall, the C-CLD algorithm reliably reproduces the radar measurements corresponding to the precipitation structure, and despite being designed for stratiform rain, it performs relatively well even for convective profiles characterized by moderate precipitation conditions. The largest differences between the simulations and the retrieval within the stratiform rain systems are observed around the along-track distance of 3450 km close to the ground. The problem that affects these profiles is a misclassification in the C-TC product of pixels with a mixture of ice and rain as pure-snow columns. This leads to erroneously large D_m estimates and failure of the algorithm. The worst performance of the algorithm is observed for the retrieval of D_m in snow. Due to limited information content about the density of ice in the instrument illuminated volume, different degrees of riming tested by the algorithm can provide comparable cost function values. Therefore, the choice of the final solution may not be entirely accurate. Future work on the algorithm should focus on including such ambiguities in the final uncertainty estimates of the state vector.

Statistics on the retrieval accuracy based on all the scenes are presented in Fig. 5. The results for snow microphysical parameters combine the solid-phase part of cold-rain profiles and pure-snow profiles. These statistics correspond to radar reflectivity values in excess of −21 dBZ, where Doppler velocity is considered reliable, and the retrieval shows a full potential of Z and U_D measurements. The snow MC retrieval is strongly correlated with the model output, with a slight tendency to underestimate. The reported RMSE of 0.23 B corresponds to a fractional uncertainty of 53 %. The retrieval of D_m is more ambiguous (which is reflected in higher values of NRMSE) due to the limited variability in Doppler measurements with snow size, especially in the case of low-density ice (Fig. A2). Moreover, due to non-Rayleigh effects the reflectivity is not a monotonic function of the size, which additionally hampers the retrieval. This results in a moderate correlation coefficient of 0.68. As expected, the sizing retrieval in rain has a higher correlation and lower RMSE values than in ice due to the tighter relationship between the size and mean Doppler velocity. Like in the case of warm rain, the algorithm underestimates sizes above 0.6 mm and underestimates the MC values. However, in cold rain these differences are more pronounced because, in addition to NUBF effects, they are amplified by differences between the forward model used in retrieval and the one used in GEM simulations (not shown). Our forward model provides higher Doppler velocity for sizes exceeding 0.7 mm and smaller velocities below this size, which explains differences in the retrieved size. These differences propagate further into the MC retrieval. For D_m<0.7 mm the radar reflectivity increases with size, so an overestimate in D_m causes negative bias in the MC retrieval. When D_m>0.7 mm, reflectivity decreases, and thus the MC is underestimated also for large raindrops.

Figure 5 The algorithm performance histograms based on the three GEM scenes. The x axis represents the model values, while the y axis corresponds to the retrieval. Panels (a) and (b) show the rain water content and rain characteristic diameter. Panels (c) and (d) show the snow water content and snow characteristic diameter. The profiles with significant contributions of graupel and hail and the regions at cloud tops where the measurements are not very well constrained (large U_D error) are excluded from the analysis.

3.3 Stability and sensitivity of the optimal estimation biases in the measurements, the forward model and the a priori estimate

The calibration of radar systems and correct assumptions on microphysics are paramount for the accuracy of remote sensing retrievals. This is presented in Fig. 6, where the precision of the C-CLD algorithm in rain with various error sources is tested. The quality of the retrieval is quantified in terms of the NRMSE. First, the sensitivity of the retrieval to biases in the measurements is tested by adding a constant offset in the forward model to the radar reflectivity and the Doppler velocity. Note that this is equivalent to adding a bias with an opposite sign to the measurements; thus the calibration errors and model biases are tested simultaneously. As expected, the retrieval of the MC is mainly affected by the biases in the reflectivity, which is manifested in the valley-like shape of a local minimum, with the RMSE changing mainly along log ₁₀MC direction. Having said that, some compensation effects are also observed; i.e. the RMSE shows little variability if Z and U_D are simultaneously increased or decreased according to the direction given by the valley shape as shown in Fig. 6a. This is due to the characteristics of the forward model, namely the fact that the reflectivity depends on both the size and the mass content of rain. Thus, for a fixed mass, deviations in the reflectivity are compensated for by changes in D_m corresponding to changes in the Doppler velocity. The accuracy of the retrieval of log ₁₀D_m is also driven mainly by one variable, the mean Doppler velocity. As for the other unknowns, biases in the Doppler velocity measurements can be, at least to a certain degree, compensated by an offset in the reflectivity. However, due to a very constrained relation between D_m and U_D, the compensation is not as effective as for the MC retrieval and is mainly driven by the Z–D_m relationship used for the a priori estimate.

Figure 6 The accuracy of the C-CLD algorithm applied to the Halifax scene quantified by the NRMSE of D_m and MC. This evaluation takes into account the impact of (a, b, c) radar simulator bias and (d, e, f) a priori assumption bias in MC and D_m estimates. Panels (c) and (f) display the average NRMSE values for MC and D_m. The cross symbolizes the location of the minimum NRMSE.

The position of a local minimum of the average of the NRMSE in log ₁₀D_m and log ₁₀MC indicates the “reference” point that provides the best possible retrieval. As can be seen, the minimum is shifted away from the origin, which indicates differences between the forward model used in the retrieval and the one used for simulations. An offset of 0.4 m s⁻¹ and −1 dB in U_D and Z, respectively, would improve the reported retrieval uncertainties, but we decided not to alter our radar simulator as there is no evidence of the GEM model assumptions being superior to those used in the retrieval.

A similar analysis was performed to quantify the effect of the a priori assumption on the quality of the retrieval. As expected, the retrieval of D_m is mainly affected by its a priori estimate, and the same applies to the retrieval of MC. The error in the MC estimation resulting from differences between forward models (or calibration errors) can be reduced by changing the a priori assumptions, and the optimal retrieval is obtained if log ₁₀MC is increased by approximately 0.5 (i.e. a factor of 3.2). This indicates that for a given reflectivity value in rain, the mean mass content in the GEM model is larger than our a priori estimate. It raises the question of whether the Z–MC relationship, based on the PSD measurements at the ground (which typically fail in detecting small raindrops and low rain rates), which provides the basis of our forward model, is applicable in the low-precipitation-rate regime (Z<10 dBZ) that constitutes the majority of the profiles tested here.

3.4 The added value of the Doppler measurements

The EarthCARE radar mission is a follow-up of the highly successful CloudSat spaceborne radar mission. A number of studies on clouds and precipitation properties were conducted based on the CloudSat measurements (Stephens et al., 2018; Luo et al., 2008; Matrosov and Heymsfield, 2008; Tourville et al., 2015). The EarthCARE CPR is more sensitive (5–6 dB); has better vertical and horizontal sampling; has a smaller instantaneous field of view; and, most of all, has Doppler measurement capabilities. With all of these assets, it is vital to determine what improvement in the understanding of the properties of precipitation and clouds will be brought by the new mission.

The analysis presented here focuses on the Doppler velocity measurements value and their impact on the retrieval. The evaluation is based on comparison of the retrieval statistics with and without Doppler information assimilation. For this purpose, the C-CLD algorithm is applied once more to the Halifax scene, this time assuming no Doppler measurements. As expected, this results in reduced quality of the characteristic-size estimate. In rain, the correlation coefficient between D_m from the model and the retrieved one drops from 0.79 (for the original algorithm) to 0.47 in the no-Doppler setting. Similarly, the RMSE is approximately doubled for the reduced-input retrieval. This decreased confidence in the size estimate propagates to the retrieval of the rain water content, and it results in an RMSE increase from 0.29 to 0.44 B. The correlation drops from 0.92 to 0.79. Importantly, the lack of velocity measurements has no effect on the accuracy of the retrievals, with the mean error being almost non-affected.

The restriction of the measurement vector to radar reflectivity only has a small effect on the retrieval in the snow/ice. The RMSE and the ME statistics are virtually unchanged. Having said that, the correlation coefficient of the snow characteristic size decreases from 0.68 to 0.36 for no-Doppler retrieval. All the retrieved sizes oscillate in a narrow range of values between 0.2 and 0.3 mm that corresponds to the a priori estimates for the range of the observed reflectivities. This shows that Doppler measurements are relevant to estimating the size of snowflakes. Insignificant differences in the RMSE values between the size retrievals with and without velocity observations are due to the relative uncertainty in the velocity observations; i.e. the assumed measurement uncertainty of 0.2 m s⁻¹ gives a large fractional uncertainty in snow, where falling velocity often does not exceed 1 m s⁻¹.

One of the advantages of OE algorithms is the ability to quantify the amount of information provided by an individual measurement. This is achieved by comparing the state vector uncertainties before and after the measurements are assimilated by the algorithm (Shannon and Weaver, 1949). In geometrical terms, the information content of an observation is defined as the ratio between the volume enclosed by 1 standard deviation of the prior and posterior probability density function of X. For Gaussian distributions, this can be computed as follows:

$$\begin{array}{}\text{(33)}& H_{\mathrm{s}}=\mathrm{0.5}\ln \left|{\mathbf{R}}_x^{-\mathrm{1}}{\mathbf{R}}_{\mathrm{a}}\right|,\end{array}$$

where $\left|\right|$ denotes the determinant of a matrix and R_x and R_a are posterior and prior covariance matrices defined in Sect. 2.2.3 (see Eq. 2.73 in Rodgers, 2000). The computation of ${\mathbf{R}}_x^{-\mathrm{1}}$ (Eq. 18) for different instrument configurations does not require multiple runs of the computationally expensive algorithm. Instead, once the retrieval has converged, the diagonal elements of the matrix ${\mathbf{R}}_{\mathrm{a}}^{-\mathrm{1}}$ that correspond to selected measurements can be set to 0 to mimic instrument shutdown. This allows for the quantification of information content for all measurements together, for just radar reflectivity or just Doppler velocity or even for a single measurement at a given height in the column as shown in Fig. 7b and c.

To ensure a fair comparison among various regimes, the subsequent analysis focuses solely on quantifying the information content of the measurement in relation to estimates of mass content and characteristic size. Factors such as a reduction in uncertainty in melting layer attenuation or liquid cloud water content are not taken into consideration, as these variables are not present in all OE retrievals.

The amount of information provided by EarthCARE CPR measurements varies depending on the size and type of the hydrometeor being observed. In general, radar reflectivity provides more information for ice and snow, while mean Doppler velocity is more informative for rain. This trend is particularly noticeable in cold-rain columns, where the information content of reflectivity decreases from 3.2 to 2.5 nat (natural unit of information) as the hydrometeor transitions from a solid to liquid phase. In contrast, the information content of Doppler velocity increases from 2 to 2.5 nat during the same transition. The Doppler velocity measurements are useful for decreasing uncertainty in precipitation size estimation, particularly in rain, where the information content can surpass 2.5 nat. In snow, the observed sedimentation velocities have a smaller dynamic range, which results in reduced information content. High information content of radar reflectivity in ice can be attributed to an effective reduction in the uncertainty in the ice water content.

The amount of information is not uniform for a given measurement or hydrometeor type, as it depends also on the precipitation size. This non-uniformity in the information content is most apparent in warm-rain profiles, but it is also evident if rain pixels from Halifax simulations are compared with the “Cabo Verde” scene, where the retrieved sizes tend to be larger. In the case of reflectivity in warm rain, the information content is the highest at the top of the precipitation column, where particles are smaller than approximately 0.8 mm. Then, it reaches a minimum at the raindrop size that is the most efficient for the backscattering radar reflectivity signal (see the radar reflectivity maximum in Fig. A1a). For sizes larger than 0.8–0.9 mm, the information content is approximately 2 nat. A similar behaviour is observed for Doppler velocity but with a less pronounced minimum at 1 mm, where a reduction in the slope of U_D is observed. The maximum in the information content is observed for D_m=0.4 mm.

The total amount of information available from EarthCARE CPR measurements ranges from 3 to 4.4 nat, depending on factors such as hydrometeor type, ice-to-rain layer thickness ratio and particle size. Upon analysis of individual components, it is evident that, in snow and cold rain, the radar reflectivity profile (with PIA) provides the most information out of all measurements considered. It is followed by the Doppler velocity profile, as demonstrated in Fig. 7d. In contrast, in warm rain, both reflectivity and velocity measurements carry a similar amount of information. This trend also applies when considering only the liquid-phase precipitation in cold-rain profiles, as can be seen in the lower portion of Fig. 7b and c.

Figure 7 Information content for different measurements. Panel (a) shows the radar reflectivity for the context. Panel (b) displays the information content for individual measurements of radar reflectivity in the column. Similarly, panel (c) shows information content for individual measurements of Doppler velocity, and panel (d) shows the information content of an entire profile of measurements, as indicated in the legend. To ensure consistency across the different profiles, the information content values have been normalized by the number of retrieval levels.

The analysis presented here shows that the Doppler measurements are particularly valuable in remote sensing retrievals because they offer additional information that complements the measurements of reflectivity and PIA. In fact, the information content from all measurements combined ($Z+U_D+\text{PIA}$) is typically about 1 nat greater than the information content of reflectivity and PIA alone, demonstrating the significance of Doppler observations. That said, it is important to note that the information content of Z and PIA summed with the information content of U_D is larger than the information content of all measurements. Thus, there is some overlap between the information content of Doppler velocity and reflectivity, and they are not entirely independent. Furthermore, it is interesting to note that the information content of Doppler velocity is comparable for all considered regimes, while the reflectivity measurements are more advantageous in ice and cold rain.

4 Conclusions

The cloud and precipitation microphysics (C-CLD) algorithm is an EarthCARE Level 2 (L2) data product that utilizes measurements from the EarthCARE 94 GHz Doppler Cloud Profiling Radar (CPR) to provide microphysical information about cloud and precipitation systems. The C-CLD algorithm primarily uses an optimal estimation (OE) approach that balances the information provided by the CPR with a priori knowledge on the climatology of cloud and precipitation systems. The algorithm is designed to retrieve profiles of two moments of the PSD drop size, namely the mass content and mean-mass-weighted diameter, owing to the information content provided by the CPR. When dealing with drizzle-free and lightly drizzling warm clouds, the OE framework is replaced with climatological relationships between the measured reflectivities and the microphysical parameters of interest.

A large dataset of in situ, surface-based observations is used to reduce the number of free parameters and to obtain the forward model relations with the corresponding uncertainties. To maintain water mass flux through the melting layer, only small perturbations from this condition are allowed. Additionally, a one-dimensional parameterization is proposed for representing a wide range of ice particle densities, from unrimed snowflakes to dense graupel particles.

The C-CLD retrieval framework has been applied to EarthCARE CPR simulated observations from high-resolution weather systems simulations occurring in three different climatological regimes (Donovan et al., 2023): tropical climate, humid continental climate bordering on an oceanic climate (Halifax) and mid-latitude conditions over North America (Baja). The CPR reflectivity and Doppler radar measurements provide sufficient information to retrieve, with high confidence, two moments of the PSD, especially in rain due to the added value of the Doppler measurements, which, in stratiform rain, are closely related to the raindrop fall speed and thus to its mean size. On average, the mean-mass-weighted diameter (D_m) of rain can be estimated within a precision of 23 % with negative bias reported for large sizes. As a result, the estimate of rain mass content (MC) is also captured well by matching the radar reflectivity to the observations. The uncertainty in the MC estimate is estimated to be 67 % for all the GEM simulation scenes combined. Despite more complex and ambiguous scattering properties of ice particles, errors in the ice mass content are smaller than in rain, and they are equal to 53 % for profiles including either snow-only or cold rain. This unexpected result may indicate differences in the forward model used in GEM simulations and in the retrieval of rain or difficulties in separating path-integrated attenuation (PIA) into the liquid cloud, melting layer and rain components. The retrieval of D_m of ice is the most challenging; it is characterized by the lowest correlation coefficient and the highest value of the normalized root mean square error among all the considered unknowns. The variety of snowflake morphology and the corresponding diversity in the relation between particle size and terminal velocity results in uncertainties of 23 %.

Due to the high susceptibility of W-band measurements to signal attenuation, the quality of the retrieval is strongly reduced when the path-integrated attenuation estimates are not assimilated. This is reflected in the degradation of the quality of the mass content retrieval in warm-rain conditions; i.e. the RMSE in log ₁₀MC increases from 0.27 to 0.45 B, while the correlation coefficient is reduced from 0.97 to 0.75.

The reported uncertainties are heavily dependent on the forward model accuracy and on the measurement calibration biases. The performed analysis revealed that, due to some differences in the fall velocities used in the GEM model and in the C-CLD retrieval framework, a systematic overestimation (underestimation) of small (large) raindrop sizes is present. These errors, combined with discrepancies in the reflectivity forward model, result in a negative bias of the rain mass content. The differences between the simulators are attributed to the different particle size distribution shape assumptions. Although the difference between the radar simulators was not systematic (i.e. it has a different sign depending on the rain characteristic size), the bias in the mass content retrieval was. This shows how susceptible to model and/or measurement biases the optimal estimation framework is and how important the calibration of the EarthCARE reflectivity and Doppler velocity (Battaglia and Kollias, 2014) will be.

Despite the detection of differences between the GEM simulations and our radar model, the algorithm was not fine-tuned to match model assumptions due to the lack of evidence that the model could reflect reality better than the long-term particle size distribution statistics. This paper aims at providing a physical basis for the retrieval, and so the modifications of the forward model or of a priori assumptions are left for the calibration–validation activity period after the launch of the satellite.

Thanks to its large antenna, CPR's unprecedented fine horizontal resolution minimizes the impact of two of the challenges of spaceborne radar-based precipitation remote sensing: multiple scattering (Battaglia et al., 2010; Matrosov et al., 2008; Matrosov and Battaglia, 2009) and non-uniform beam filling (NUBF; Tanelli et al., 2012). Since the horizontal resolution of the model simulations is finer than the one of the radar, errors related to NUBF are quantified and included in the reported total algorithm errors, with small biases observed. Furthermore, negligible multiple scattering effects were simulated and reported to the flag produced in C-PRO (Kollias et al., 2022b).

The CloudSat mission radar measurements and the information about clouds and precipitation derived from them provide a strong heritage for the C-CLD product development. However, this legacy poses a challenge for consistency of the retrieved parameters for the two missions. Continuity between the products is important for evaluation of the long-term trends in precipitation statistics and climate change studies. To address this issue, the information content of Doppler measurements and their impact on the retrieval must be carefully evaluated. The knowledge about the cloud and precipitation properties gained with this additional measurement can be transferred to the reflectivity-only algorithm via refinement of the a priori assumptions. The updated C-CLD algorithm (i.e. without Doppler measurements and considering differences in the instrument specifications) can be applied to the CloudSat measurements and compared to its products for an assessment of potential biases. In the case of detection of systematic differences, the CloudSat dataset can be reprocessed to provide a consistent long-term cloud and precipitation record.

Our preliminary analysis shows that the amount of information provided by EarthCARE CPR measurements varies depending on the size and type of hydrometeor observed, with reflectivity more informative for ice and snow and mean Doppler velocity for rain. The Doppler velocity measurements are particularly useful in reducing uncertainty in precipitation size estimation, especially for rain. The non-uniformity in the information content is most apparent in warm-rain profiles, where the size of particles is evolving with height. The maximum in the information content of U_D is observed for D_m=0.4 mm. The analysis reveals that the Doppler measurements complement the measurements of reflectivity and PIA, and the information content of all measurements combined is typically about 1 nat greater than the information content of reflectivity and PIA alone.

Further development of the algorithm is necessary to ensure its effectiveness in a wider range of weather conditions. The GEM simulations used in this study do not include weather systems with raindrops or snowflakes larger than 1 mm. These conditions are particularly challenging for W-band retrievals due to significant signal attenuation and saturation of radar reflectivity and Doppler measurements (Mróz et al., 2019). To improve the credibility of a priori estimates in the drizzle size particle regime, more in situ measurements are required. This aspect should be addressed during the calibration–validation activities. Moreover, characterizing the shape of the liquid cloud mass content profiles is essential for reducing uncertainties in path-integrated attenuation simulations and retrieved rain and snow mass content below the liquid cloud top. As suggested by Battaglia and Panegrossi (2020), this issue can be mitigated by the inclusion of the W-band brightness temperatures in the observables adopted in the derivation of the C-FMR product. In addition, in order to produce realistic transitions in the retrieved state vector between consecutive profiles, future algorithms could make use of the two-dimensional information provided by the radar. In addition, future algorithms could leverage the two-dimensional information from the radar to create seamless transitions between consecutive profiles. This not only preserves the horizontal continuity of the state vector but also enables accurate quantification and correction of NUBF effects. Additionally, using a two-dimensional approach can aid in detecting liquid clouds that may be present within precipitation. These clouds typically have long correlation lengths and are easier to spot outside of precipitation zones where larger hydrometeors are not masking them. Once these clouds are detected in these areas, their location within precipitation can be inferred through continuity. Finally, future work should include cross-validation with the other precipitation products, e.g. ACM-CAP (Mason et al., 2023), which provides a synergistic retrieval of the hydrometeor properties based on the full suite of sensors on board the EarthCARE satellite. This latter product should provide more accurate estimates due to the increased information content provided by the other instruments.

Appendix A: Parameterization of scattering properties in the W-band

Here, we report the parameterizations of the scattering properties at 94 GHz that are used in C-CLD. These relations link the CPR observables (reflectivity Z_e and Doppler velocity U_D) with two state vector parameters (D_m and MC) in terms of power laws. This simplifies the analytical computation of the Jacobian.

A1 Rain

The radar observables and PSD moments are approximated by polynomials in $x={\log }_{\mathrm{10}}\text{MC}\left[\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}\right]$ and y=log ₁₀D_m[cm], i.e.

$$\begin{array}{}\text{(A1)}& {\displaystyle }Z\left[\mathrm{dBZ}\right](x,y)=\mathrm{10}x+P_Z\left(y\right)=\mathrm{10}x+\sum _{i=\mathrm{0}}^n{c}_i^Z{y}^i,\text{(A2)}& {\displaystyle }{U}_D\left[\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}\right](x,y)=P_{U_D}\left(y\right)=\sum _{i=\mathrm{0}}^n{c}_i^{U_D}{y}^i,\text{(A3)}& {\displaystyle }\begin{array}{rl}\mathrm{10}{\log }_{\mathrm{10}}\left(k\right[\mathrm{dB}{\mathrm{km}}^{-\mathrm{1}}\left]\right)(x,y)& =\mathrm{10}x+P_k\left(y\right)\\ & =\mathrm{10}x+\sum _{i=\mathrm{0}}^n{c}_i^k{y}^i,\end{array}\text{(A4)}& {\displaystyle }\begin{array}{rl}\mathrm{10}{\log }_{\mathrm{10}}\left(\text{PR}\right[\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}}\left]\right)(x,y)& =\mathrm{10}x+P_{\mathrm{PR}}\left(y\right)\\ & =\mathrm{10}x+\sum _{i=\mathrm{0}}^n{c}_i^{\mathrm{PR}}{y}^i,\end{array}\end{array}$$

where the coefficients $c_i^f$ for $f=Z,k,U_D$ and PR are given in Table A1. The degree of the fitting polynomial results from its high accuracy in replicating the simulations for the gamma PSD model over a broad range of characteristic rain sizes, i.e. from 0.1 to 3.5 mm.

Figure A1 Two-dimensional histograms of the radar-observable simulations corresponding to the in situ PSD measurements at the ground. (a) Radar reflectivity factor in dBZ per 1 g m⁻³ of rain. (b) 10log ₁₀ of the (one-way) specific attenuation in dB km⁻¹ per 1 g m⁻³ of rain. (c) Mean Doppler velocity in standard atmospheres. (d) Precipitation rate in standard atmospheric conditions (15 ^∘C, 1013.25 mbar) per 1 g m⁻³ of rain. The black line shows the simulations for the gamma PSD model with $\mathit{\mu }=\mathrm{10}{D}_{\mathrm{m}}^{-\mathrm{0.8}}-\mathrm{4}$ that is used as a forward model.

Table A1 Coefficients of the polynomial representation of the forward model for rain given by formulas (A1), (A2), (A3) and (A4).

A2 Cloud attenuation coefficients

The two-way attenuation coefficient in dB km⁻¹ g⁻¹ m³ is parameterized as a quadratic function of the temperature expressed in Celsius with the zeroth-, the first- and the second-order coefficients equal to $[\mathrm{8.4979},-\mathrm{0.0062},-\mathrm{0.0022}]$. This replicates the empirically verified model at 94 GHz very well (Tridon et al., 2020, Fig. 1) with a maximum of about 8.5 dB km⁻¹ g⁻¹ m³ at 271.8 K decreasing to 7 dB km⁻¹ g⁻¹ m³ at 245.7 and at 297.8 K.

A3 Ice

As in the case of rain, the radar observables and PSD moments are approximated by polynomials in $x={\log }_{\mathrm{10}}\text{MC}\left[\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}\right]$ and y=log ₁₀D_m [cm]. These polynomials are of different degrees, and their coefficients depend on the degree of riming. Therefore, it is impractical to list all the coefficients here. Instead, these tables are freely available at https://doi.org/10.5281/zenodo.7529739 (Mroz, 2023). Figure A2 shows the scattering properties as parameterized in the forward model for five selected degrees of riming. Note that while attenuation and Doppler velocities tend to increase with melted diameter, reflectivities reach maximum values in correspondence to sizes between 0.4 and 0.8 mm.

Figure A2 Ice-scattering properties as parameterized in the forward model as a function of the equivalent melted size; different colours correspond to different degrees of riming (α). For selected values of α, the histograms in the background show the gamma PSD modelling corresponding to the rain PSD measurements collected at the ground based on the “melting-only” assumption (see Sect. 2.4.1). The line represents $\mathit{\mu }=\mathrm{10}{D}_{\mathrm{m}}^{-\mathrm{0.8}}-\mathrm{4}$. All the simulations are performed for 1 g m⁻³ of snow. (a) Radar reflectivity factor is in decibels relative to Z. (b) 10log ₁₀ of the (one-way) specific attenuation is in decibels per kilometre. (c) Mean Doppler velocity is in standard atmospheres. (d) Precipitation rate is in standard atmospheres.

Appendix B: Other parameterizations

The vertical gradient of the change in adiabatic LWC is calculated as in Rogers and Yau (1989) and is parameterized as

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}{\displaystyle \frac{\partial \text{LWC}}{\partial z}}& ={\mathit{\rho }}_{\mathrm{air}}{\displaystyle \frac{c_p}{L_{\mathrm{ev}}}}\left({\mathrm{\Gamma }}_d-{\mathrm{\Gamma }}_s\right)\\ & \approx \left(d_{\mathrm{0}}+d_{\mathrm{1}}{T}_c+d_{\mathrm{2}}{T}_c^{\mathrm{2}}\right){\left({\displaystyle \frac{p}{p_{\mathrm{90}\mathrm{kPa}}}}\right)}^{c_{\mathrm{0}}+c_{\mathrm{1}}{T}_c},\end{array}\end{array}$$

with $[d_{\mathrm{0}},d_{\mathrm{1}},d_{\mathrm{2}},c_{\mathrm{0}},c_{\mathrm{1}}]$ = $[\mathrm{1.615},\mathrm{0.0554},-\mathrm{6.287}\times {\mathrm{10}}^{-\mathrm{4}},\mathrm{0.42},\mathrm{0.015}]$.