Vertical wind and drop size distribution retrieval with the CloudCube G-band Doppler radar

Macrophysical properties of clouds are influenced by underlying microphysical processes. In practice, there is often an observational gap in bridging the two. For example, our current understanding of aerosol–cloud interaction and cloud–climate feedback is hindered by a lack of robust measurements of the distribution of drop sizes within clouds, especially for the smallest drop sizes. Doppler radar measurements have proven useful in estimating rainfall drop size distributions (DSDs) but face an intermediate challenge of requiring a correction for the presence of vertical air motion. Recent advances in millimeter-wave technology have made radar measurements at increasingly smaller wavelengths possible, allowing for analysis of particle-size-dependent scattering effects to derive estimates of vertical winds and thereby DSDs. This work demonstrates a method of deriving range-resolved DSDs using Doppler spectra at 238 GHz measured by the CloudCube ground-based G-band atmospheric Doppler radar. The observations utilized are of marine boundary layer clouds during March and April 2023 in La Jolla, CA, USA, taken as part of CloudCube's participation in the Eastern Pacific Cloud Aerosol Precipitation Experiment (EPCAPE) campaign. This method first identifies notches in the velocity spectra and compares them to the theoretical notch velocities predicted by size-dependent backscattering and terminal velocity models to estimate the range-dependent vertical wind. After removing the vertical wind, binned DSDs are retrieved from the zero-wind spectrum. Bulk properties of the precipitation are then derived, including the number concentration, liquid water content, characteristic drop size, and precipitation rate. For the case study presented here, calculated bulk properties are found to be relatively invariant to the forward-model assumptions made in the estimation of the full DSD retrieval. Validation of this method on larger volumes of data would make such retrievals useful tools in assessing physical models of drizzle.

1 Introduction

Marine boundary layer clouds represent the largest physical source of uncertainty in projections of climate sensitivity (Zelinka et al., 2020) and are central to understanding the radiative forcing of aerosol–cloud interactions (Bellouin et al., 2020). A consistent finding is the relationship between the occurrence of precipitation and the mesoscale organization of low clouds (Abel et al., 2017; Yamaguchi et al., 2017; Smalley et al., 2022), where a transition from closed-cell clouds to open-cell clouds is associated with precipitation onset. Therefore, these cloud transitions are critical in constraining both aerosol–cloud interaction and cloud–climate feedbacks. A current dilemma in climate projection is the fact that the accuracy of future projections is limited by a negative correlation between aerosol–cloud interactions and cloud–climate feedback (Gettelman et al., 2024). This anti-correlation has been clearly linked to climate model representation of the precipitation formation process (Suzuki et al., 2013).

Accurate measurements of drizzle and light rain are essential to improve the understanding of the microphysical processes in boundary layer clouds and constrain both the aerosol–cloud interactions and the cloud–climate feedback. Understanding how the size distribution of drizzle drops evolves in both time and space can provide insight into coalescence, breakup, and evaporation processes that shape cloud macrophysical properties. Current methods of directly measuring DSDs include using either ground-based or airborne disdrometers, devices which directly measure drop sizes. Ground-based disdrometers only measure drops which fall all the way to the surface and have a limited capability to observe the smallest precipitation drops (Wang and Bartholomew, 2023). Airborne measurements are more likely to capture data at several elevations; however these measurements are sparse.

Continuous observations from remote sensing measurements are necessary to fill the gaps in in situ sampling. Radar is the optimal tool for remote observations of drizzle and light rain from ground-based or airborne platforms. The most straightforward method to derive drizzle parameters is by assuming a reflectivity drizzle–rate (Z–R) relationship (Comstock et al., 2004). However, in practice, Z–R relationships have primarily been used operationally for satellite cloud radar observations where reliable Doppler observations are not available (Lebsock and L'Ecuyer, 2011; Mroz et al., 2023). The widespread proliferation of millimeter-wave Doppler cloud radars has enabled a new class of retrieval of drizzle and light rain that combines radar reflectivity and higher Doppler moments. For example, Frisch et al. (1995) combine Ka-band Doppler spectral moments with an assumption of zero vertical wind and an assumed drop size distribution shape to derive the vertical profile of drizzle parameters. O'Connor et al. (2005) built on this approach by combining W-band Doppler moments with lidar backscatter and a method to correct for turbulent broadening. This multi-sensor approach has subsequently been used to make novel observations of drizzle in stratocumulus clouds (Ghate and Cadeddu, 2019). Galloway et al. (1999) invert an airborne W-band Doppler spectrum to derive a binned drizzle DSD without assuming a DSD shape while retaining the zero mean wind assumption.

One common shortfall of the Doppler-based methods mentioned above is the difficulty in accounting for the effect of the vertical air motion on the mean Doppler. In this respect, Mie scattering in millimeter-wave radars can be useful in constraining the vertical air motion when the size of drops is similar to the observing wavelength of a radar system. Specifically, the backscattering efficiency at a particular observing wavelength contains several peaks and valleys as a function of drop radius, as seen in Fig. 1. Lhermitte (1987) first proposed the technique of using full W-band Doppler spectra that show similar oscillatory shapes along with information about the theoretical backscattering efficiency in each velocity bin to simultaneously retrieve information about vertical air velocity and rain DSD. The difference between the theoretical and observed locations of any minima seen in the Doppler spectrum would yield information about vertical air motion in the scene, while the relative heights of any maxima seen in the spectrum would yield information about the DSD of the scattering particles. Kollias et al. (2002) demonstrated vertical air motion retrievals derived from W-band Doppler spectrum observations of stratiform precipitation. Giangrande et al. (2010) was the first to demonstrate the full utility of this technique in analyzing W-band Doppler spectra. This work retrieved measurements of both vertical winds and best-fit parameters to a Marshall–Palmer log-linear DSD; however, useful spectra can only be captured at the W-band for storms with large droplets greater than ∼0.8 mm. As seen in Fig. 1, the first backscattering minimum occurs at a larger drop radius for lower observing frequencies. Thus, to be sensitive to measurements of drops as small as drizzle (smaller than a diameter of 0.5 mm), it is necessary to make Doppler spectrum measurements at a higher frequency. For the G-band in particular, the first minimum is located at a small enough radius that its location is insensitive to the parameterization of the drop aspect ratio, permitting high-accuracy quantification of the vertical air motion in all but the lightest liquid-phase precipitation. The utility of G-band observations was theorized by Battaglia et al. (2014) and first demonstrated by Courtier et al. (2024) in the addition of G-band Doppler spectra to a multi-frequency DSD retrieval.

Figure 1 Backscattering efficiency as a function of particle radius for three different radar bands. The solid lines represent backscattering calculated with the T-matrix method, which assumes some oblateness of the drops. The dotted lines represent Mie backscattering, which assumes the drops are spherical. The horizontal dotted lines represent lines of efficiency equal to 0 and 1 for each band. Calculations are described in Sect. 2.3.

This work explores the capability of making retrievals based only on G-band spectra to profile liquid-phase precipitation DSDs in marine boundary layer clouds. This class of precipitation is ideally suited for the G-band for two reasons: (1) the preponderance of small drops means that the Mie resonance (or notch) will frequently not be observed in W-band spectra but will be observed at the G-band, and (2) the liquid water content is small, and thus the hydrometeor attenuation tends to be small. The latter fact means that attenuation correction can be performed without incurring the large errors that are common in rainfall retrievals at frequencies used to observe more highly attenuated conditions (e.g., Ka-band or W-band; Hitschfeld and Bordan, 1954). To demonstrate these capabilities, this paper uses the first operational data from the deployment of a newly developed G-band Doppler cloud radar to a large field deployment at a coastal site with frequent marine boundary layer clouds and validates the results against ancillary observations.

2 Instrument and data overview

2.1 CloudCube instrument

CloudCube is a modular triple-frequency (Ka-band, W-band, and G-band) atmospheric radar instrument developed at the Jet Propulsion Laboratory. Its use of both a fully solid-state design and direct up/down conversion between the baseband and RF allows it to have a uniquely compact architecture, ideal for deployment in the field. This paper focuses specifically on the G-band channel, currently the only CloudCube channel with full Doppler spectral resolution. The observing frequency, 238.8 GHz, was strategically chosen to take advantage of an intersection between an allowed frequency allocation and a trough in the atmospheric absorption curve. It also lies close to the limits at which transmit sources of sufficient power are available. A summary of the instrument parameters are shown in Table 1, and more detail on the instrument can be found in Socuellamos et al. (2024a).

Table 1 Summary of KaZR and CloudCube G-band radar specifications at EPCAPE. FMCW – frequency-modulated continuous wave.

2.2 EPCAPE campaign

The data presented in this paper were collected as part of CloudCube's participation in the Eastern Pacific Cloud Aerosol Precipitation Experiment (EPCAPE), a campaign of the U.S. Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) user facility (Russell et al., 2021). The main goal of EPCAPE was to better understand marine stratocumulus clouds and their effect on Earth's radiation budget. CloudCube measured cloudy and lightly raining cumulus and drizzling stratocumulus over several days during March and April 2023 from atop Scripps Pier in La Jolla, CA, stationed adjacent to the ARM Mobile Facility (AMF). Data from all three bands of CloudCube were saved during this time. Details on the post-processing of the data can be found in Socuellamos et al. (2024c), and the datasets are made publicly available in Socuellamos et al. (2023). The majority of the precipitation events during this deployment period were not observed by CloudCube because at that time the instruments did not have radomes and had to be covered during periods of surface precipitation to protect the radars.

This paper also uses data taken from several ARM instruments to both aid and supplement the presented analysis. Notably, our retrievals rely on temperature, pressure, and humidity profiles collected by radiosondes (Holdridge, 2020) for determining the correct values of particle backscattering efficiencies, fall speeds, and gaseous attenuation. The retrievals are validated with the ARM Ka-band radar (KaZR; Widener et al., 2012; see the radar parameters in Table 1) and the 2-D video disdrometer (VDIS; Bartholomew, 2020) instruments.

2.3 Scattering properties

The single-scattering properties of liquid precipitation drops are calculated with the T-matrix method (Mishchenko and Travis, 1998), using the Python wrapper of Leinonen (2014). The aspect ratio of drops is modeled using the equation $\frac{b}{a}=\mathrm{1.055}-\mathrm{0.0653}$ D, where D is the drop diameter (mm), valid in the range 1.5–8 mm (Thurai and Bringi, 2005), and $b/a$ is the axis ratio of the spheroids. The aspect ratio is equal to 1 for the smallest drops (smaller than a diameter of 0.84 mm), for which this formula produces aspect ratios larger than unity. A lookup table is created with the drop single-scattering properties in 1 µm increments in radius and 1 K increments in temperature. The temperature-dependent refractive index is taken from Warren (1984). At the CloudCube observing frequency (238 GHz), the first minimum is located at a drop radius of 0.33 mm.

3 Data filtering and minima finding

CloudCube collected around 51 h of data over 13 separate days during its deployment. However, only a few instances spread over 2 d of this dataset contained spectra resolving at least one backscattering minimum for at least 0.5 continuous kilometers and 100 continuous seconds to perform robust retrievals. These values are chosen arbitrarily – future work may investigate selecting these thresholds adaptively. The spectrum in Fig. 2b, referred to hereafter as “the example spectrum”, is used to demonstrate our methods for the remainder of this paper.

Figure 2 (a) Reflectivity curtain for several minutes of G-band radar data collected on 30 March 2023. The pink line shows the location of the cloud base, as measured by the ARM laser ceilometer. (b) Mean Doppler velocity corresponding to the reflectivity curtain in (a). (c) A full Doppler spectrum from data collected during the time marked by the dotted black line in (a). The dashed red line shows the elevation for the example 1-D spectrum on the right. (d) The black line represents the original measured spectrum, whereas the pink line shows the Wiener smoothed spectrum. The thicker dashed vertical lines represent minima in the spectrum that are considered to be truly correlated with minima in the backscattering function, whereas the thin dashed lines represent other candidate minima that the find_peaks function found without enforcing any checks.

To provide context for this example spectrum, Fig. 2a shows the G-band reflectivity curtain and marks the time at which the example spectrum was collected. Also shown is the cloud base height as measured by the ARM infrared laser ceilometer (Morris, 2016). Figure 2c shows a characteristic Doppler spectrum of a precipitating postfrontal shallow cumulus ice-phase cloud. This cloud happens to have a cloud base near the melting level, which can be very clearly observed in Fig. 2b as the region near an altitude of 1.4 km, below which the narrow Doppler spectrum with small Doppler velocity rapidly broadens as large ice crystals melt into falling raindrops with significantly increased Doppler velocities. Note also how the melting layer appears as a weak bright band in panel (a) and a rapid change in the mean Doppler velocity in panel (b). In the liquid precipitation below the melting layer, there are multiple resonances in the Doppler spectra, corresponding to the Mie notches that can be exploited to find the vertical air motion. As the scattering properties of liquid drops and the resulting Doppler spectra are relatively straightforward to calculate compared to the complexity of modeling ice, this work limits itself to analyzing regions where only liquid water is present. An investigation of ice retrievals is saved for future work.

To begin, the spectra are smoothed in the range direction to filter out the smallest scale of vertical turbulence, improving the signal to noise in the minima detection. The spectra originally have a range resolution of 10 m, and a 1-D Gaussian blurring kernel is used to smooth to an effective range resolution of 50 m. Following the smoothing step, each 1-D velocity power spectrum is analyzed to identify any minima.

To discern true backscattering minima from Rayleigh-distributed noise fluctuations, we search for minima using a smoothed version of each 1-D spectrum. The spectrum is smoothed using a Wiener filter, and we use the difference between the smoothed and original spectra to estimate the standard deviation of the noise fluctuations. Next, we utilize the Python function scipy.find_peaks (Virtanen et al., 2020), which identifies a peak as any point where its two neighbors are of a lower value – using sign-inverted spectra allows the minima to show up as peaks. To separate true minima from spurious ones, a few criteria are imposed. We mandate that for true detections, the depth of the minima (measured as the peak to trough distance between the minimum and the nearest peak) must be at least 5 times the standard deviation of the noise fluctuations. Based on the spacing between successive minima seen in the T-matrix backscattering calculation, it is enforced that the spacing between minima found in the spectra must be separated by at least 1.15 m s⁻¹. For sets of minima found with smaller spacings, minima with lower signal to noise ratios (SNRs) are preferentially filtered out. Our final check is to ensure that we are not erroneously identifying points in the noise floor of the spectrum as minima by ensuring that for true detections there are continuous data points for at least 1 m s⁻¹ on each side of the minimum. A demonstration of this minima retrieval for a single elevation is shown in Fig. 2c.

After minima at all elevations have been identified in a Doppler spectrum profile from a single time, they are classified based on which minimum in the backscattering function they correlate to. The simplest method of doing this would be to draw boundaries of fixed width in velocity space and assume all points within each section are correlated to the same backscattering notch. However, this requires an a priori assumption of a mean vertical wind value with relatively low variance across elevation. To mitigate the risk of poor retrievals resulting from incorrect initial assumptions, a Gaussian mixture model (GMM) is used to cluster points associated with the same backscattering notch. This method assumes that all points in a given dataset are drawn from one of N multivariate Gaussian distributions, each with their own means and covariance matrices. We utilize the Python implementation of GMM in sklearn.mixture.GaussianMixture (Pedregosa et al., 2011). The algorithm initially looks for data points in the 1 m s⁻¹ regions around the locations of the first three theoretical backscattering minima for T=270 K. The number of components for the GMM to sort data into is determined by how many of these initial regions have data within them. For example, Fig. 3a shows a spectrum with these divisions overlaid. Since data are present in all three of the divisions, a GMM with three components is used. The initial locations of the three Gaussian components are decided using the mean elevation and velocity of the data within each division. The shape of each Gaussian is initialized by a covariance matrix derived from the standard deviations of the elevations and velocities of the data. An example of these initial guesses is also shown in Fig. 3a.

Figure 3 (a) A 2-D Doppler spectrum with all identified minima overlaid. The dashed black lines represent the boundaries used to define the initialization of the classifier. The blue stars represent the means of the points within the boundaries, and the colored ellipses represent the initial covariances in both the height and the velocity directions. (b) A 2-D Doppler spectrum with all identified minima classified according to their correlation to notches in the backscattering function. The cyan points are minima that were filtered out as being erroneous, and the ellipses show the probability spaces of the Gaussian mixture model components. (c) Vertical winds calculated from the classified minima.

The final components are fit by adjusting the parameters of the Gaussians until the likelihood of all points being drawn from one of the distributions is maximized. The final classifications are adjusted if necessary by ensuring that each height only contains one point from each distribution. Any outlier points that have anomalously low likelihoods are masked. An example of the final distributions can be seen in Fig. 3b. Note the cyan point to the far left of the spectrum that was identified as a minimum but is masked as an outlier due to being many standard deviations away from any of the three Gaussian components. Within each cluster, any remaining outlier points are identified using the LocalOutlierFactor (LOF) method in sklearn, which measures the local densities of points and identifies any points that have anomalously low local densities. In this example case, the LOF method identified the cyan point near 250 m, close to the purple points, as being an outlier point. However, it did not identify the two yellow points that seem to be, upon visual inspection, outside the distribution. Tuning this method exactly can be a challenge.

4 Vertical wind retrieval

After classifying the minima, the vertical winds are retrieved from their measured locations. The Gaussian component with the mean velocity closest to zero is assumed to contain the points that correspond to the first backscattering minimum. For each elevation, the temperature measured by the radiosonde is used to select the correct temperature-dependent backscattering efficiency function. As the backscattering efficiency is a function of diameter and our spectrum power is a function of velocity, drop diameters are transformed to drop velocity by assuming that all drops fall at terminal velocity. For drop diameters greater than 100 µm, we linearly interpolate between the data points presented in Gunn and Kinzer (1949) to calculate terminal velocity. For smaller drops, we use Stokes' law, $v_{\mathrm{t}}=\frac{\mathrm{1}}{\mathrm{4}}k{D}^{\mathrm{2}}$ ($k=\mathrm{1.19}\times {\mathrm{10}}^{\mathrm{8}}$ is a constant, D is in m, and v_t is terminal velocity in m s⁻¹). The effects of air density are corrected for by multiplying the terminal velocity by a correction factor, $C=({\mathit{\rho }}_{\mathrm{0}}/\mathit{\rho }{)}^m$, where ρ₀=1.204 kg m⁻³ (density for standard temperature and pressure) and $m=\mathrm{0.375}+(\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{5}})D$ (Beard, 1985). The air density, ρ, as a function of elevation is calculated using the temperature and pressure values measured by radiosondes. Once backscattering efficiency is transformed to be a function of velocity, the measured minimum value is subtracted from the theoretical value to retrieve vertical wind as a function of height: $v_{\text{wind}}\left(h\right)=v_{\text{meas}}\left(h\right)-v_{\text{theo}}\left(h\right)$. An example of this retrieval using each of the minima is shown in Fig. 3c. The colors of the vertical wind curves correspond to the colors of the minima used to derive the wind speeds.

There are small inconsistencies between wind speeds at the same height calculated from different minima. It is difficult to determine the exact cause of this inconsistency, but it is likely that a combination of uncertainty in the Gunn–Kinzer terminal velocity relationship and the drop obliquity parameterization are the largest contributors to this discrepancy. For a monotonically varying DSD, the details of the DSD should not affect the locations of the minima. For the rest of the analyses presented in this paper, only the vertical winds derived from the first minimum are considered (corresponding to the red points in Fig. 3). Recall that the location of the first minimum can be assumed to be insensitive to the drop obliquity parameterization.

5 Drop size distribution retrieval

As described in Kollias et al. (2011), the measured Doppler spectra can be described by the following equation:

$$\begin{array}{}\text{(1)}& S(v+v_{\text{wind}}{)}_{\text{obs}}=(A+{\mathit{ϵ}}_{\mathrm{a}}\left)\right[S(v{)}_{\mathrm{Q}}\cdot g({\mathit{\sigma }}_{\text{v,turb}}\left)\right]+{\mathit{ϵ}}_{\mathrm{s}},\end{array}$$

where v is the true particle velocity, v_wind is the vertical wind speed, A is attenuation, ϵ_a is the attenuation error, S(v)_Q is the quiet-air spectrum (no turbulence), g(σ_v,turb) is the convolution kernel that describes spectral broadening due to turbulence, and ϵ_s represents error in the measured spectral power. Retrieval of the quiet-air spectrum from the measured spectrum would enable measurement of the drop size distribution as a function of drop radius, N(r), in units of ${\mathrm{m}}^{-\mathrm{3}}{\mathrm{m}}^{-\mathrm{1}}$, using the following relationship:

$$\begin{array}{}\text{(2)}& S(v{)}_{\mathrm{Q}}={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{5}}\left|K\right(\mathit{\lambda }){|}^{\mathrm{2}}}}{\mathit{\sigma }}_{\text{bck}}(r\left)N\right(r){\displaystyle \frac{\mathrm{d}r}{\mathrm{d}{v}_{\mathrm{t}}}},\end{array}$$

where λ is the observing wavelength in mm, $\left|K\right(\mathit{\lambda }){|}^{\mathrm{2}}$ is the dielectric factor at the observing wavelength, and σ_bck(r) is the backscattering cross section in mm² as a function of particle size. S(v)_Q has units of ${\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}(\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{)}^{-\mathrm{1}}$. As the CloudCube spectra are saved in units of dBZ, the spectra are transformed to linear units using the relations dBZ $=\mathrm{10}{\log }_{\mathrm{10}}(Z/Z_{\mathrm{0}})$, Z₀=1 mm⁶ m⁻³, and Z=S(v) dv.

The radius resolution of the DSD retrieval is defined by the velocity resolution of CloudCube (0.06 ms⁻¹) and varies according to $\frac{\mathrm{d}r}{\mathrm{d}{v}_{\mathrm{t}}}$. The radius resolution varies greatly as a function of drop radius. For radii less than 50 µm, the resolution starts off coarse as $\frac{\mathrm{d}r}{\mathrm{d}{v}_{\mathrm{t}}}$ is larger. This value decreases initially as the fall speed of the smallest drops is given by Stokes' law. The resolution value increases with larger drop sizes. Plots of the terminal velocity relationships and the radius resolution at standard temperature and pressure are shown in Fig. 4. The variations seen in the radius resolution plot are likely due to measurement errors in the Gunn–Kinzer points.

Figure 4 (a) The terminal velocity relationships plotted along with interpolated radius values corresponding to the measured CloudCube spectrum velocities. The transition between the Stokes' law regime and the Gunn–Kinzer interpolated points occurs at a radius of 50 µm. (b) The radius resolution plotted as a function of the radii interpolated from the spectrum velocities.

5.1 Turbulence-free assumption

To begin with the simplest retrieval scheme, we first assume that spectra were captured in a turbulence-free environment. Ignoring attenuation error for now as well, the measured spectrum can then be modeled simply as

$$\begin{array}{}\text{(3)}& S(v+v_{\text{wind}}{)}_{\text{obs}}=\text{AS}(v{)}_{\mathrm{Q}}+{\mathit{ϵ}}_{\mathrm{s}}.\end{array}$$

With this simplification, once attenuation is corrected for, the above equation can be inverted to solve for N(r). Calculation of the spectrum error is carried out according to the analysis presented in the appendix of Hogan et al. (2005). Based on CloudCube's observing wavelength and the scale of wind speeds, both measured directly by the sondes and calculated from the spectra, we can assume that each collected Doppler spectrum, which is sampled once every 36 ms, is fully independent from the previous spectrum. Then, the spectrum error can be written as

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{s}}}{S}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{M}}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}}{\text{SNR}}}+{\displaystyle \frac{\mathrm{1}}{{\text{SNR}}^{\mathrm{2}}}}\right)},\end{array}$$

where M is the number of averaged samples in our spectra, and SNR is the signal to noise ratio of each of the points in the spectrum. For CloudCube, M=30 spectra were averaged together before being saved to disk. For attenuation, both water vapor attenuation and hydrometeor attenuation should be considered. Elevation-dependent water vapor attenuation is derived using the temperature and relative humidity measured by the radiosondes (Rosenkranz, 1998). As shown in Fig. 5a, the hydrometeor attenuation is calculated and accumulated for each elevation where data to retrieve the DSD are available as per

$$\begin{array}{}\text{(5)}& A_{\text{tot}}\left(R_N\right)=\sum _{n=\mathrm{0}}^N\mathrm{\Delta }A\left(R_n\right).\end{array}$$

Figure 5 (a) G-band attenuation contributions by both water vapor and hydrometeors. (b) Diagram visualizing how hydrometeor attenuation is accumulated over elevation.

Attenuation accumulation is determined by assuming that the DSD stays constant between successive range bins. The extinction coefficient can then be assumed to be constant between R(n) and R(n+1), making the two-way optical depth $\mathit{\tau }=\mathrm{2}\int k\left(R\right)\mathrm{d}R=\mathrm{2}k\mathrm{\Delta }R$ (the factor of 2 accounts for the round-trip distance made by a radar echo). To find k, the extinction coefficient, we need to integrate over extinction contributions from all particle sizes in the measured DSD: $k=\int k\left(r\right)\mathrm{d}r=\int N\left(r\right){\mathit{\sigma }}_{\text{ext}}\left(r\right)\mathrm{d}r$. Here, σ_ext(r) represents the extinction cross section of the particles, which is computed from the procedures described in Sect. 2.3. The incremental attenuation contribution is then determined by

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }A\left(R_n\right)\left(\text{in dB}\right)=\mathrm{10}{\log }_{\mathrm{10}}\left(e^{-\mathit{\tau }}\right)=\mathrm{10}{\log }_{\mathrm{10}}\left(e^{-\mathrm{2}k\mathrm{\Delta }R}\right).\end{array}$$

The relative importance of each of these attenuation contributions at the G-band is shown in Fig. 5b. Because the precipitation in this case is light, the attenuation is dominated by the water vapor. Figure 6a shows the 2-D DSD retrieved from the example spectrum. There are lines of increased particle number density coinciding with the radii where backscattering minima occur. These are likely unwanted artifacts in the retrieval due to some combination of not taking turbulence into account, errors in the radius–terminal velocity relationship, or errors in the radius–backscattering efficiency relationship. To retrieve a 2-D DSD that mitigates these artifacts without robust knowledge of the sources and magnitude of error, we implement a forward-modeling approach.

Figure 6 (a) DSD for each elevation calculated with the no-turbulence assumption. Dotted gray lines represent particle radii that are at minima in the backscattering efficiency. (b) The pink curve shows the DSD at 1 km, with the shaded region representing the 1σ instrument error, ϵ_s. The vertical dotted lines represent the backscattering minima, and the dashed black curve represents the limiting DSD values that CloudCube would have been able to detect at this elevation. CloudCube's sensitivity is −50 dBZ at 1 km. N(r) is in units of ${\mathrm{m}}^{-\mathrm{3}}{\mathrm{m}}^{-\mathrm{1}}$.

5.2 Forward modeling of turbulent spectrum

The forward-modeling approach attempts to retrieve a vector that best represents the DSD. The model first uses an initial vector N that represents the DSD to create an idealized spectrum using Eq. (2). A log-linear best fit is applied to the DSD calculated under the turbulence-free assumption for initialization. This spectrum is then smoothed with a blurring kernel that represents the effect of the turbulence to compute a spectrum that can be best compared with the measured spectrum. We minimize a loss function to find the most likely vector N. A diagram of this is depicted in Fig. 7.

Figure 7 Block diagram explaining the basic steps of the forward model.

Our forward model needs an estimate of the turbulence scale at every height. To do this, the framework described in O'Connor et al. (2005), which uses large-scale turbulence to estimate smaller-scale turbulence as such, is utilized:

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{\text{v,turb}}^{\mathrm{2}}={\mathit{\sigma }}_{\text{v,air}}^{\mathrm{2}}\left({\displaystyle \frac{L_{\text{small}}^{\mathrm{2}/\mathrm{3}}}{L_{\text{large}}^{\mathrm{2}/\mathrm{3}}-L_{\text{small}}^{\mathrm{2}/\mathrm{3}}}}\right).\end{array}$$

The term ${\mathit{\sigma }}_{\text{v,air}}^{\mathrm{2}}$ is the variance in the vertical wind speeds, while the terms L_small and L_large represent the small and large length scales of the turbulence, respectively. The L terms are dependent on the horizontal wind in the line of sight of the observation (U), the range of interest (R), the beamwidth of the radar (θ), and the averaging time (t):

$$\begin{array}{}\text{(8)}& L=Ut+\mathrm{2}R\sin \left({\displaystyle \frac{\mathit{\theta }}{\mathrm{2}}}\right).\end{array}$$

The short timescale, t_small, is the time between successive spectra (1.13 s for CloudCube). The long timescale, t_large, is the total time over which the variance of the vertical winds is calculated. The variation in the vertical winds in the time vicinity of the example spectrum is shown in Fig. 8a. This figure also shows the horizontal wind speeds captured by the sondes at the same time as the example spectrum and the final turbulence derived from those values. The turbulence scale is very similar to the velocity resolution of the spectra, so smoothing will have a relatively small effect on the DSD retrieval.

Figure 8 (a) Depiction of how the vertical wind varies with time for a few different elevations. The dotted lines represent ±0.25 m s⁻¹ for each elevation. Plots of (b) the horizontal wind speeds measured by the sondes and (c) the final calculated turbulent broadening, in units of m s⁻¹.

The most obvious choice of loss function for our forward model would be least squares. However, as we noted in the previous section, factors beyond turbulence correction lead to unphysical artifacts in the retrieved DSDs. One way to retrieve a smooth DSD would be to impose a functional form for the DSD, such as the modified gamma distribution (Deirmendjian, 1969). We take another approach. To encourage the retrieval of smoother and more physically realistic DSDs, we utilize a regularized least squares loss function:

$$\begin{array}{}\text{(9)}& L=\sum \left[{\displaystyle \frac{(S_{\text{calc},i}-S_{\text{meas},i}{)}^{\mathrm{2}}}{{\mathit{ϵ}}_{\mathrm{s},i}^{\mathrm{2}}}}+\mathit{\alpha }(N_{i+\mathrm{1}}-N_i{)}^{\mathrm{2}}\right].\end{array}$$

The first part of the loss function is a classic least squares loss. The second term represents the regularization. We use the total squared variation (TSV) regularizer, represented by $(N_{i+\mathrm{1}}-N_i{)}^{\mathrm{2}}$. This regularizer was first introduced by Kuramochi et al. (2018) as a way to enforce smoothness in 2-D imaging retrievals. The same principle applies to 1-D vectors, as penalizing the squared difference between adjacent points in a DSD favors a smoothly varying vector. The term α represents the regularizer weight, which determines how strictly we want to enforce vector smoothness. A large amount of regularization, meaning a larger value for α, will retrieve highly smoothed DSD vectors. We show a demonstration of this in Fig. 9.

Figure 9 (a) DSDs retrieved from the 1-D spectrum shown in Fig. 2d using five different regularizer weights. The grey curves plotted underneath the colored plots show the DSD computed from the turbulence-free assumption for comparison. The dotted vertical lines represent particle radii that are at minima in the backscattering efficiency. N(r) is in units of ${\mathrm{m}}^{-\mathrm{3}}{\mathrm{m}}^{-\mathrm{1}}$. (b) The smoothed spectra derived from each of the DSDs. The grey curves underneath show the measured spectrum.

We see that the pink curve (representing the highest amount of regularization we explored) produces a very smooth DSD and, in turn, a very smooth final spectrum. In Fig. 10, we show a 2-D DSD for the example spectrum, retrieved using the highest regularizer weight we tested, 10⁴. Compared to the turbulence-free DSD retrieval, we see a significant reduction in sharp gradients, though they are not completely eliminated. A primary issue with using regularized least squares, however, is that we currently have no way to validate our choices of regularization. As the problem is over-constrained, different regularizer choices can yield final results with similar least squares errors. Without any ground truth data to train for the correct regularizer weight, we can only place confidence in the general shape and statistical properties derived from the DSD.

Figure 10 DSD for each elevation calculated using a regularized forward model and taking into account turbulent broadening. Dotted grey lines represent particle radii that are at minima in the backscattering efficiency. N(r) is in units of ${\mathrm{m}}^{-\mathrm{3}}{\mathrm{m}}^{-\mathrm{1}}$.

6 Discussion

6.1 Sources of uncertainty

Thus far in this work, we have not propagated errors to derive a final estimate of the uncertainty for the retrieved DSD values. While we acknowledge the importance of robust uncertainty analysis, many of the sources of error present are unable to be parameterized for easy propagation through all of the steps of this retrieval. Therefore, this section lists each potential source of error and qualitatively comments on the sensitivity of the final DSD to this error.

The most significant source of uncertainty is the estimation of the vertical wind. Confidence in the vertical wind value is dependent on the confidence in the minimum finding. The Gaussian mixture model used in minimum classification does provide an in-distribution probability for each classified point; however, there exists no clean way to relate this probability value to an uncertainty in m s⁻¹. Using minimum depth as a proxy for uncertainty is also not feasible, as there is no guarantee that the deepest minimum in a spectrum actually corresponds to a backscattering notch. The terminal velocity is seen to affect the quiet-air estimate of the first minimum; however, several relationships are seen to provide a consistent theoretical first minimum velocity. Hence, utilizing only the first minimum to compute vertical wind reduces uncertainty in the final wind estimate.

Offsets in the radiosonde release time and the measurement time of the G-band spectra may introduce a small uncertainty as well. Important properties either directly measured by the sonde or computed from sonde measurements include air density, air temperature, gaseous attenuation, and horizontal wind. The largest source of uncertainty likely arises from temporal variation in the gaseous attenuation profile, which may fluctuate by up to a few decibels if the humidity changes significantly. The horizontal wind measurement directly affects the turbulent broadening calculation. For the case of EPCAPE, the turbulent broadening is minimal, so temporal variations in the horizontal wind likely would not have a large effect. However, for measurements taken at much windier sites with large variations in the horizontal winds, a confident estimate on the turbulent broadening may be difficult to calculate. Air temperature and density are used to adjust the scattering functions and the terminal velocity functions, respectively, though these only change by a few percent over the course of 1–2 km. Without a priori knowledge of the timescales of the variability of these quantities, it is impossible to estimate an uncertainty value to be included in the retrieval process.

Another source of uncertainty that is challenging to quantify is the retrieval uncertainty. The uncertainty in the retrieved DSD introduced by regularization choices and initial condition choices cannot be parameterized. The measurement uncertainties for the spectrum data are known, but part of the retrieval uncertainty depends on the solution space for each model. The scipy.optimize function utilized in this retrieval provides an estimated covariance matrix; however, the high dimensionality of this problem likely leads to a complex solution space that is likely not explored thoroughly enough for good covariance estimates. Monte Carlo optimization techniques perform a more thorough search of the solution space generally, but such methods are too computationally expensive for high-dimensionality problems such as ours. Thus, a full error estimation for this forward model is highly non-trivial and is beyond the scope of this work.

6.2 Validation with co-observing instruments

Despite not fully understanding the uncertainty in the final retrieval, a first-order understanding of the accuracy of the retrieved DSD can be found by utilizing data from both the ARM KaZR and the VDIS instruments. Using the DSD depicted in Fig. 10 along with Eq. (2), we can calculate a predicted Ka-band spectrum for each elevation. We use T-matrix scattering coefficients at the Ka-band to calculate both the backscattering and the extinction cross sections. We correct for water vapor and hydrometeor attenuation in our theoretical spectrum before integrating across velocity to compute a single reflectivity value for each elevation. This reflectivity can be directly compared to the reflectivity measured by KaZR, as shown in Fig. 11a. We can see that for the example spectrum, the predicted Ka-band reflectivity matches the KaZR fairly well, generally within a few decibels. The presence of larger drops, which G-band instruments are not as sensitive to but Ka-band instruments are, may affect the accuracy of the predicted Ka-band reflectivity. Uncertainties in the hydrometeor attenuation also increase with elevation, potentially leading to higher inaccuracies in the Ka-band predictions as well. Still, the general consistency between the two curves gives us some confidence in the quality of our retrievals.

Figure 11 (a) Ka-band reflectivity for each elevation calculated using the retrieved regularized DSD plotted with the Ka-band data collected by the ARM KaZR instrument. (b) Comparison between the lowest-elevation DSDs retrieved from the G-band spectra (dashed lines) and DSDs measured by the VDIS at the same time (solid line). The dashed dark blue line represents a measurement of the DSD at 200 m above the surface, while the dashed light blue line represents a measurement 100 m above the surface.

Figure 12 Visualization of estimation uncertainty in precipitation properties due to errors in the vertical wind retrieval and choices of regularization in the DSD retrieval. The solid lines represent parameters derived with the winds shown in Fig. 3 and with a regularizer weight of α=10⁴. The dashed line shows properties derived with the same winds but with very little regularization (α=1). The shaded region represents properties retrieved with α=10⁴ but assuming a ±0.12 m s⁻¹ deviation from the measured vertical wind. This represents two spectrum bins away from the best estimate.

We can also use direct measurements of the DSD taken by VDIS and compare them to the lowest-elevation retrieved DSD we have available for the same time. The video disdrometer measures number densities in 0.1 mm radius increments, with a limiting drop size of 0.05 mm. However, the accuracy of the VDIS measurements below 0.1 mm is reduced due to the instrument struggling more to distinguish between smaller drop sizes. Additionally, the VDIS only saves data in 1 min increments, and unfortunately we only have a few minutes of data for which DSDs are able to be retrieved in the times adjacent to the example spectrum. Thus, there are only two coincident times between the CloudCube measurements and the VDIS measurements, plotted in Fig. 11. We see a similarity between the retrieved DSD and the VDIS-measured DSD, with discrepancies being the highest at the smallest drop sizes. Because of the significant fall time of the small droplets from the lowest DSD elevation (typically around 50 m) to the ground, comparing measurements with the same time stamps may involve comparing slightly different drop populations. However, with a very limited amount of G-band data of good enough quality to retrieve DSDs (typically only of the order of a few minutes) and the slow sampling time of the VDIS, it is challenging to compare measurements with a sufficient lag time to account for the fall time.

Figure 13 Compilation of statistics that can be calculated from retrieved DSDs, shown for five G-band spectra spaced 35 s apart. Times shown in the left column are seconds elapsed since 30 March 2023 at 17:10:18 UTC. Gaps in the DSD are due to minima corresponding to the first backscattering minimum not being available at that elevation. In the bottom row in particular, this can be seen for elevations where the minimum is actually below the noise floor and is therefore not detected by our algorithm. Spectra are shown in dBZ, using the same color bar as used previously. N(r) is in units of ${\mathrm{m}}^{-\mathrm{3}}{\mathrm{m}}^{-\mathrm{1}}$, using the same color bar as used previously.

6.3 Estimating bulk precipitation properties

Despite the various uncertainties described above, we can derive bulk properties of the DSDs, which are easier to use and should be more robust to uncertainties. Here we derive four bulk properties of the distribution. Figure 12 shows the plots of the mass-weighted mean radius, total number density, liquid water content, and precipitation rate derived from the DSD of the example spectrum. The mass-weighted mean radius, R_m, is calculated as

$$\begin{array}{}\text{(10)}& R_{\mathrm{m}}={\displaystyle \frac{\int N\left(R\right)M\left(R\right)R\mathrm{d}R}{\int N\left(R\right)M\left(R\right)\mathrm{d}R}},\end{array}$$

where $M\left(R\right)={\mathit{\rho }}_{\mathrm{w}}V\left(R\right)={\mathit{\rho }}_{\mathrm{w}}(\mathrm{4}/\mathrm{3})R^{\mathrm{3}}$ is the mass of a water droplet with radius R, where ρ_w is the density of liquid water and V(R) is the volume of a drop with radius R. The total number density is simply calculated as $N_{\text{tot}}=\int N\left(R\right)\mathrm{d}R$, the liquid water content is calculated as $\text{LWC}=\int N\left(R\right)M\left(R\right)\mathrm{d}R$, and the precipitation rate is calculated as $P=\int N\left(R\right)V\left(R\right)v\left(r\right)\mathrm{d}R$. Profiles of these properties are shown for the lowest and highest explored values of the regularization weight, α, and for a ±0.12 m s⁻¹ error in the vertical wind speed (representing ± two bins in the measured Doppler spectrum). We see that while these errors may affect the details of the binned DSDs, the bulk precipitation properties are relatively insensitive to the choice of regularization weight and fairly robust to vertical air motion uncertainties.

Figure 13 shows the plots of these precipitation properties derived from the DSDs of five different spectra, spaced 35 s apart. This figure highlights the rapid timescale of variability present in these drizzling systems. Coarse sampling times in measurements of the precipitation properties are at risk of missing important details in the cloud and precipitation processes.

7 Conclusions

We have presented a retrieval methodology to derive vertical wind and the precipitation DSD in light rainfall from a nadir-pointing G-band Doppler spectrum. This work extends the methods developed for the W-band to lighter rainfall than has been possible to date. The G-band retrievals work well for light precipitation because the first Mie notch occurs near a radius of 334 µm, thereby enabling accurate estimation of the wind speed for very light precipitation rates. Furthermore, the precipitation water contents are very small, so the attenuation from condensed water is insignificant relative to the gaseous attenuation. As pointed out by Courtier et al. (2024), the method demonstrated here would optimally be combined with multi-frequency W- and K-band Doppler spectra (e.g., Tridon and Battaglia, 2015) to seamlessly extend from the lightest precipitation events to the heaviest precipitation events.

There are residual uncertainties in the binned DSD due to inaccuracies in the droplet fall velocity and drop obliquity relationships, which appear as ripples near the location in the spectrum where Mie notches are present. We expect some uncertainty contribution from factors such as errors in the minimum-finding routine, errors in the estimation of turbulent broadening effects, errors in computed gaseous attenuation, and discrepancies in air temperature and density, which propagate into errors in the DSD retrieval. Nevertheless, the bulk statistics of the DSD, such as the water content, number concentration, precipitation rate, and mass-weighted mean size, are relatively robustly derived. We note that the retrievals presented here are for high-SNR cases with clearly defined notches in the spectrum. Future work will refine the algorithm presented here to work for spectra with less clearly defined spectral features.

With the growing number of G-band radar observations (including CloudCube G-band's ongoing participation in the Cloud And Precipitation Experiment at kennaook; Mace et al., 2023), the Doppler spectral retrieval method offers the unique potential to provide profiles of light rainfall and drizzle in stratocumulus and shallow cumulus clouds relative to approaches centered on radar reflectivity.