Spaceborne differential absorption radar water vapor retrieval capabilities in tropical and subtropical boundary layer cloud regimes

Differential absorption radar (DAR) near the 183 GHz water vapor absorption line is an emerging measurement technique for humidity profiling inside of clouds and precipitation with high vertical resolution, as well as for measuring integrated water vapor (IWV) in clear-air regions. For radar transmit frequencies on the water line flank away from the highly attenuating line center, the DAR system becomes most sensitive to water vapor in the planetary boundary layer (PBL), which is a region of the atmosphere that is poorly resolved in the vertical by existing spaceborne humidity and temperature profiling instruments. In this work, we present a high-fidelity, end-to-end simulation framework for notional spaceborne DAR instruments that feature realistically achievable radar performance metrics and apply this simulator to assess DAR's PBL humidity observation capabilities. Both the assumed instrument parameters and radar retrieval algorithm leverage recent technology and algorithm development for an existing airborne DAR instrument. To showcase the capabilities of DAR for humidity observations in a variety of relevant PBL settings, we implement the instrument simulator in the context of large eddy simulations (LESs) of five different cloud regimes throughout the trade-wind subtropical-to-tropical cloud transition. Three distinct DAR humidity observations are investigated: IWV between the top of the atmosphere and the first detected cloud bin or Earth's surface; in-cloud water vapor profiles with 200 meter vertical resolution; and IWV between the last detected cloud bin and the Earth's surface, which can provide a precise measurement of the sub-cloud humidity. We provide a thorough assessment of the systematic and random errors for all three measurement products for each LES case and analyze the humidity precision scaling with along-track measurement integration. While retrieval performance depends greatly on the specific cloud regime, we find generally that for a radar with cross-track scanning capability, in-cloud profiles with 200 m vertical resolution and 10 %–20 % uncertainty can be retrieved for horizontal integration distances of 100–200 km. Furthermore, column IWV can be retrieved with 10 % uncertainty for 10–20 km of horizontal integration. Finally, we provide some example science applications of the simulated DAR observations, including estimating near-surface relative humidity using the cloud-to-surface column IWV and inferring in-cloud temperature profiles from the DAR water vapor profiles by assuming a fully saturated environment.

1 Introduction

Cloud morphology and precipitation depend sensitively on the three-dimensional distributions of water vapor and temperature, especially within the planetary boundary layer (PBL) where most convective initiation occurs. Existing spaceborne sensors have limited ability to sample water vapor and temperature in the PBL with high spatial resolution, with increased difficulty inside of cloudy and precipitating volumes for passive infrared and microwave sounders (Wulfmeyer et al., 2015; Stevens et al., 2017; Sahoo et al., 2015). The 2017 Decadal Survey for Earth Science and Applications from Space (NASEM, 2018) has recommended the incubation of technologies that enable improved spaceborne measurements of PBL thermodynamics during the current decade. Active humidity sounding approaches, including differential absorption lidar (DIAL) and differential absorption radar (DAR), offer new potential spaceborne solutions for providing high-vertical-resolution water vapor profiles in clear-sky and cloudy regions, respectively (Nehrir et al., 2017), with a typical vertical resolution target of 200 m. The advantage of the active DIAL/DAR approaches is to constrain the vertical variability in water vapor desired by the decadal survey (NASEM, 2018)​​​​​​​. It is anticipated that these observations would exist within a larger observing system where the requisite horizontal variability is observed by passive infrared and microwave systems. However, because DAR is a relatively new measurement approach with limited instrument deployment (Roy et al., 2020, 2021; Cooper et al., 2020) and simulation (Lebsock et al., 2015; Millán et al., 2016; Battaglia and Kollias, 2019; Millán et al., 2020; Schnitt et al., 2020) heritage, there is a need to critically assess the measurement capabilities of notional spaceborne DAR systems with detailed radar simulations from an orbital altitude, in addition to continued assessment of observational capabilities from airborne platforms.

In this work, we expand significantly on previous spaceborne DAR instrument simulation efforts to provide a detailed assessment of DAR retrieval capabilities in the context of different regimes throughout the trade-wind marine cloud transition from stratocumulus to deep convection. First, we utilize five different large-eddy simulations (LESs) that represent distinct cloud regimes throughout the transition, and we perform forward radar simulations that include detailed effects such as realistic electromagnetic scattering from ice particles, non-uniform beam filling leveraging the high spatial resolution of the LES, and multiple scattering. Second, we build on retrieval approaches developed previously for real DAR observations (Roy et al., 2020) and implement a new least-squares retrieval algorithm that mitigates retrieval bias using an improved spectral fitting function and retrieves the entire vertical humidity profile as part of a single optimization procedure. The new formulation employs an improved humidity interpolation function that allows for a flexible retrieval of water vapor between arbitrarily spaced cloud range bins, which is important given the inherently sparse sampling of radar reflectivity profiles. The result is a seamless retrieval of both high-vertical-resolution (200 m) in-cloud profiles and integrated water vapor measurements in clear-air columns (e.g., between the cloud base and the surface). Importantly, this retrieval does not require specification of a prior distribution for the humidity profile, nor does it impose any regularizing constraints on the form of the retrieved profile as in Roy et al. (2020). The linearity of the retrieval algorithm's forward model in the state vector allows for a straightforward transformation of random measurement error to retrieved state uncertainty, which is then analyzed as a function of along-track averaging distance to assess the achievable humidity precision in the different cloud scenarios.

For these simulations, we place emphasis on prescribing instrument performance metrics that are technologically feasible for a spaceborne radar deployment in the next decade and will therefore provide an honest assessment of the platforms capabilities. In fact, only recently has it become possible to develop cloud radars at frequencies above 100 GHz with useful measurement sensitivity for atmospheric science studies. The first in this new line of G-band cloud radars, the Jet Propulsion Laboratory's Vapor In-cloud Profiling Radar (VIPR) (Cooper et al., 2020; Roy et al., 2020, 2021; Lamer et al., 2021), achieves a noise-equivalent reflectivity, a number which roughly dictates the minimum detectable cloud signal, of −40 dBZ at 1 km range using a range resolution of 15 m. If an identical system were deployed in low-earth orbit with an altitude of 400 km and range resolution of 200 m, the resulting noise-equivalent reflectivity would become about 0 dBZ. This sensitivity level is 16 dB worse than the equivalent figure for CloudSat's 94 GHz Cloud Profiling Radar and, when combined with the increased attenuation at G-band in clear air and within clouds and precipitation, suggests that such a system would have minimal cloud sampling ability, especially in the PBL. Thus, significant improvements to the VIPR transceiver sensitivity are necessary for spaceborne implementation.

The radar performance metrics assumed here are consistent with a system that employs a long-duration pulse with linear frequency modulation and pulse compression with very large time-bandwidth product, similar to that in the recent successful Radar in a CubeSat (RainCube) mission (Peral et al., 2019). We assume that the instrument is deployed on a traditional medium-sized satellite bus, similar in scale to the CloudSat platform, that can afford a large, solid reflector antenna of diameter D_a=2 m and a high-power, high-duty-cycle transmitter with a peak output power of 200 W. Such transmitter performance is consistent with a vacuum-electronics-based, traveling-wave tube amplifier that is currently under development, and which is based on a heritage implementation at 231–235 GHz for a synthetic aperture radar application (Basten et al., 2016). The three transmit frequencies located at 155.5, 168.0, and 174.8 GHz are carefully chosen to lie in bands that do not feature international transmission restrictions (NTIA, 2015) and correspond to bands that are currently used in the VIPR system. The middle DAR channel used here differs slightly from the 167 GHz lower frequency channel used in past two-frequency DAR measurements with VIPR (Roy et al., 2020). This new location is chosen to minimize covariance between fit parameters in the three-frequency retrieval (see a₂ and a₃ in Eqs. 6 and 7). Finally, while the desired vertical humidity profile resolution of 200 m suggests that the radar could have a range resolution of 200 m, we find that it is necessary to oversample the range dimension by a small factor in order to remove the humidity measurement bias that arises from non-uniform filling of range bins by cloud and precipitation particles. Therefore, we assume a range resolution of 50 m, which can be achieved with a moderate radar pulse modulation bandwidth of around 6 MHz. We note that the approach implemented here does not include specific aspects of the radar range weighting function (RWF) and is instead equivalent to assuming an RWF that is a top-hat function of width 50 m.

2 Methods

2.1 LES case studies

The synthetic observations used in this study were produced by means of large-eddy simulation (LES), which is a well-known method for realistic high-resolution modeling of multi-phase three-dimensional non-hydrostatic turbulent atmospheric flows (Stevens and Lenschow, 2001). We choose five canonical cases covering the transition from the subtropics to the tropics over ocean. The simulated convective regimes were a focus of the previous field campaigns: DYCOMS-II and VOCALS (both on stratocumulus), BOMEX (non-precipitating shallow Cu), RICO (shallow Cu), and GATE (deep convection). Results of those campaigns served as the initial and boundary conditions for the LES. Details of the field campaigns and the modeling setups are provided in Table 1 and the references therein.

Stevens et al. (2005)Berner et al. (2011)Siebesma et al. (2003)vanZanten et al. (2011)Gentine et al. (2016)Matheou and Chung (2014)Khairoutdinov and Randall (2003)Matheou and Chung (2014)Khairoutdinov and Randall (2003)Skamarock et al. (2008)(Morrison et al., 2005)(Morrison et al., 2005)(Lin et al., 1983)

Table 1 Summary of LES parameters for the five cases.

^* These ranges for GATE correspond to values at 24 and 6 h into the simulation, respectively, with a monotonic decrease over this time window.

For each simulation, the cascade of turbulent motions develops within the domain, with the boundary layer dynamics driven by surface fluxes of latent and sensible heat and the vertical shear of horizontal wind, modified by radiation (either interactive for stratocumulus or prescribed for other cases) and large-scale advective tendencies (for shallow Cu cases). For the shallow precipitating cases, two-moment microphysics provides the information about both mixing ratios and droplet number concentrations of the precipitating and non-precipitating forms of water. For deep convection, bulk one-moment microphysics is used that represents six classes of water: vapor, cloud liquid, cloud ice, warm rain, snow, and graupel. Horizontal resolutions range from 5 m for DYCOMS-II to 100 m for GATE. Three-dimensional LES outputs are combined with one-dimensional MERRA-2 reanalysis data (Gelaro et al., 2017) aloft, producing a full-depth atmospheric column that is the input for the radiative transfer model. We utilize the grid point from MERRA-2 that is closest to the model domain, where the MERRA-2 latitudinal and longitudinal resolutions are 0.5^∘ and 0.625^∘, respectively, and the thermodynamic profiles used here are given on 42 vertical levels.

2.2 Radar forward model

We implement the radar forward model on a given LES atmospheric state in three distinct steps: (1) calculation of volumetric and surface scattering parameters for each of the DAR channel frequencies at the LES model resolution; (2) interpolation in the vertical dimension to an equally spaced grid that is much finer than the radar range resolution, and subsequent computation of the time-dependent radiative transfer solution including the effects of multiple scattering; and (3) processing of the ideal, high-resolution radar backscatter quantities to produce observations consistent with the assumed radar instrument parameters. The important elements of the single-scattering calculations in step (1) are described in Appendix A, while steps (2) and (3) are described below.

The single-scattering properties of the atmospheric gases, cloud and precipitation hydrometeors, and ocean surface (see Appendix A) at the LES spatial resolution serve as inputs to the multiple-scattering simulator of step (2), for which we utilize the time-dependent, two-stream implementation of Hogan and Battaglia (2008). In order to assess systematic error in the DAR retrieval resulting from non-uniform range-bin filling, we interpolate the scattering property fields from the LES grid resolution – which has a variable spacing that increases with height – to an equally spaced vertical grid with 10 m resolution. This ensures that the instrument range resolution of Δr=50 m is sufficiently oversampled in the forward simulator throughout the LES domain. At this stage, the only instrument parameter that is required to be specified for the multiple-scattering calculation is the angular width of the radar beam, which we calculate from the antenna diameter D_a assuming a Gaussian beam shape and an 11 dB beam taper. The output of the multiple-scattering calculation is the observed, or apparent, volume backscatter coefficient ${\mathit{\eta }}_{\mathrm{obs}}(f,x,y,z)$ as a function of transmit frequency at the LES horizontal resolution (Δx,Δy) and interpolated vertical resolution (10 m).

Given the forward-simulated, observed, volume backscatter coefficient, ${\mathit{\eta }}_{\mathrm{obs}}(f,x,y,z)$, we begin step (3) by averaging over the instrument horizontal beam footprint and range resolution. Here it is assumed that the along-track direction corresponds to the y dimension. Note that in this step of resolution degradation we are leveraging the fine spatial resolution of LES to provide realistic representation of non-uniform beam-filling effects both in terms of the horizontal footprint and within single range bins. Given an antenna footprint $\left|F\right(f,\mathit{\theta },\mathit{\varphi }){|}^{\mathrm{2}}=|F(f,\mathit{r}){|}^{\mathrm{2}}$, which we express in the interpolated Cartesian coordinate system, we calculate the footprint-averaged, observed backscatter after an along-track integration time of T_int from

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{obs}}(f,x_i,y_j,z_k)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }r{R}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{eff}}{T}_{\mathrm{int}}}}\underset{V_{ijk}}{\int }{d}^{\mathrm{3}}{\mathit{r}}^{\prime }{\mathit{\eta }}_{\mathrm{obs}}(f,{\mathit{r}}^{\prime })\\ & \times \underset{\mathrm{0}}{\overset{T_{\mathrm{int}}}{\int }}\mathrm{d}{t}^{\prime }\left|F\right(f,x^{\prime },y^{\prime }-v_g{t}^{\prime }){|}^{\mathrm{4}}.\end{array}\end{array}$$

Here the orbital altitude R appears via the solid-angle normalization factor, where we have substituted $\mathrm{d}{\mathrm{\Omega }}^{\prime }=\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime }(R-z_k{)}^{-\mathrm{2}}\approx \mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime }{R}^{-\mathrm{2}}$ throughout the LES domain for simplicity since the true range from the satellite R−z_k varies by only a few percent. The integration volume V_ijk centered on the observation point $(x_i,y_j,z_k)$ is defined by $-\mathrm{\Delta }{x}_{\mathrm{3}\mathrm{dB}}/\mathrm{2}\le x^{\prime }-x_i\le \mathrm{\Delta }{x}_{\mathrm{3}\mathrm{dB}}/\mathrm{2}$, $-\mathrm{\Delta }{y}_{\mathrm{3}\mathrm{dB}}\left(T_{\mathrm{int}}\right)/\mathrm{2}\le y^{\prime }-y_j\le \mathrm{\Delta }{y}_{\mathrm{3}\mathrm{dB}}\left(T_{\mathrm{int}}\right)/\mathrm{2}$, and $-\mathrm{\Delta }r/\mathrm{2}\le z^{\prime }-z_k\le \mathrm{\Delta }r/\mathrm{2}$, where Δx_3 dB is the 3 dB full width of the two-way beam pattern $\left|F\right(f_{\mathrm{0}},\mathit{r}){|}^{\mathrm{4}}$ at the lowest frequency, and Δy_3 dB(T_int) is the 3 dB full width of the time-averaged beam footprint in the along-track direction. Finally, the effective two-way solid angle Ω_eff is the integral of the time-averaged two-way beam pattern over the restricted angular domain defined by V_ijk. Averaging the observed normalized radar cross section (NRCS) over the satellite footprint to get ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{obs}}^{\mathrm{0}}(f,x_i,y_j)$ is performed by making the replacement ${\mathit{\eta }}_{\mathrm{obs}}(f,{\mathit{r}}^{\prime })\to \mathrm{\Delta }r\mathit{\delta }\left(z^{\prime }\right){\mathit{\sigma }}_{\mathrm{obs}}^{\mathrm{0}}(f,x^{\prime },y^{\prime })$ in Eq. (1), where δ(z^′) is the Dirac delta function.

We evaluate Eq. (1) for each footprint centroid (x_i,y_j) by performing a discrete summation at the interpolated model resolution. The radar sampling horizontal grid is defined by $x_{i+\mathrm{1}}-x_i=\mathrm{\Delta }{x}_{\mathrm{3}\mathrm{dB}}$ and $y_{i+\mathrm{1}}-y_i=v_g{T}_{\mathrm{int}}$, implying that adjacent footprints in the cross-track dimension contain no common LES pixels, while sequential along-track footprints feature realistic overlap. Note that we do not assume that the radar can measure multiple, simultaneous footprints in the cross-track dimension but instead are maximizing use of the rectangular LES domain to gather as many statistical samples as possible for eventual retrieval performance evaluation. Finally, for comparison of radar-retrieved humidity profiles with those in the LES, we average the model humidity field over the time-averaged, two-way antenna pattern as well.

Next, we calculate the observed reflectivity from $Z_{\mathrm{obs}}\left(f\right)={\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{obs}}\left(f\right)c^{\mathrm{4}}{f}^{-\mathrm{4}}{\mathit{\pi }}^{-\mathrm{5}}\left|K_{\mathrm{w}}\right(f){|}^{-\mathrm{2}}$, where c is the speed of light, $K_{\mathrm{w}}\left(f\right)=\left[{\mathit{\epsilon }}_{\mathrm{w}}\right(f)+\mathrm{2}{]}^{-\mathrm{1}}[{\mathit{\epsilon }}_{\mathrm{w}}\left(f\right)-\mathrm{1}]$, and ε_w(f) is the dielectric constant of pure water evaluated at 280 K. At this point we must specify the instrument parameters that determine the minimum detectable signal and random measurement uncertainty. These values are listed in Table 2. As previously discussed, we assume a transmitter with a pulse-averaged power of 200 W and duty cycle of 25 %. Then, in order to determine the single-pulse radar sensitivity, we must specify the pulse duration τ_p and compare this with the measurement decoherence time. Generally speaking for a pulse compression radar that is not performing Doppler velocity estimation, it is desirable to have τ_p as large as possible while maintaining coherent measurement integration to maximize the signal-to-noise ratio (SNR). However, for a spaceborne radar in low-earth orbit, the fast transit of the beam across a region of interest combined with the fact that the measurement relative error is at minimum ${\mathit{\sigma }}_Z{Z}_{\mathrm{obs}}^{-\mathrm{1}}=N_{\mathrm{i}}^{-\mathrm{1}/\mathrm{2}}$, where N_i is the number of statistically independent pulses transmitted at a given frequency, implies that τ_p should not be so large that insufficient incoherent averaging occurs.

Table 2 High-level airborne VIPR system parameters.

^* The minimum detectable reflectivity is evaluated at the surface range.

Because we expect measurement decorrelation to be dominated in general by broadening of the returned signal spectrum due to the satellite-motion-induced Doppler effect, we can quantify the decoherence timescale by analyzing the resulting pulse-to-pulse correlations. The relevant equations are discussed in several references (Tanelli et al., 2002; Hogan et al., 2005) and are revisited in Appendix B, with the end result that the pulse-to-pulse decoherence timescale for radar signal power measurements, sometimes referred to as the “time to independence”, is given by ${\mathit{\tau }}_{\mathrm{i}}=(\mathrm{2}\sqrt{\mathit{\pi }}{\mathit{\sigma }}_f{)}^{-\mathrm{1}}$, where ${\mathit{\sigma }}_f=v_{\mathrm{sat}}{\mathit{\theta }}_{\mathrm{0}}{\mathit{\lambda }}^{-\mathrm{1}}$ is the standard deviation of the Doppler-effected spectral density for a stationary target that fills the radar beam. However, because of the inherent scaling ${\mathit{\theta }}_{\mathrm{0}}\propto \mathit{\lambda }/D_{\mathrm{a}}$, the Doppler spectral width in frequency space is independent of the radar frequency, and for our assumed hardware parameters the coherence time becomes ${\mathit{\tau }}_{\mathrm{i}}\approx \mathrm{0.40}D{v}_{\mathrm{sat}}^{-\mathrm{1}}=\mathrm{104}$ µs. Given N_p sequential pulses with an inter-pulse spacing of T_p, the number of statistically independent pulses for radar signal power estimation N_i is determined from (Doviak and Zrnić, 1993)

$$\begin{array}{}\text{(2)}& \mathit{\xi }({\mathit{\tau }}_{\mathrm{i}},T_{\mathrm{p}})\equiv {\displaystyle \frac{N_{\mathrm{p}}}{N_{\mathrm{i}}}}=\mathrm{1}+\mathrm{2}\sum _{m=\mathrm{1}}^{N_{\mathrm{p}}-\mathrm{1}}\left(\mathrm{1}-{\displaystyle \frac{m}{N_{\mathrm{p}}}}\right)\exp \left(-{\displaystyle \frac{m^{\mathrm{2}}{T}_{\mathrm{p}}^{\mathrm{2}}}{{\mathit{\tau }}_{\mathrm{i}}^{\mathrm{2}}}}\right).\end{array}$$

Because the radar cycles between widely separated transmit frequencies after each pulse, the relevant inter-pulse spacing in this case is T_p=N_fT_rep, where N_f=3 is the number of frequency channels. Evaluating Eq. (2) using the parameters in Table 2, we find that ξ≈1 or N_p=N_i. Next, we calculate the relative uncertainty in the reflectivity measurement due to random error according to (Papoulis, 1965; Doviak and Zrnić, 1993; Torres, 2001)

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{{\mathit{\sigma }}_Z}{Z_{\mathrm{obs}}}}={\displaystyle \frac{\mathrm{1}}{\sqrt{N_{\mathrm{p}}}}}{\left[\mathit{\xi }({\mathit{\tau }}_{\mathrm{i}},T_{\mathrm{p}})+{\displaystyle \frac{\mathrm{2}}{\mathrm{SNR}}}+{\displaystyle \frac{\mathrm{1}}{{\mathrm{SNR}}^{\mathrm{2}}}}\right]}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where SNR is the measurement signal-to-noise ratio, which is easily determined in our simulations by comparing the calculated radar echo power with the prescribed receiver noise $P_{\mathrm{n}}=k_{\mathrm{B}}{T}_{\mathrm{B}}F{\mathit{\tau }}_{\mathrm{p}}^{-\mathrm{1}}$. Here k_B is Boltzmann's constant, T_B is the scene brightness temperature, assumed to be 280 K, and F is the noise factor related to the noise figure through $F={\mathrm{10}}^{\mathrm{NF}/\mathrm{10}}$. Note that Eq. (3) assumes that the noise power contribution to the detected signal can be removed without inflating the measurement variance. The noise-equivalent reflectivity Z_NE and surface NRCS ${\mathit{\sigma }}_{\mathrm{NE}}^{\mathrm{0}}$ are defined as the respective observed values that correspond to SNR=1, while the minimum detectable reflectivity and NRCS improve upon this figure as a result of incoherent averaging according to $\mathrm{dB}{Z}_{\min }=\mathrm{dB}{Z}_{\mathrm{NE}}-\mathrm{10}{\log }_{\mathrm{10}}\left(\sqrt{N_{\mathrm{i}}}\right)$ and ${\mathit{\sigma }}_{\min }^{\mathrm{0}}\left(\mathrm{dB}\right)={\mathit{\sigma }}_{\mathrm{NE}}^{\mathrm{0}}\left(\mathrm{dB}\right)-\mathrm{10}{\log }_{\mathrm{10}}\left(\sqrt{N_{\mathrm{i}}}\right)$. In addition to thermal noise, surface clutter can potentially limit the radar detection sensitivity in the lowest layers of the atmosphere. In this case, we find that the large pulse bandwidth of 6 MHz and pulse duration that is well below the decoherence time result in a surface clutter profile that is largely confined to the lowest few hundred meters above the surface, with the surface equivalent reflectivity falling below −40 dBZ at a height of 200 m. Therefore, we do not expect this clutter source to limit detection sensitivity. We note that this contrasts with the experience of RainCube, where measurement sensitivity is degraded in the lowest 0.5–1 km of the atmosphere because of surface clutter (Battaglia et al., 2020). However, key differences between the two systems that explain this difference are the relative chirp bandwidths and the ratio of pulse duration to coherence time, which in the case of RainCube is larger than 1.

2.3 DAR measurement methodology

The physical basis of the DAR measurement has been discussed in a number of previous works. Here, we focus on the complete treatment first detailed in Battaglia and Kollias (2019) and subsequently in Roy et al. (2020) that quantifies the impact of frequency-dependent hydrometeor scattering and attenuation on the DAR humidity retrieval, and we simply recount the main results. First, we define the local, observed extinction coefficient,

$$\begin{array}{}\text{(4)}& {\mathit{\beta }}_{\mathrm{obs}}(r,f)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial }{\partial r}}\ln \left(Z_{\mathrm{obs}}(r,f)\right).\end{array}$$

Anticipating that differential extinction will be measured over a baseline of N_b range bins, we form the finite difference version of Eq. (4),

$$\begin{array}{}\text{(5)}& {\mathit{\beta }}_{\mathrm{obs}}(r,f)\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{N}_{\mathrm{b}}\mathrm{\Delta }r}}\ln \left({\displaystyle \frac{Z_{\mathrm{obs}}\left(r^{-},f\right)}{Z_{\mathrm{obs}}\left(r^{+},f\right)}}\right),\end{array}$$

where $r^{\pm }=r\pm N_{\mathrm{b}}\mathrm{\Delta }r/\mathrm{2}$. For small changes in range N_bΔr and assuming negligible multiple scattering, the appropriate fitting function for β_obs can be shown to be (Battaglia and Kollias, 2019; Roy et al., 2020)

$$\begin{array}{}\text{(6)}& {\widehat{\mathit{\beta }}}_{\mathrm{obs}}(r,f_i)=a_{\mathrm{1}}+a_{\mathrm{2}}(f_i-f_{\mathrm{0}})+a_{\mathrm{3}}{\mathit{\kappa }}_{\mathrm{v}}(r,f_i)+{\mathit{\beta }}_{\mathrm{g},\mathrm{d}}(r,f_i),\end{array}$$

where we assume that measurements are made at a discrete set of frequencies {f_i}, the values a_j are the regression coefficients, κ_v(r,f) is the water vapor absorption cross section per unit mass, and ${\mathit{\beta }}_{\mathrm{g},\mathrm{d}}(r,f)$ is the absorption coefficient for dry air. Defining the unattenuated reflectivity $Z(r,f)=Z_{\mathrm{obs}}(r,f)\exp \left(\mathrm{2}\mathit{\tau }\right(r,f\left)\right)$, where τ is the one-way optical depth, the differential backscatter function $\mathit{\alpha }(r,f)=\left(\mathrm{2}{N}_{\mathrm{b}}\mathrm{\Delta }r{)}^{-\mathrm{1}}\ln \right[Z(r^{-},f)/Z(r^{+},f)]$, and the hydrometeor extinction coefficient β_h(r,f), the regression coefficients in Eq. (6) can be shown to relate to the following physical quantities:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{a}_{\mathrm{1}}\leftrightarrow \mathit{\alpha }(r,f_{\mathrm{0}})+{\mathit{\beta }}_{\mathrm{h}}(r,f_{\mathrm{0}})\\ {\displaystyle }& {\displaystyle }{a}_{\mathrm{2}}\leftrightarrow {\displaystyle \frac{\partial }{\partial f}}{\left.\left(\mathit{\alpha }(r,f)+{\mathit{\beta }}_{\mathrm{h}}(r,f)\right)\right|}_{f=f_{\mathrm{0}}}\\ \text{(7)}& {\displaystyle }& {\displaystyle }{a}_{\mathrm{3}}\leftrightarrow {\mathit{\rho }}_{\mathrm{v}}\left(r\right).\end{array}$$

In these relations, it is understood that all functions of r are approximated by their average value between r⁻ and r⁺.

With three regression parameters, the model in Eq. (6) requires observations at three or more frequencies in order to perform a least-squares fit. For N_f=3, the solution is independent of the measurement covariance matrix and can be written as $\widehat{\mathit{x}}=K^{-\mathrm{1}}(\mathit{y}-\mathit{b})$, where $\widehat{\mathit{x}}=[{\widehat{a}}_{\mathrm{1}},{\widehat{a}}_{\mathrm{2}},{\widehat{a}}_{\mathrm{2}}{]}^{\mathrm{T}}$, K_i0=1, $K_{i\mathrm{1}}=f_i-f_{\mathrm{0}}$, $K_{i\mathrm{2}}={\mathit{\kappa }}_{\mathrm{v}}(r,f_i)$, $b_i={\mathit{\beta }}_{\mathrm{g},\mathrm{d}}(r,f_i)$, and $y_i={\mathit{\beta }}_{\mathrm{obs}}(r,f_i)$. It is instructive to use an example LES hydrometeor profile to compare humidity profiles estimated using this three-frequency model and those derived using the two-frequency differential absorption approach (e.g., Lebsock et al., 2015; Millán et al., 2016; Roy et al., 2018). Specifically, in the two-frequency approach, one assumes that the contributions from a₂ and β_g,d to observed differential extinction in Eq. (6) are negligible and therefore can estimate the humidity as ${\widehat{\mathit{\rho }}}_{\mathrm{v}}\left(r\right)=\left[{\mathit{\beta }}_{\mathrm{obs}}\right(r,f_{\mathrm{1}})-{\mathit{\beta }}_{\mathrm{obs}}(r,f_{\mathrm{0}}\left)\right]\left[{\mathit{\kappa }}_{\mathrm{v}}\right(r,f_{\mathrm{1}})-{\mathit{\kappa }}_{\mathrm{v}}(r,f_{\mathrm{0}}){]}^{-\mathrm{1}}$. Figure 1 compares the humidity profiles retrieved using these two approaches for two different atmospheric profiles from the RICO LES case with different precipitation characteristics. To highlight the fact that the two-frequency DAR measurement can exhibit biased retrievals even for a perfect radar measurement, here we use reflectivity profiles calculated within the single-scattering approximation at the LES model resolution (i.e., not averaged over the radar footprint). These results show the added utility of a third DAR frequency channel over a wide transmitter bandwidth to eliminate hydrometeor-scattering-induced biases in the humidity retrieval.

Figure 1 Comparison of two-frequency and three-frequency DAR humidity retrievals (see Sect. 2.3) using atmospheric profiles from the RICO LES. An ideal reflectivity measurement is assumed, with no multiple scattering, arbitrarily high range resolution and measurement sensitivity, and zero random measurement error. Retrievals are performed using a differential absorption baseline equal to the model vertical resolution of $r^{+}-r^{-}=\mathrm{40}$ m. (a) Humidity retrievals vs. LES truth (left subpanel) for a lightly precipitating grid cell; (right subpanel) Liquid water mixing ratio profiles and corresponding observed reflectivity at f₁=168 GHz. (b) Same as panel (a) but for a moderate precipitation case. In both cases, performing the three-parameter fit using the model function in Eq. (6) removes the humidity biases that originate from frequency dependence of hydrometeor backscatter and extinction.

2.4 Water vapor retrieval algorithm

The retrieval algorithm employed in this work follows a similar formalism to that developed in Roy et al. (2020), but with some improvements and important modifications to accommodate the N_f≥3 DAR approach and additional regression parameter at each retrieval height corresponding to the a₂ term that is linear in frequency in Eq. (6). In essence, the task is to relate the radar observations y to a forward model F(x,c) that is a generalization of the fitting function in Eq. (6) to model the entire vertical profile. Here x is the state vector consisting of quantities that will be retrieved as part of the inverse problem and c contains parameters that must be assumed to compute the forward model. We note at the outset that the inverse problem will be formulated in such a way that the forward model will be a linear transformation of the state vector, or $\mathit{F}=K\mathit{x}+\mathit{b}$, where b is the part of the forward model that does not depend on the state vector, and it is understood that K and b depend on the forward model assumptions c.

We define the observation vector as $\mathit{y}=[{\mathit{y}}_{\mathrm{0}},{\mathit{y}}_{\mathrm{1}},{\mathit{y}}_{\mathrm{2}}]$, where $[{\mathit{y}}_j{]}_i=\ln (Z_{\mathrm{obs}}(r_i,f_j)Z_{\mathrm{0}}^{-\mathrm{1}})$ for $\mathrm{0}\le i\le M-\mathrm{2}$, $[{\mathit{y}}_j{]}_{M-\mathrm{1}}=\ln ({\mathit{\sigma }}_{\mathrm{obs}}^{\mathrm{0}}\left(f_j\right))$, and $\mathit{r}=[r_{\mathrm{0}},r_{\mathrm{1}},\mathrm{\dots },r_{M-\mathrm{1}}]$ contains the range-bin positions for which the radar signals exceed the instrument sensitivity threshold at all three frequencies. Given the full range vector ${\mathit{r}}^{\prime }=[r_{\mathrm{TOA}},r_{\mathrm{TOA}}+\mathrm{\Delta }r,\mathrm{\dots },r_{\mathrm{s}}-\mathrm{\Delta }r,r_{\mathrm{s}}]$ of length M^′, where r_TOA and r_s are the top-of-atmosphere and surface ranges, respectively, we define the measurement projection matrix P_r according to $\mathit{r}=P_{\mathrm{r}}{\mathit{r}}^{\prime }$. In cases of severe hydrometeor beam attenuation, the surface return may be below the sensitivity threshold ${\mathit{\sigma }}_{\min }^{\mathrm{0}}$, in which case the final element of each y_j is replaced by the last cloud bin with dBZ_obs>dBZ_min. The measurement covariance matrix is defined in block-diagonal form as $S_y=\mathrm{diag}(S_{y_{\mathrm{0}}},S_{y_{\mathrm{1}}},S_{y_{\mathrm{2}}})$, where $[S_{y_j}{]}_{ik}=({\mathit{\sigma }}_Z{Z}_{\mathrm{obs}}^{-\mathrm{1}}{|}_{r_i,f_j}{)}^{\mathrm{2}}{\mathit{\delta }}_{ik}$, and the relative error ${\mathit{\sigma }}_Z{Z}_{\mathrm{obs}}^{-\mathrm{1}}$ is evaluated at range r_i and frequency f_j using Eq. (3). Note that in general y has dimension N_f×M and that this formalism is easily extended to more than three frequencies.

Next, we define the state vector according to $\mathit{x}=[{\mathit{x}}_{\mathrm{0}},{\mathit{x}}_{\mathrm{1}},{\mathit{x}}_{\mathrm{2}}]$, where $[{\mathit{x}}_{\mathrm{0}}{]}_i=\ln (Z_r(r_i,f_{\mathrm{0}})Z_{\mathrm{0}}^{-\mathrm{1}})$ is a vector of dimension M, Z_r(r,f₀) is the retrieved reflectivity at f₀ including attenuation from hydrometeors, ${\mathit{x}}_{\mathrm{1}}=[{\mathit{\gamma }}_{\mathrm{0}},{\mathit{\gamma }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\gamma }}_{M-\mathrm{1}}]$, where γ_i is a range-dependent slope vs. frequency and is analogous to a₂ in Eq. (6), and x₂ contains the N retrieved humidity values and requires careful specification given an arbitrary measurement projection P_r. The forward model is then defined in block form as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathit{F}(\mathit{x},\mathit{c})& =\left(\begin{array}{ccc}{I}_M& (f_{\mathrm{0}}-f_{\mathrm{0}})I_M& -\mathrm{2}{T}_{\mathrm{0}}\\ I_M& (f_{\mathrm{1}}-f_{\mathrm{0}})I_M& -\mathrm{2}{T}_{\mathrm{1}}\\ I_M& (f_{\mathrm{2}}-f_{\mathrm{0}})I_M& -\mathrm{2}{T}_{\mathrm{2}}\end{array}\right)\left(\begin{array}{c}{\mathit{x}}_{\mathrm{0}}\\ {\mathit{x}}_{\mathrm{1}}\\ {\mathit{x}}_{\mathrm{2}}\end{array}\right)\\ & -\mathrm{2}\left(\begin{array}{c}{\mathit{\tau }}_{\mathrm{g},\mathrm{d}}\left(f_{\mathrm{0}}\right)\\ {\mathit{\tau }}_{\mathrm{g},\mathrm{d}}\left(f_{\mathrm{1}}\right)\\ {\mathit{\tau }}_{\mathrm{g},\mathrm{d}}\left(f_{\mathrm{2}}\right)\end{array}\right),\end{array}\end{array}$$

where I_M is the M×M identity matrix, T_i is a matrix of size M×N that computes the water vapor absorption contribution to the total optical depth profile given the values in x₂, and τ_g,d(f_i) is the optical depth profile at f_i due to dry air absorption. At this point, we must assume vertical profiles of temperature T(r^′) and pressure P(r^′) at the radar range resolution, which are elements of c, in order to calculate the gas absorption quantities τ_g,d(f_i) and T_i. For this work, we do so by averaging these thermodynamic profiles over the whole LES domain and smoothing the result by performing a convolution with a box of length 2 km. Note that due to the relatively weak dependence of the water vapor absorption cross section on pressure and temperature, the systematic error associated with these assumptions is not a leading order effect (Roy et al., 2018). Then, the dry air optical depth vector is easily defined by the recursive relation $\left[{\mathit{\tau }}_{\mathrm{g},\mathrm{d}}\right(f_j){]}_i=[{\mathit{\tau }}_{\mathrm{g},\mathrm{d}}\left(f_j\right){]}_{i-\mathrm{1}}+\mathrm{\Delta }r{\mathit{\beta }}_{\mathrm{g},\mathrm{d}}(r_{i-\mathrm{1}}^{\prime },f_j)$ and the initial condition $\left[{\mathit{\tau }}_{\mathrm{g},\mathrm{d}}\right(f_j){]}_{\mathrm{0}}=\mathrm{0}$.

The first step in defining T_i, and therefore the structure and interpretation of x₂, is to prescribe a maximum vertical resolution for the retrieved humidity profile, which we call Δz. Furthermore, we restrict Δz to values that satisfy the condition that $\mathrm{\Delta }z/\mathrm{\Delta }r=O$ is an integer, and we define the new vertical axis ${\mathit{z}}^{\prime }=[z_{\mathrm{TOA}},z_{\mathrm{TOA}}-\mathrm{\Delta }z,\mathrm{\dots },\mathrm{\Delta }z+\mathrm{\Delta }r,\mathrm{\Delta }r]$, where the lowest height position corresponds to the first range bin above the surface r_s−Δr, and z_TOA corresponds to the range $r^{\prime }=r_{\mathrm{TOA}}+\mathrm{\Delta }z-\mathrm{\Delta }r$. In this work, we use Δz = 200 m, resulting in a range-bin oversampling factor of O=4. To determine which vertical humidity values will be retrieved given a radar measurement projection P_r, we define a new projection matrix P_z and axis z, with $\mathit{z}=P_{\mathrm{z}}{\mathit{z}}^{\prime }$, and P_z projects onto the non-zero elements of the following downsampling vector d:

$$\begin{array}{}\text{(9)}& [\mathit{d}{]}_i=\sum _{j=\mathrm{0}}^{M-\mathrm{1}}\sum _{k=Oi}^{Oi+O-\mathrm{1}}[P_{\mathrm{r}}{]}_{jk}.\end{array}$$

Note that this definition for z ensures that a humidity value at a certain height z will only be retrieved if there is at least one radar observation within the interval $\left[z-\mathrm{\Delta }r,z+(O-\mathrm{1})\mathrm{\Delta }r\right)$.

Next, we write the matrix T_i as the product of three matrices, $T_i=P_{\mathrm{r}}{T}_i^{\prime }A$, where A ($M^{\prime }\times N$) first interpolates the humidity state vector x₂ to the r^′ full range space, then $T_i^{\prime }$ ($M^{\prime }\times M^{\prime }$) computes the full optical depth profile, and finally P_r projects this profile onto the measurement range vector r. In contrast with the humidity interpolation approach utilized in Roy et al. (2020), here we interpolate to the Δr resolution between retrieval heights in z using an exponential profile with fixed scale height H_ρ. For situations where sequential elements for z are closely spaced, the choice of interpolation function is irrelevant, since the absorption line shape properties change negligibly and thus the retrieval is only sensitive to the integrated water vapor (IWV) between the two retrieval heights. However, for situations with widely spaced elements of z, including total column water vapor (TCWV) retrieval scenarios where no clouds are detected, the exponentially interpolated humidity profile provides a realistic weighting of the vertical absorption cross section. We note that the assumption of the exponential profile is not necessary. We use this assumption here to demonstrate how a radar-only retrieval might perform in the absence of any additional information. In practice, in an operational environment a more detailed profile shape would be taken from a weather analysis system or from coincident observations from passive sounders. In fact, we show below that the derived column water vapor measurements below are significantly sensitive to the exponential scale height motivating the future use of ancillary data in the DAR retrievals. Finally, the matrix $T_i^{\prime }$ can be shown to be a lower triangular matrix defined by

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& T_i^{\prime }=\\ & \mathrm{\Delta }r\left(\begin{array}{llllll}{\mathit{\kappa }}_{\mathrm{v}}(r_{\mathrm{0}}^{\prime },f_i)& \mathrm{0}& \mathrm{\cdots }& & & \mathrm{0}\\ {\mathit{\kappa }}_{\mathrm{v}}(r_{\mathrm{0}}^{\prime },f_i)& {\mathit{\kappa }}_{\mathrm{v}}(r_{\mathrm{1}}^{\prime },f_i)& \mathrm{0}& \mathrm{\cdots }& & \mathrm{⋮}\\ {\mathit{\kappa }}_{\mathrm{v}}(r_{\mathrm{0}}^{\prime },f_i)& {\mathit{\kappa }}_{\mathrm{v}}(r_{\mathrm{1}}^{\prime },f_i)& {\mathit{\kappa }}_{\mathrm{v}}(r_{\mathrm{2}}^{\prime },f_i)& \mathrm{0}& \mathrm{\cdots }& \mathrm{⋮}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{⋮}\\ & & & & \mathrm{\ddots }& \mathrm{0}\\ {\mathit{\kappa }}_{\mathrm{v}}(r_{\mathrm{0}}^{\prime },f_i)& {\mathit{\kappa }}_{\mathrm{v}}(r_{\mathrm{1}}^{\prime },f_i)& {\mathit{\kappa }}_{\mathrm{v}}(r_{\mathrm{2}}^{\prime },f_i)& \mathrm{\cdots }& \mathrm{\cdots }& {\mathit{\kappa }}_{\mathrm{v}}(r_{M^{\prime }-\mathrm{1}}^{\prime },f_i)\end{array}\right).\end{array}\end{array}$$

An example of the humidity profile grid selection and interpolation procedure is given in Appendix C.

To proceed towards retrieving the best estimate of the state vector x, we define the cost function to equal the weighted least-squares sum $C\left(\mathit{x}\right)=(\mathit{y}-K\mathit{x}-\mathit{b}{)}^{\mathrm{T}}{S}_y^{-\mathrm{1}}(\mathit{y}-K\mathit{x}-\mathit{b})$. Note that, because of the introduction of the model parameters in x₁ corresponding to the linear frequency fit term, the DAR retrieval no longer requires regularization as was implemented in Roy et al. (2020), and there is no need to introduce a systematic error covariance matrix, because the current model is designed to eliminate the hydrometeor-scattering-induced bias. Furthermore, it worth highlighting here that the only a priori information necessary for this retrieval is specification of representative T and P profiles, which exemplifies the sharp contrast between this inversion approach and that of optimal estimation. The value of x that minimizes C(x), as well as the estimated state's covariance matrix, is then calculated as follows:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\widehat{S}}_x& =(K^{\mathrm{T}}{S}_y^{-\mathrm{1}}K{)}^{-\mathrm{1}}\\ & \widehat{\mathit{x}}={\widehat{S}}_x{K}^{\mathrm{T}}{S}_y^{-\mathrm{1}}(\mathit{y}-\mathit{b}).\end{array}\end{array}$$

For much of the analysis in Sect. 3, it is more appropriate to analyze the IWV for individual partial columns of the retrieval (i.e., corresponding to a single column of the matrix A as shown in Fig. C1b) than it is to analyze water vapor densities at specific heights, as this provides a common framework to discuss both in-cloud profiling retrievals with high vertical resolution and partial and total column retrievals encompassing large cloud-free volumes. It is easy to show that converting the estimated humidity state vector x₂ to a partial-column IWV variable amounts to a simple linear transformation c=Lx₂, where $L_{ii}=\mathrm{\Delta }r{\sum }_j{A}_{ji}$ is a diagonal matrix where each entry is the sum of all elements in a column of A. The corresponding IWV covariance matrix is then calculated according to $S_c=L{S}_{x_{\mathrm{2}}}{L}^{\mathrm{T}}$, where $S_{x_{\mathrm{2}}}$ is the square block from ${\widehat{S}}_x$ that pertains to the retrieved humidity values x₂.

3 Results

In this section we detail the results of case-by-case forward simulations and retrievals for the five LES scenarios described in Sect. 2.1. Because the analysis framework is identical for all five cases, we provide the most detail in Sect. 3.1 in describing the GATE simulations, with more limited discussion in the sections for the remaining four cases.

Before presenting the results, however, it is important to decide on a framework for conveying the independent aspects of (1) random uncertainty, which stems from reflectivity measurement error, and (2) systematic uncertainty, or bias. Under the assumption of small relative measurement error ${\mathit{\sigma }}_Z{Z}_{\mathrm{obs}}^{-\mathrm{1}}\ll \mathrm{1}$, which is ensured by our specification of the minimum detectable reflectivity, the errors in the observation space (y) are Gaussian distributed. Furthermore, since the forward model is linear in the state vector, this Gaussianity is preserved in going from observation to state space (x). Therefore, the estimators for the state vector and its covariance matrix in Eq. (11) describe a multivariate normal distribution from which a real observing system would sample. Note that because our inverse algorithm does not impose a prior as is done in the optimal estimation approach, the estimated state vector and covariance matrix are not biased by the assumed prior distribution. One approach often used in instrument simulator studies to assess the effect of random error on the inferred atmospheric state and its uncertainty is to inject random noise into the observations according to an error model like that in Eq. (3). However, this approach would be redundant in this case because of the detailed uncertainty quantification given by ${\widehat{S}}_x$ and the Gaussianity properties described above. Furthermore, by not adding random noise to the observed reflectivity profiles, we can utilize the humidity retrieval mean values in $\widehat{\mathit{x}}$ to assess measurement bias on a grid-cell-by-grid-cell basis.

For these reasons, we proceed by simulating observations without the injection of random noise and subsequently analyzing the simulated retrievals by separately addressing the issues of systematic and random error. While the variability of humidity at a given height within the LES is small, the retrieval estimate of that quantity in the state vector $\widehat{\mathit{x}}$ can vary widely from grid cell to grid cell due to changing hydrometeor fields. Therefore, by analyzing the variability of retrieved state mean values, we can assess the level of systematic error in the DAR measurement due to effects including non-uniform beam-filling, multiple scattering, incorrect humidity interpolation assumptions, and insufficiency of the linear-in-frequency fit term for mitigating the hydrometeor scattering biases. Then, we separately examine the retrieval random uncertainty and its scaling with along-track averaging distance by combining sequential along-track grid cells using a weighted mean and variance approach. Furthermore, to explore the systematic error that comes from assuming an exponential humidity interpolation function, we perform retrievals for each case using two different scale heights of H_ρ=1.5 and 2.5 km.

Finally, it is important to recognize that while we can calculate from the LES a 3D snapshot of observed reflectivity values, the notional DAR system is not intended to be able to measure such a 3D field in a single along-track overpass. Instead, we treat the simulated 3D reflectivity field as providing a rich set of statistically independent vertical profiles that can be used to assess the accuracy and precision of measurements from a DAR that would sample similar clouds over an along-track distance that is typically far greater than the linear size of the LES domain. To this end, we will examine two different beam sampling strategies in the following analysis: one using a fixed, nadir-pointing beam that executes sequential along-track overpasses in a raster-scan fashion across the LES domain; and a second approach where we assume that the radar possesses intelligent pointing capabilities and can arbitrarily target a single pixel in the cross-track (x) dimension for each along-track (y) position. We envision a cross-track scanning capability on the order ±3^∘, which from 400 km altitude provides a 42 km range in cross-track sampling location. This relatively limited scanning capability would provide significant increases to cloud sampling opportunities in shallow convective cloud regimes while also limiting the technical complexity of a wide-angle scan. We note that we do not model the angular dependence of the NRCS in this study, which has little variation over this modest angular range (Roy et al., 2021). For this latter approach, we envision that the spacecraft would be outfitted with an additional passive sensor, for instance a microwave imager, that can provide an estimate of cloud liquid water path in advance of the radar measurement and can therefore feed-forward this information to the intelligently scanning antenna. In this work, we assume that for each along-track pixel the antenna scans to the cross-track location with maximum liquid water path. While such a capability is highly idealized, it provides a useful benchmark for assessing the limits of an intelligently scanned system.

3.1 Deep convection: GATE

Forward simulated DAR observations for the GATE deep convection case are presented in Fig. 2. In general for the different LES cases, results are shown for the online frequency f₂=174.8 GHz because it experiences the most attenuation and therefore is the limiting factor for the three-channel DAR system in terms of sampling. For instance, as shown in the vertical reflectivity profiles for all three frequencies in Fig. 2d, the two offline channels f₀ and f₁ penetrate the entire cloud column down to the surface, while the online channel is fully extinguished at a height of about 2 km. The observed ocean surface NRCS in Fig. 2a reveals the horizontal structure of this deep convection scene, showing how regions near plumes contain a significant hydrometeor burden that prohibits detection of surface echoes at f₂. Nevertheless, there is still substantial coverage in both the horizontal and vertical (see Fig. 2c) near these deep convective clouds. Finally, Fig. 2d provides a clear picture of the water vapor differential absorption signature in the observed multi-frequency reflectivity profiles with continuous sampling at all frequencies from roughly 2 to 12 km in height.

Figure 2 Selected results from GATE forward simulations. (a) Observed surface NRCS at f₂=174.8 GHz. Surface echoes with SNR below the detection threshold appear as white and generally correlate with convectively active regions with large liquid and ice water content (LWC and IWC) values. (b) Two-dimensional cross section of model condensed water profiles (ice and liquid) for the along-track path (dashed line) in panel (a). (c) Corresponding observed reflectivity profiles at 174.8 GHz. Gray regions in lower-tropospheric convective cores indicate regions where the multiply scattered signal is at least 3 dB higher than the single-scattering reflectivity. We assume that these grid points can be filtered out before performing water vapor retrievals. (d) Vertical reflectivity profiles at all three DAR frequencies corresponding to the dashed line in panel (c).

The results of the GATE retrievals are summarized in Fig. 3, where panels a and b pertain only to systematic measurement error, and panel c shows the scaling of random uncertainty, or precision, with along-track averaging for the two radar sampling approaches, fixed nadir and intelligently scanned. We note that the GATE case is the only one for which multiple scattering features noticeably in the forward simulations, and we assume that range bins for which the multiply scattered signal is 3 dB higher than the single-scattering reflectivity can be identified and filtered out before retrieval. The retrieval humidity values x₂ are classified according to three different types of DAR measurement columns: the “top column” between the radar and the first detected range bin; the “in-cloud” columns that include cloud returns at both the near and far ranges; and the “cloud-to-surface” column between the final cloud bin and the surface. In Fig. 3a, the DAR retrieval (blue circles) corresponds to the average at each height of the in-cloud profiling retrievals that have a vertical resolution of 200 m. That is, partial-column retrievals that are still in cloud but have a column lengths greater than 200 m are excluded as these do not represent high-resolution, local humidity measurements. Importantly, the blue error bars in Fig. 3a do not represent random measurement error but instead are the standard deviations of the mean retrieval values at each height for the entire LES domain. Therefore, these error bars provide a measure of the spread of in-cloud profiling biases, while the close agreement between the mean values and the conditionally sampled truth profile (open black circles) reveals that the mean bias is small.

Figure 3 GATE DAR retrieval results for systematic (a, b) and random (c) error assessments. (a) Lower-tropospheric mean-retrieved water vapor profile (blue solid circles) compared with model values, including the full domain mean profile (solid black line), range of model humidity at each height (gray shaded region), and the conditionally sampled truth profile (black circles), which is calculated using only grid cells for which the DAR retrieved a 200 m resolution profile value. Note that the DAR retrieval error bars correspond to the standard deviation of mean-retrieved values (i.e., elements of $\widehat{\mathit{x}}$) at each height and do not depict measurement uncertainty from random error. (b) Two-dimensional histograms of IWV bias from retrievals using H_ρ=2.5 km as a function of truth IWV, for the three distinct DAR column retrieval types (top: lowest column between surface and last cloud return; middle: partial columns within/between cloud layers; bottom: uppermost column between top-of-atmosphere and first detected range bin). Each vertical histogram is normalized to have a peak value of 1 for display purposes. Solid black lines correspond to the median bias when using H_ρ=2.5 km, while dashed lines are the same for H_ρ=1.5 km. (c) Retrieval precision scaling with along-track averaging, where the in-cloud curves (top panel) correspond to a profiling height of 3.9 km. For each LES case, this profiling height is chosen to lie well within the average cloud and precipitation volume to avoid sparse sampling near the lower and upper boundaries. Dashed lines correspond to an intelligently scanned radar (see text). Solid circles represent the humidity precision scaling for a hypothetical along-track path that pieces together nadir-pointing transects for x=21.1, 42.3, 63.1, and 84.3 km.

The results of the IWV bias assessment for all three retrieval column types are shown in Fig. 3b. Here the histograms correspond to retrievals performed using a humidity scale height of H_ρ=2.5 km, with the median biases shown by the solid black lines. The black dashed lines show the median IWV biases that result from retrievals performed using the other scale height of H_ρ=1.5 km. For this analysis, we include all columns of the in-cloud type and do not restrict them to be of 200 m length. The magnitude of the median bias never exceeds 2 mm IWV for either value of H_ρ, and the results for H_ρ=2.5 km show the exceptional accuracy of this DAR measurement method in capturing both in-cloud and clear-air water vapor columns, with biases rarely exceeding 0.5 mm.

We assess the scaling of in-cloud, 200 m profiling precision with along-track averaging distance (Fig. 3c) by performing a weighted average of retrievals at a specific height of 3.9 km using the both fixed-nadir (orange circles) and intelligently scanned (orange dashed line) approaches. The result is that after 100 km of along-track averaging, these sampling approaches yield random humidity uncertainties of 2 and 1 g m⁻³, respectively. It is important to note that the humidity profiling resolution and random uncertainty can be traded off in a linear fashion, meaning we would expect profiles with 400 m resolution to exhibit half the uncertainty shown here. Because the cloud-to-surface and top columns are typically measured over a much larger differential absorption baseline, the amount of averaging necessary to reach a specified relative uncertainty is reduced relative to the in-cloud profiling. Here, retrieval precisions for both column types are at the 1 mm level after only 10 km of along-track averaging. These column IWV retrieval capabilities could provide unprecedented observations of horizontal water vapor variability in deep convective cold pools.

3.2 Precipitating shallow cumulus: RICO

Forward simulated DAR observations for the RICO case with precipitating shallow cumulus clouds are presented in Fig. 4. Because of the high horizontal resolution of the LES relative to the radar footprint, we compare the simulated observed ocean surface NRCS both before (panel a) and after (panel b) averaging over the instrument beam pattern. In this case, the reduced hydrometeor burden relative to the deep convection regions of the GATE simulation results in nearly ubiquitous sampling of the surface at f₂. In Fig. 4c and d, we show the model cloud and rain liquid water content (LWC) profiles and corresponding observed reflectivity for the fixed-nadir path shown by the dashed line in panels a and b. Even at the highly absorbing online frequency f₂, the radar is able to sample almost a large fraction of the LWC field within the PBL. Furthermore, we highlight using gray shading in Fig. 4d the along-track positions for which the cloud-to-surface column will be retrieved.

Figure 4 Selected results from RICO forward simulations. Simulated observed ocean surface NRCS at 174.8 GHz before (a) and after (b) filtering for instrument horizontal resolution and detection sensitivity. (c) LWC profiles for the along-track transect indicated by the dashed line in panels (a) and (b). (d) Corresponding observed reflectivity profiles at 174.8 GHz. Positions where cloud-to-surface partial-column water vapor retrievals can be performed are highlighted with gray shading.

The retrieval results for the RICO case are presented in Fig. 5 in a nearly identical format to the GATE simulations of the previous section. Because of the reduced sampling of IWV bias vs. truth relative to GATE, which is a result of the decreased cloud and precipitation amount, we present the bias results as scatter plots instead of histograms for both the H_ρ=1.5 km (Fig. 5b) and H_ρ=2.5 km (Fig. 5c) cases. An important feature of the in-cloud profiling retrievals to highlight is the ability of the DAR to resolve the inversion layer thickness and capture the water vapor profile near the top of the PBL, where there is a rapid change with height of the hydrometeor fields and therefore strong potential for systematic error. In fact, in separate end-to-end simulations run using a radar range resolution of 200 m (not shown), we find a systematic, large negative bias in this regime, which is due to non-uniform hydrometeor filling of range bins near the cloud top. Because the retrieval algorithm interprets reflectivity signals as representing an average throughout the bin, this non-uniformity tends to diminish the differential absorption between the first two cloud bins, causing the negative humidity bias. However, by oversampling the humidity retrieval resolution by a factor of O=4 here (see Sect. 2.4), this non-uniform range-bin-filling issue is mitigated, and biases are significantly reduced.

Figure 5 RICO DAR retrieval results. (a) Comparison of mean-retrieved and model in-cloud humidity profiles, depicted as in the GATE case (Fig. 3a). (b–c) IWV retrieval bias vs. truth for the three distinct column types and two different humidity retrieval scale heights of 1.5 and 2.5 km. (d–e) Retrieval precision scaling with along-track averaging, similar to Fig. 3c, except sequential along-track transects are executed by cycling through all possible cross-track (x) positions.

By comparing the retrieval biases for H_ρ=1.5 and 2.5 km in Fig. 5b and c, we see that the humidity scale height choice mostly affects the top-column retrievals. This is because the in-cloud and cloud-to-surface columns are of such short length that the choice of humidity interpolation function is unimportant, and the retrieval is only sensitive to bulk IWV in the partial columns. Furthermore, even if one uses an inappropriate retrieval scale height of 2.5 km, the top-column biases remain fairly small, with a clear-sky TCWV bias of at most 2 mm, or about 5 % (see the rightmost cluster of top-column retrievals in Fig. 5c).

The retrieval precision analysis in Fig. 5d and e reveals similar performance as in GATE, where in this case we present in-cloud precision results at a profiling height of 1.3 km. However, because in RICO the clouds are limited to the lowest 2.5 km of the atmosphere, the top-column retrievals are always associated with a column that is nearly the depth of the atmosphere, resulting in reduced uncertainty relative to GATE. In this case, it is only necessary to average in the along-track direction for 1 km to achieve a top-column precision of 1 mm IWV and for 20–30 km to realize the same precision for the cloud-to-surface column, depending on the scanning strategy.

3.3 Non-precipitating shallow cumulus: BOMEX

The BOMEX case of non-precipitating shallow cumulus clouds is a particularly challenging one for in-cloud DAR observations, and for cloud radar in general, due to the relatively low volume cloud fraction and therefore limited sampling. This fact is highlighted by the forward simulation results in Fig. 6, where the 2D map of cloud liquid water path (LWP) in panel a suggests that the radar will primarily encounter clear-air columns down to the ocean surface, especially without intelligent scanning capabilities. In Fig. 7c and d, we show the cloud LWC profiles and observed reflectivities at 155.5 and 174.8 GHz for the dashed line path in panel a.

Figure 6 Selected results from BOMEX forward simulations. (a) Liquid water path (LWP) map for the shallow cumulus scene highlighting the low cloud area fraction. (b) Simulated observed ocean surface NRCS at 174.8 GHz at the radar resolution. (c) Cloud LWC profiles corresponding to the dashed line in panel (a). (d) Observed reflectivity profiles at 155.5 (left) and 174.8 (right) GHz.

Figure 7 BOMEX DAR retrieval results. (a) Comparison of mean-retrieved and model in-cloud humidity profiles, depicted as in the GATE case (Fig. 3a). (b–c) Retrieved IWV bias vs. truth for assumed humidity scale heights of 1.5 and 2.5 km. (d–e) Retrieval precision scaling with along-track averaging, where the in-cloud results correspond to a retrieval height of 1.5 km.

The BOMEX retrieval results are summarized in Fig. 7 in a format identical to that for RICO. Unsurprisingly, the limited in-cloud sampling in this case results in increased variance of the profiling bias shown in panel a. Furthermore, the in-cloud profiling precision is diminished, even for the intelligently scanned system, which reaches a precision of 2 g m⁻³ only after an along-track averaging distance of 300 km. We conclude from these results that any hope of obtaining in-cloud profiles in non-precipitating shallow cumulus requires an intelligent scanning capability. Furthermore, even a scanned measurement is unable to resolve the PBL inversion as can be done in the RICO case due to the weak cloud scattering from within the inversion layer. Nevertheless, the notional DAR can retrieve cloud-to-surface and top-column IWV values with good accuracy and precision, especially for the intelligently scanned system, with biases below 0.5 mm for H_ρ=1.5 km and averaging distances of 30 and 3 km necessary to achieve 1 mm IWV precision, respectively.

3.4 Drizzling stratocumulus: VOCALS

The VOCALS case represents the ubiquitous light precipitation and drizzle that is found in marine stratocumulus, for example in areas of open cellular convection. Because of the increase in area cloud fraction relative to shallow cumulus, and general horizontal homogeneity of the cloud field, both VOCALS and DYCOMS-II simulations will reveal little difference between a fixed-nadir and intelligently scanned radar system. Furthermore, the position of the cloud layer atop a shallow PBL with strong temperature inversion and hydrolapse permits a natural division of the atmospheric column into separate regions for the DAR retrieval. Specifically, the top column between the radar and the stratocumulus cloud top always provides a bulk measure of the amount of water vapor in the dry, free troposphere, while the cloud-to-surface column of variable depth ≤1 km quantifies the IWV in the moist PBL.

The forward simulation results for the VOCALS case are detailed in Fig. 8. In panel a, we show the model LWC at the peak cloud height of 1.25 km, highlighting the high area cloud fraction in this case, with resulting observed reflectivity map at the same height in panel b. Even at the online frequency f₂, all horizontal columns contain a significant cloud detection. Figure 8c and d show profiles of model cloud and drizzle LWC profiles and corresponding observed reflectivity for the along-track path indicated in panel a with a dashed line. Because of the reduced water vapor burden in the free troposphere relative to the previous cases of deep and shallow convection, there is greatly reduced atmospheric attenuation at f₂ permitting increased cloud sampling in the PBL. Note that as with the GATE case, we apply a multiple-scattering mask to the simulated observations, where we assume that range bins with a multiply scattered signal that is 3 dB higher than the single-scattering value can be filtered out. Regions where this mask has been applied are indicated with gray shading in Fig. 8d and are likely the result of multiple scattering from drizzle droplets within the cloud layer.

Figure 8 Selected results from VOCALS forward simulations. (a) Model LWC field for z=1250 m. (b) Corresponding observed reflectivity map at 174.8 GHz for the range bin centered at 1250 m height. (c) Vertical LWC profiles for the along-track path indicated by the dashed line in panel (a), with drizzle clearly visible below cloud base, and (d) corresponding 174.8 GHz observed reflectivity profiles. The gray shaded region in panel (d) indicates where the multiple-scattering mask has been applied, using the same criteria as in the GATE case.

The VOCALS retrieval results are presented in Fig. 9. While in-cloud profiling is not expected to be a major strength of DAR for stratocumulus clouds, the retrievals reveal quite promising bias and precision capabilities in this regard, with random uncertainty at the 1 g m⁻³ level after 50 km of along-track averaging (see Fig. 9c). Furthermore, the free tropospheric column and cloud-to-surface column results, with or without intelligent scanning, suggest that the notional DAR can provide very useful IWV estimates in this drizzling stratocumulus scene that could, for instance, be used to constrain near-surface humidity.

Figure 9 VOCALS DAR retrieval results. Panels (a), (b), (d), and (e) are as in Figs. 5 and 7, while panel (c) shows retrieval bias histograms for H_ρ=1.5 km (solid lines) and H_ρ=2.5 km (dashed lines). In panels (d) and (e), note the insensitivity of the precision to choice of scanning strategy (circles: fixed nadir, dashed line: intelligently scanned) because of the horizontal uniformity of the cloud field.

3.5 Non-drizzling stratocumulus: DYCOMS-II

The forward simulation results for the final case, DYCOMS-II, involving a non-drizzling stratocumulus cloud layer of thickness ≈200 m are shown in Fig. 10. As for the VOCALS case, we plot the model cloud LWC field in Fig. 10a for the peak cloud height of 820 m and the corresponding radar reflectivity map in Fig. 10b, revealing that the radar detects cloud in each sampled atmospheric column. From the observed reflectivity profiles in Fig. 10d, we can see that each retrieval will feature three distinct water vapor quantities: the free tropospheric column between the radar and cloud top, a single in-cloud humidity measurement of 200 m vertical resolution, and a cloud-to-surface column of length 600 m representing the sub-cloud bulk humidity in the PBL.

Figure 10 Selected results from DYCOMS-II forward simulations. (a) Model cloud LWC field for z=820 m. (b) Corresponding observed reflectivity map at 174.8 GHz for the range bin centered at 800 m height. (c) Vertical LWC profiles for the along-track path indicated by the dashed line in panel (a) (note the reversed color scale with respect to a) and (d) corresponding 174.8 GHz observed reflectivity profiles.

The results of the retrieval are presented in Fig. 11 in an identical format to the VOCALS case, where we again see that the added utility of an intelligently scanned radar is minimal due to the high area cloud fraction. However, due to the reduced SNR of the cloud targets for the non-drizzling stratocumulus case, the in-cloud profiling and sub-cloud IWV precision feature a less favorable scaling with along-track distance. In realistic measurement scenarios the required >100 km averaging distance may suffice due to the typically large horizontal extents of stratocumulus cloud structures.

Figure 11 DYCOMS-II DAR retrieval results. Because of the limited vertical extent of cloud volume, every horizontal grid cell features just a single in-cloud profiling retrieved value and identical upper- and lower-tropospheric column lengths. All panels have identical descriptions to those in Fig. 9.

4 Discussion

In the previous section, we presented the basic outputs of the end-to-end simulations and analyzed them from a statistical perspective with a focus on disentangling systematic and random retrieval uncertainty. Here, we provide example science case studies for utilizing the DAR observations and begin the discussion with a summary of the DAR retrieval capabilities within the context of the subtropical to tropical cloud transition.

4.1 Retrieval performance summary

In order to synthesize a coherent picture of DAR retrieval capabilities for all five LES cases, we generate the following three quantities from the basic retrieval product of partial-column IWV and its associated uncertainty, assuming an intelligently scanned radar system and variable along-track averaging distance. First, we transform the cloud-to-surface IWV retrieval for each along-track position into an estimate of the surface humidity, assuming a well-mixed boundary layer for all cases except GATE, where the partial columns often extend from the anvil cloud bottom to the surface and therefore would not be well captured by a well-mixed assumption. Therefore, for the GATE case we assume an exponential humidity distribution with a scale height of 2.5 km. Second, we again filter the in-cloud retrievals for only those with 200 m vertical resolution. Third, we transform the top-column IWV estimate into a representative estimate of the humidity 500 m above the highest in-cloud profiling position, for which we use the same exponential profile and humidity scale height as was used in the retrieval. Note that it should be expected that this third quantity does not agree with the model humidity value at the same height, because the upper-tropospheric humidity profile shape is typically not well captured by an exponential profile. These three measured quantities and their uncertainties are then combined using a weighted mean and variance approach for sequential along-track grid cells to produce average observed profiles with associated random uncertainty.

Figure 12 shows the results of this analysis, with an overlaid schematic depiction of the marine trade-wind subtropical to tropical cloud transition. The variable along-track averaging distances, indicated by the respective Δy values in the figure, are determined based on the following criteria. For the surface and upper-tropospheric humidity estimates, the averaging distance that produces a 10 % relative measurement error is used, with a minimum distance of 10 km imposed to ensure that potential systematic error is also accounted for. Then, for the in-cloud profiles, we again use the along-track distance that produces a specified relative uncertainty for at least one profiling height but vary this number depending on the case to account for differences in cloud sampling and SNR. The minimum relative uncertainty criteria for the five cases in Fig. 12 are as follows: 15 % (GATE), 10 % (RICO), 20 % (BOMEX), 10 % (VOCALS), and 20 % (DYCOMS-II).

Figure 12 Summary of humidity retrieval capabilities for an intelligently scanned DAR system displayed within the context of the subtropical-to-tropical cloud transition regimes. For each LES case, we present the in-cloud profiles and surface humidity estimates calculated using the cloud-to-surface partial column with error bars representing 1σ random error. In-cloud estimates in the lower troposphere with relative uncertainty greater than 100 % are excluded. The along-track averaging distances for the different retrieval quantity estimates are listed above. The top-column representative water vapor density value (orange circles) is compared to the model truth value (open black circles). While the retrieval biases are nearly identical for the fixed-nadir case with large along-track averaging, the precision is degraded and can be assessed by examining the precision scaling plots in Figs. 3d, 5d, 7d, 9d, and 11d.

These results constitute a concise summary of the extensive simulations for the notional DAR instrument and what its humidity measurement capabilities would be in these various cloud scenarios. It is important to note that both the high-resolution in-cloud profiles and the high-precision sub-cloud humidity observations in Fig. 12 are inaccessible with other sensor platforms, with the possible exception of GNSS radio occultation (RO) techniques. Specifically, GNSS RO systems can provide high-resolution vertical profiles of refractivity in clouds but must partition that refractivity signal between temperature and water vapor. Furthermore, GNSS RO retrievals struggle to profile below strong inversions due to the associated hydrolapse and in general measure an average over clear-air and cloudy regions with 100 km scale horizontal resolution. DAR, on the other hand, is not affected by sharp thermodynamic profile gradients and can utilize the surface echoes to probe water vapor all the way to the surface.

4.2 Constraining near-surface RH

The surface relative humidity (RH) is a critical quantity in determining PBL structure. It sets the lifting condensation level and thus the cloud base and governs the magnitude of the surface moisture fluxes. As an example application of the DAR observations, we examine the ability to constrain the near-surface RH using the surface humidity estimates from Fig. 12. Recall that these estimates themselves are derived from the cloud-to-surface IWV column retrievals by assuming a well-mixed boundary layer. Therefore, two independent potential sources of bias in estimating near-surface RH are the retrieval IWV biases for the cloud-to-surface column and deviations of the atmospheric profile from the well-mixed state. The results are presented in Table 3, where we use the DAR-inferred surface humidity ρ_v,DAR and assume perfect knowledge of the surface temperature from the model to derive the DAR RH estimate. Then, the sensitivity of this RH estimate to the surface temperature assumption is quantified by computing $\partial \mathrm{RH}/\partial T$. While there is clearly RH sensitivity in the DAR cloud-to-surface column measurement, the relatively large uncertainty of ${\mathit{\sigma }}_{\mathrm{v},\mathrm{DAR}}/{\mathit{\rho }}_{\mathrm{v},\mathrm{DAR}}\approx \mathrm{0.1}$ and noticeable biases compared to the model values, even when assuming perfect knowledge of temperature, imply that the sub-cloud DAR moisture observations should be applied judiciously to the problem of constraining surface RH. An obvious approach to improve the constraint from the sub-cloud measurements on RH would be to relax the assumption of the well-mixed sub-cloud layer; however, this is beyond the scope of this paper.

Table 3 Estimating near-surface water vapor density (ρ_v,DAR) and RH using the cloud-to-surface partial-column DAR humidity measurement. Reported DAR RH values assume perfect knowledge of surface temperature, and sensitivity to this assumption is quantified by the term $\partial \mathrm{RH}/\partial T$. Note that for the GATE case, the partial column often extends well into the middle and upper troposphere, and therefore the inferred surface RH is highly dependent on the assumed sub-cloud humidity profile shape. For all other cases, a well-mixed PBL is assumed for the shallow sub-cloud column.

4.3 Inferring in-cloud temperature

The second example application of DAR retrievals explored here involves utilizing the high-resolution profiles of water vapor inside of clouds and precipitation to estimate vertical temperature profiles by assuming that the volume is saturated with respect to liquid water. In the low latitudes it is observed that the spatial variability in temperature is significantly smaller than the variability in water vapor due to the large Rossby radius of deformation (Schneider, 1977). Therefore a constrain on the in-cloud temperature profile in the subtropical marine boundary layer will provide a strong constraint on the domain mean. For this application, we focus on the retrievals from the RICO case, which involves only liquid clouds and provides roughly 2 km of continuous humidity profiling within the PBL. Generally, we expect an RH of 100 % inside of cloud layers, but for precipitating volumes this can deviate significantly from saturation. Of course, from the perspective of the DAR, it is difficult to unambiguously distinguish between radar returns coming from range bins that are mostly cloud versus mostly rain. Furthermore, within a single radar footprint, it is likely that both cloudy and rainy regions are averaged together at a single height, in which case the rain signal will likely dominate due to the much larger average particle size. Therefore, by using the average DAR humidity profile acquired after a large along-track averaging distance, we are likely to be combining fully saturated and sub-saturated volumes and therefore potentially complicating the conversion from water vapor density to temperature.

In the fully saturated case, the temperature T inferred from the DAR measurement at each height is simply the one that satisfies the relationship ρ_v(z)=ρ_s(T), where ρ_v(z) is the retrieved humidity value and ${\mathit{\rho }}_{\mathrm{s}}\left(T\right)=e_{\mathrm{s}}\left(T\right)R_{\mathrm{v}}^{-\mathrm{1}}{T}^{-\mathrm{1}}$ is the saturated absolute humidity. Here e_s(T) is the saturation vapor pressure curve for liquid water determined by the Clausius–Clapeyron relation, and R_v=461.5 K⁻¹ kg⁻¹ is the water vapor gas constant. The results of this analysis using the RICO LES and simulated DAR retrievals are presented in Fig. 13, where we perform the conversion of humidity to temperature for three different profiles: the DAR retrieval, including random uncertainty, from Fig. 12; the conditionally sampled DAR truth profile from Fig. 5a; and the model humidity profile conditioned on grid cells with LWC > 0.01 g m⁻³. Note that the temperature profiles inferred from the DAR retrieval and truth humidity profiles incorporate non-uniform beam-filling effects, both in the horizontal and range dimensions, while the conditional model profile is at the model spatial resolution. Furthermore, the DAR truth and conditional model temperature profiles represent two different ideal measurements, with the former corresponding to a perfectly unbiased measurement made using the notional radar parameters in Table 2 and the latter a hypothetical measurement at the LES resolution requiring range and horizontal resolutions that are practically achievable.

Figure 13 Using DAR humidity observations to estimate temperature profiles in RICO. (a) LES relative humidity profiles conditioned on LWC > 0.01 g m⁻³, including the conditional mean (black line) and full model conditional range (gray region). (b) Comparison of model temperature profiles, including the conditional (LWC > 0.01 g m⁻³) mean and full range, with the values inferred from the following profiles assuming saturation: (dashed line) the model humidity profile conditioned on LWC < 0.01 g m⁻³, (red circles) the conditionally sampled DAR truth profile from Fig. 5a, and (blue circles) the DAR retrieval from Fig. 12.

For comparison with the humidity-inferred temperature profiles, we present in Fig. 13a the LES mean RH profile, as well as the range at each height, both conditioned on LWC > 0.01 g m⁻³​​​​​​​. Below the lifting condensation level around 500 m, we do not expect even the ideal humidity-inferred temperature measurements to agree with the model profile, as the liquid within this region is purely precipitation falling through a sub-saturated volume. Furthermore, while the mean conditional model RH profile remains near 100 % between the lifting condensation level (0.5 km) and cloud top (2.5 km), the conditional RH range is quite wide in this region, with the minimum conditional RH decreasing rapidly with height. We infer from this that there exist in the LES grid cells that are horizontally displaced from the saturated updraft region, especially near the cloud top, that consist of rain, but not cloud, liquid water. This is likely the result of vertical wind shear within the model forcing a horizontal separation between updrafts and downdrafts and has important consequences for the DAR measurement. Specifically, because the plumes are not highly resolved by the radar footprint, there will be many instances in which a single footprint encompasses updraft and downdraft regions and therefore both sub-saturated and fully saturated volumes. However, because the radar backscatter from precipitation particles is so much larger than for cloud droplets, the measured signal is heavily weighted towards the sub-saturated region within the footprint. This non-uniform beam-filling bias towards lower humidity implies a negative temperature bias when assuming that the local environment is fully saturated, leading to the divergence of the DAR retrieval and truth profiles from the model values above 1 km in Fig. 13b.

While it is theoretically possible to mitigate this negative temperature bias by decreasing the radar footprint size significantly, such a solution would require an unrealistically large antenna. Therefore, for practical applications, it would be necessary to incorporate knowledge of the degree of non-uniform beam-filling and precipitation amount in order to adjust the RH = 100 % assumption as a function of height within the cloud layer. Furthermore, this approach in general requires knowledge of the lifting condensation level – which could in principle be estimated from the DAR cloud-to-surface column IWV – in order to identify the sub-saturated region below the cloud base. Note, however, that we do not expect that similar non-uniformity biases will occur in non-precipitating clouds because all of the radar reflection will be coming from saturated volumes. Therefore, much care must be taken in inferring in-cloud temperature from DAR observations to account for the presence of precipitation. In fact, such studies would benefit greatly from airborne DAR observations of real precipitating liquid clouds with coincident, in situ observations of thermodynamic profiles from, for instance, dropsondes. In this case, the very high horizontal resolution offered by the airborne system would allow for systematic study of this non-uniform beam-filling effect as a function of horizontal averaging by combining along-track observations to form a larger effective footprint of variable size.

5 Conclusions

Spaceborne DAR holds considerable promise to fill gaps in the existing observing system for water vapor measurements, especially in the presence of clouds and precipitation. Through extensive LES and end-to-end radiative modeling, we have provided a high-fidelity assessment of the spatial coverage and systematic and random uncertainties that can be expected from a G-band DAR with realistically attainable performance metrics. Basic scaling arguments based on existing airborne G-band radar technology, as well as output from the radar simulator, suggest that a radar transmit power of at least 200 W and antenna diameter of 2 m are necessary for achieving adequate sampling of all cloud types explored in this work. Despite a fine range resolution of 50 m, which is necessary for mitigating non-uniform range-bin-filling biases from the water vapor retrievals, a radar with these parameters is able to probe well into deep convection, as well as sensitively detect non-drizzling marine stratocumulus cloud layers with a thickness of only 200 m. Future investigations focused on specific instrument architectures and radar waveform parameters are necessary to fully assess the impact of range-bin-to-range-bin correlations introduced by the range weighting function.

The notional DAR frequency channels are most sensitive to PBL humidity, with high-value water vapor retrieval products including in-cloud vertical profiles at 200 m resolution, sub-cloud IWV utilizing the strong surface reflection, and an upper-tropospheric IWV column between the radar and first cloud target. Additionally, the system provides a high precision measurement of total column water vapor when clouds are absent, which can be performed over land or ocean. As a general summary of the DAR capabilities for the five different LES cases, we find that in-cloud profiles with 10 %–20 % relative uncertainty can be retrieved for horizontal resolutions of 100–200 km, depending on the cloud type, and that upper-tropospheric and sub-cloud IWV measurements with 10 % uncertainty require along-track averaging on the scale of tens of kilometers. For all retrieval products, we find the systematic uncertainty to be generally much smaller than the random error, thanks to the improved spectral fitting routine that incorporates a third DAR frequency and new fit parameter (Battaglia and Kollias, 2019).

In addition to a detailed retrieval performance assessment, we have outlined several potential science applications of DAR observations, including the estimation of near-surface relative humidity using the sub-cloud IWV column, and inference of in-cloud temperature profiles using the DAR humidity profiles and assumption about the average saturation state. While much work remains to explore the limits of DAR for studies like these, the results demonstrate useful signals in the DAR observations for these applications. Beyond these examples and the basic utility of high-vertical-resolution water vapor information for modeling and assimilation applications, there are exciting opportunities to exploit DAR observations for probing the horizontal variability of water vapor in deep convective cold pools below the cloud anvil. Both continued simulation studies and instrument deployment on high-altitude airborne campaigns with in situ validation are essential for understanding the applicability of DAR observations to these science questions.

While this study chose the subtropical-to-tropical marine cloud transition to evaluate the potential capabilities of a spaceborne DAR, we can speculate as to the expected performance in different cloud regimes. We expect two primary limitations. First, we expect significant degradation of the radar sampling due to attenuation in the most heavily precipitating storms (e.g., mesoscale convective systems, hurricanes, and strong midlatitude frontal systems. Second, relative to the subtropical ocean many land regions will have insufficient scattered cumulus to acquire significant samples to perform vapor retrievals even with a scanning system. On the other hand, we expect that DAR retrievals will be quite useful in the midlatitude storm track regions over ocean in regions of post-frontal convection and stratocumulus which are frequently multi-layered and lightly precipitating, thereby providing ideal targets. In addition, cloudy polar and high-latitude scenarios will be particularly well suited to the DAR technique because of the relatively small amount of attenuation that the radar would experience.

Appendix A: Single-scattering properties

This appendix summarizes the approaches used to model all relevant single-scattering properties for the radar forward simulations.

Miles et al. (2000)Abel and Boutle (2012)Cox (1988)Wood (2011)Field et al. (2019)

Table A1 Summary of microphysical parameterizations for liquid- and ice-phase hydrometeors.​

^* The dendrite growth model uses the Leinonen and Szyrmer (2015) implementation of the algorithm presented in Reiter (2005).

Figure A1 Comparison of DDA, T-matrix, and Mie backscattering cross-section calculations for the two solid-phase hydrometeor species described in this work. The ice (a) and snow (b) species are treated as hexagonal columns and pristine dendrites for DDA calculations and as cylinders and oblate spheroids for the T-matrix approach, respectively. For the snow T-matrix calculations, a maximum diameter greater than 1.5 mm results in a spheroid aspect ratio that is beyond the convergence criteria for the T-matrix routine.

• Atmospheric gases. At millimeter-wave frequencies in the troposphere, the only relevant beam interaction with atmospheric gases is absorption by molecular oxygen and water vapor. We calculate the absorption properties using the millimeter-wave propagation model from the Earth Observing System Microwave Limb Sounder, which includes both line-by-line and continuum absorption contributions for each species (Read et al., 2004, 2006).

• Liquid hydrometeors. We treat the two LES liquid hydrometeor species, cloud and rain, as dielectric spheres and calculate the single-scattering properties using Mie scattering theory (Bohren and Huffman, 2004). The dielectric constant of pure water is calculated using the parameterization described in Liebe et al. (1991), which is an updated version of the classic model by Ray (1972). The particle size distributions (PSDs) for all hydrometeor species are parameterized using the modified gamma distribution,

  $$\begin{array}{}\text{(A1)}& N\left(D\right)={\displaystyle \frac{N_{\mathrm{0}}}{\mathrm{\Gamma }\left(\mathit{\nu }\right)}}{\left({\displaystyle \frac{D}{D_n}}\right)}^{\mathit{\nu }-\mathrm{1}}{\displaystyle \frac{\mathrm{1}}{D_n}}{e}^{-D/D_n},\end{array}$$

  where N₀ is the number concentration with units of particle number per unit volume, ν is the distribution shape parameter which can vary with species, and D_n is a characteristic diameter. Because the LES models used in this work include both single-moment and double-moment microphysics schemes, our radar forward simulator features two different parameterizations for calculating volume scattering quantities by integrating the single-particle, single-scattering quantities over a given PSD. For the two-moment schemes, the prognostic species number concentration and mass mixing ratio uniquely determine D_n once ν is prescribed and therefore permit straightforward calculation of the PSD-integrated quantities. For models employing one-moment schemes, however, an additional microphysical constraint is needed. In this work, all such constraints take on the form

  $$\begin{array}{}\text{(A2)}& N_{\mathrm{0}}=x_{\mathrm{1}}{D}_n^{x_{\mathrm{2}}}{e}^{-x_{\mathrm{3}}(T-T_{\mathrm{0}})},\end{array}$$

  where the parameters x₁, x₂, and x₃ have been determined from studies of comprehensive in situ microphysical observations, and T₀=273.15 K. All particle species except for cloud have one-moment parameterizations that assume an exponential PSD, or ν=1. For the cloud species, we assume a constant N₀ that is consistent with observations (Miles et al., 2000) and use the same shape parameter as for the two-moment scheme. The shape parameters and parameterization details for the cloud and rain species are detailed in Table A1.

• Ice hydrometeors. While the majority of convective regimes studied in this work include liquid-phase particles only, the case of deep convection necessitates modeling of radar scattering from three species of ice-phase particles: ice, snow, and graupel. Unlike cloud and rain droplets, the ice hydrometeors feature non-trivial relationships between particle mass m and maximum linear dimension D, so-called mass-dimensional relations:

  $$\begin{array}{}\text{(A3)}& m\left(D\right)=a_{\mathrm{m}}{D}^{b_{\mathrm{m}}}.\end{array}$$

  The mass-dimensional relation and PSD shape parameters for the two-moment implementation described here are taken from the Regional Atmospheric Modeling System (RAMS) two-moment microphysics parameterization (Walko et al., 1995; Meyers et al., 1997; Oue et al., 2020) and are presented in Table A1. While this implies that the microphysical parameterizations in the radar simulator differ from those in the LES models, this difference should have little impact on the fidelity of the forward simulated observations, since the LES models provide only bulk parameters describing the PSD.

  Before calculating single-particle scattering quantities, we must specify details of the three-dimensional shape and composition of the ice species. For graupel, we assume a homogeneous spherical geometry with reduced density relative to solid ice that is determined by the mass-dimensional relation, and we calculate effective dielectric and scattering properties using the Maxwell Garnet effective medium approximation and Mie scattering theory, respectively. For the ice and snow species, we assume particle shapes of hexagonal columns and dendrites, respectively, and generate discretized, three-dimensional crystal models using the approach of Leinonen and Szyrmer (2015). In this case, the mass-dimensional relations given in Table A1 provide the necessary constraint to fix the size-dependent aspect ratios of the columnar and dendritic crystals. We then calculate the scattering properties for ice and snow using the discrete dipole approximation (DDA) as implemented in the Amsterdam DDA (ADDA) solver of Yurkin and Hoekstra (2011). Since we expect the DDA results to converge to those of the T-matrix method (Mishchenko et al., 1996) for sufficiently small particle size, we specify an alternative bulk representation of the ice and snow species as homogeneous cylinders and spheroids with effective dielectric properties calculated again using the Maxwell Garnet approximation, and with the same size dependence of the aspect ratio as the 3D models. T-matrix calculations are then performed using the solver of Mishchenko and Travis (1998) as implemented in the Python package PyTMatrix by Leinonen (2014). Once convergence of DDA to T-matrix is established at some minimum particle size, the T-matrix values are used for all subsequent smaller particle sizes. For all ice hydrometeor species, the frequency-dependent dielectric constant values of solid ice are taken from Warren and Brandt (2008).

  Comparisons of the size-dependent backscattering cross section calculated using the Mie, T-matrix, and DDA methods for ice and snow are shown in Fig. A1. We note that DDA calculations for a particular particle size are performed for a single orientation that is chosen randomly for each particle realization and therefore do not correspond to orientation-averaged quantities. Instead, it is assumed that the integration over the PSD to derive volume scattering coefficients provides sufficient averaging over particle orientation.

• Ocean surface. Scattering from the ocean surface is treated within the geometrical optics (GO) approximation, in which backscatter is dominated by quasi-specular reflection from normally oriented wave facets (Kodis, 1966; Barrick, 1968; Valenzuela, 1978). At nadir incidence within the GO approximation, the ocean surface normalized radar cross section (NRCS) σ⁰ depends only on the mean square slopes of the ocean surface in the upwind and crosswind directions and the transmit frequency via the seawater dielectric constant. Here we use the classic parameterization of Cox and Munk (1954) to express the mean square slopes in terms of the near-surface wind and the millimeter-wave dielectric constant model of Meissner and Wentz (2004), for which we assume a constant salinity of 35 ppt. Lastly, we apply a fixed offset of +1.5 dB to the GO prediction for all frequencies and wind speeds in order to match recent airborne observations of σ⁰ at low wind speeds using the VIPR instrument (Roy et al., 2021).

Appendix B: Pulse-to-pulse correlations

In this appendix, we recount the derivation of the Doppler frequency, or equivalently velocity, spectrum and associated pulse-to-pulse correlations that result from the combined effects of the finite antenna beam width and fast satellite motion. Our approach is similar to previous treatments (e.g., Tanelli et al., 2002) but with modified notation to elucidate only the essential points. To begin, we examine the power spectral density observed by the spaceborne receiver due to radar backscatter from a stationary, volumetric target at a specific range r, assuming a transmit frequency of f₀:

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}S(r,f)& \propto \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{r-\mathrm{\Delta }r/\mathrm{2}}{\overset{r+\mathrm{\Delta }r/\mathrm{2}}{\int }}\mathrm{d}{r}^{\prime }\mathrm{d}{y}^{\prime }\mathrm{d}{x}^{\prime }{\mathit{\eta }}_{\mathrm{obs}}(f_{\mathrm{0}},{\mathit{r}}^{\prime })\\ & \left|F\right(f_{\mathrm{0}},x^{\prime },y^{\prime },r^{\prime }\left){|}^{\mathrm{4}}\mathit{\delta }\right(f-f_{\mathrm{0}}-\mathit{\delta }{f}_{\mathrm{D}}(y^{\prime },r^{\prime })).\end{array}\end{array}$$

Here we have again used the relation $z^{\prime }\approx R-r^{\prime }$ due to the narrow width of the beam and have assumed that the range weighting function can be approximated by a box of width Δr. Note that satellite-motion-induced Doppler shifts are introduced through the integration over the along-track distance variable y^′ using a delta function in frequency, where the shift itself is given by

$$\begin{array}{}\text{(B2)}& \begin{array}{rl}\mathit{\delta }{f}_{\mathrm{D}}(y^{\prime },r^{\prime })& ={\displaystyle \frac{\mathrm{2}{v}_{\mathrm{sat}}}{{\mathit{\lambda }}_{\mathrm{0}}}}\sin \mathit{\theta }(y^{\prime },r^{\prime })\\ & \approx {\displaystyle \frac{\mathrm{2}{v}_{\mathrm{sat}}{y}^{\prime }}{\mathit{\lambda }{r}^{\prime }}}.\end{array}\end{array}$$

Next, we assume a Gaussian beam shape with one-way intensity pattern

$$\begin{array}{}\text{(B3)}& \left|F\right(f_{\mathrm{0}},x^{\prime },y^{\prime },r^{\prime }){|}^{\mathrm{2}}=\exp \left(-{\displaystyle \frac{{x^{\prime }}^{\mathrm{2}}+{y^{\prime }}^{\mathrm{2}}}{{r^{\prime }}^{\mathrm{2}}{\mathit{\theta }}_{\mathrm{0}}^{\mathrm{2}}\left(f_{\mathrm{0}}\right)}}\right),\end{array}$$

where ${\mathit{\theta }}_{\mathrm{0}}=C_{\mathit{\theta }}{\mathit{\lambda }}_{\mathrm{0}}/D_{\mathrm{a}}$ and C_θ is a constant of order 1 that depends on the specific antenna implementation. Here, we assume a beam with C_θ=0.7. Then, the power spectral density in Eq. (B1) becomes, up to a constant multiplicative factor,

$$\begin{array}{}\text{(B4)}& S(r,f)=S(r,f_{\mathrm{0}})\exp \left(-{\displaystyle \frac{(f-f_{\mathrm{0}}{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_f^{\mathrm{2}}}}\right),\end{array}$$

where the spectrum width is given by

$$\begin{array}{}\text{(B5)}& {\mathit{\sigma }}_f={\displaystyle \frac{v_{\mathrm{sat}}{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\lambda }}_{\mathrm{0}}}}={\displaystyle \frac{v_{\mathrm{sat}}{C}_{\mathit{\theta }}}{D_{\mathrm{a}}}}.\end{array}$$

Equation (B5) reveals the important fact that the spectral width resulting from differential Doppler shifts in the along-track dimension for a fixed antenna size D_a is independent of the radar frequency. To get the corresponding width of the Doppler velocity spectrum, we use the relationships $\mathit{\delta }{f}_{\mathrm{D}}=\mathrm{2}v/{\mathit{\lambda }}_{\mathrm{0}}$ and ${\mathit{\sigma }}_v^{\mathrm{2}}=(\partial v/\partial \mathit{\delta }{f}_{\mathrm{D}}{)}^{\mathrm{2}}{\mathit{\sigma }}_f^{\mathrm{2}}$ to get ${\mathit{\sigma }}_v=v_{\mathrm{sat}}{\mathit{\theta }}_{\mathrm{0}}/\mathrm{2}$.

With the form of the power spectral density in hand, we calculate the signal autocorrelation function by invoking the Wiener–Khinchin theorem:

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}R\left(t\right)& =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}S(r,f)e^{\mathrm{2}\mathit{\pi }ift}\mathrm{d}f\\ & =R\left(\mathrm{0}\right)e^{-\mathrm{2}\mathit{\pi }({\mathit{\sigma }}_f^{\mathrm{2}}{t}^{\mathrm{2}}-i{f}_{\mathrm{0}}t)},\end{array}\end{array}$$

where t is the time lag between two sample sequences, and the power spectral density and autocorrelation functions in Eq. (B6) are intended to represent statistical ensemble averages. Finally, recognizing the signal voltage autocorrelation coefficient as ${\mathit{\rho }}_{\mathrm{s}}\left(t\right)=R\left(t\right)/R\left(\mathrm{0}\right)$, we calculate the power autocorrelation coefficient according to ${\mathit{\rho }}_{\mathrm{p}}\left(t\right)=\left|{\mathit{\rho }}_{\mathrm{s}}\right(t){|}^{\mathrm{2}}=\exp (-\mathrm{4}\mathit{\pi }{\mathit{\sigma }}_f^{\mathrm{2}}{t}^{\mathrm{2}})$ and obtain an estimate of the number of statistically independent radar samples N_i from (Papoulis, 1965; Doviak and Zrnić, 1993)

$$\begin{array}{}\text{(B7)}& N_{\mathrm{i}}=\sum _{m=-(N_{\mathrm{p}}-\mathrm{1})}^{N_{\mathrm{p}}-\mathrm{1}}{\displaystyle \frac{N_{\mathrm{p}}-\left|m\right|}{N_{\mathrm{p}}^{\mathrm{2}}}}{\mathit{\rho }}_{\mathrm{p}}\left(m{T}_{\mathrm{p}}\right),\end{array}$$

where N_p is the number of pulses transmitted and T_p is the pulse repetition interval. Note that the power autocorrelation coefficient provides a natural timescale for signal decorrelation, referred to as the time to independence, ${\mathit{\tau }}_{\mathrm{i}}=(\mathrm{2}\sqrt{\mathit{\pi }}{\mathit{\sigma }}_f{)}^{-\mathrm{1}}$. Then, Eq. (2) follows directly from Eq. (B7) and the definition of ρ_p(t).

Appendix C: Retrieval humidity profile interpolation example

In this appendix, we utilize a single-footprint measurement from the GATE case to clarify details of the interpolation-based humidity retrieval algorithm detailed in Sect. 2.4. These ideas are explored in Fig. C1. The essential task of the grid-selection and interpolation procedure is to take as an input profiles of observed reflectivity at the three DAR frequencies and calculate in an algorithmic fashion the resulting measurement axis $\mathit{r}=P_{\mathrm{r}}{\mathit{r}}^{\prime }$ and humidity retrieval axis $\mathit{z}=P_{\mathrm{z}}{\mathit{z}}^{\prime }$. To determine the measurement axis, we identify range-bin locations for which the observed reflectivity at all three frequencies is above the instrument detection limit dBZ_min(f) (solid black circles in Fig. C1b and c). Then, we must decide on the humidity profile vertical spacing, which we force to be an integer multiple of the radar range resolution. Recall that in Sect. 2.4 we defined this integer as the range-bin oversampling factor O=4. Given the set of possible humidity retrieval positions z^′, we identify actual positions for humidity retrieval (open black circles in Fig. C1b and c) as those for which there is at least one element of r in the list $[z^{\prime }-\mathrm{\Delta }r,z^{\prime },z^{\prime }+\mathrm{\Delta }r,z^{\prime }+\mathrm{2}\mathrm{\Delta }r]$. Note the offset of size Δr of the first humidity position in Fig. C1c from the surface position. This is due to the fact that the forward model must compute absorption for a given level using information from levels above it, and therefore the lowest level in the humidity profile is one range bin displaced above the actual surface.

Figure C1 Humidity interpolation procedure for the DAR retrieval, using an example from the GATE case. (a) Vertical profiles of observed reflectivity at the three DAR frequencies for the y=100 km position in Fig. 2c. In this case, the surface NRCS (not shown) is above the sensitivity threshold for all three frequencies. (b) Graphical depiction of the procedure for determining the projected measurement axis ($\mathit{r}=P_{\mathrm{r}}{\mathit{r}}^{\prime }$ – solid black circles) and humidity retrieval axis ($\mathit{z}=P_{\mathrm{z}}{\mathit{z}}^{\prime }$ – open black circles), as well as the humidity interpolation functions, which are the columns of the matrix A (differently colored, partial exponential profiles). (c) Close-up of panel (b) for clarity. (d) Resulting retrieved humidity profile (uncertainty not shown), including the N retrieved humidity elements of the state vector (x₂ – black circles) and the corresponding interpolated humidity profile Ax₂ (blue circles).