Vertical air motion retrievals in deep convective clouds using the ARM scanning radar network in Oklahoma during MC3E

Abstract

The US Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) program's Southern Great Plains (SGP) site includes a heterogeneous distributed scanning Doppler radar network suitable for collecting coordinated Doppler velocity measurements in deep convective clouds. The surrounding National Weather Service (NWS) Next Generation Weather Surveillance Radar 1988 Doppler (NEXRAD WSR-88D) further supplements this network. Radar velocity measurements are assimilated in a three-dimensional variational (3DVAR) algorithm that retrieves horizontal and vertical air motions over a large analysis domain (100 km $\times$ 100 km) at storm-scale resolutions (250 m). For the first time, direct evaluation of retrieved vertical air velocities with those from collocated 915 MHz radar wind profilers is performed. Mean absolute and root-mean-square differences between the two sources are of the order of 1 and 2 m s${}^{-\mathrm{1}}$, respectively, and time–height correlations are of the order of 0.5. An empirical sensitivity analysis is done to determine a range of 3DVAR constraint weights that adequately satisfy the velocity observations and anelastic mass continuity. It is shown that the vertical velocity spread over this range is of the order of 1 m s${}^{-\mathrm{1}}$. The 3DVAR retrievals are also compared to those obtained from an iterative upwards integration technique. The results suggest that the 3DVAR technique provides a robust, stable solution for cases in which integration techniques have difficulty satisfying velocity observations and mass continuity simultaneously.

Introduction

The representation of deep convection at cloud-resolving and global climate model scales (CRMs and GCMs) remains a serious challenge (Lin et al., 2006; Jakob, 2010). Part of the challenge can be attributed to the lack of comprehensive observations of dynamics and microphysics in these vigorous cloud systems (Ferrier, 1994; Milbrandt and Yau, 2005; Mrowiec et al., 2012; Donner et al., 2016). In particular, cloud dynamical insights may provide necessary guidance for improving these simulations at convection allowing scales and act as a basis for improving convective parameterizations at GCM scales (Lang et al., 2007; Wu et al., 2009; Nicol et al., 2015).

Despite the importance of vertical velocity measurements in deep convection, such measurements are difficult to acquire. Aircraft penetration of convective clouds offers the most direct method to measure these vertical air motions (Lenschow, 1976). However, practical hazards and operational costs have resulted in a valuable, but limited, dataset (e.g., Byers and Braham, 1948; LeMone and Zipser, 1980; Donner et al., 2001). Recent studies using profiling Doppler radars have suggested an ability to retrieve vertical velocities in convective clouds with an uncertainty of the order of 1–2 m s${}^{-\mathrm{1}}$, thus offering a viable substitute for in situ aircraft measurements (Jorgensen and LeMone, 1989; Cifelli and Rutledge, 1994; May and Rajopadhyaya, 1999; Williams, 2012; Heymsfield et al., 2010; Giangrande et al., 2013a; Kumar et al., 2015; Giangrande et al., 2016). Furthermore, profiling radars provide a high degree of detail of convective clouds in time and height, and can sample even the most intense convective cores. However, profiling radars potentially have a limited role for direct (rather than statistical) cloud model constraint due to their narrow view of these large, three-dimensional systems.

Scanning Doppler radars and the use of multi-Doppler retrieval techniques may help overcome known in situ aircraft and profiling radar sampling limitations. In addition to improved spatial representation of deep convection, this measurement approach offers an ability to document the three-dimensional structure of updrafts and downdrafts. Unfortunately, multi-Doppler retrieval applications are not straightforward. First, the number of radars in the network and their respective locations has a direct impact on the retrieval quality. Second, distributed Doppler radar networks, including mobile radar deployments, are not widely available or standardized. Operational radar networks often provide inadequate coverage throughout the depth of deep convective clouds, particularly at cloud top, necessary to constrain traditional vertical velocity retrievals at higher altitudes. Finally, since retrieval methodologies evolve, it is often difficult to establish a consensus pick among the many versions of multi-Doppler retrieval techniques that have been proposed.

In simpler terms, these retrieval techniques can be categorized as either “iterative” or “simultaneous”, based on their treatment of mass continuity. Iterative techniques solve the integral mass continuity equation throughout the column, contingent on known vertical velocity boundary conditions at the bottom level (i.e., upwards integration) and/or top level (i.e., downwards integration) (e.g., O'Brien, 1970; Ray et al., 1980; Protat and Zawadzki, 1999). This requires an estimate of horizontal wind divergence at each level made in a previous step, hence the non-simultaneity of iterative techniques (Dowell and Shapiro, 2003; Potvin et al., 2012a). By their nature, iterative upwards/downwards integration techniques propagate information in one direction, thus errors in horizontal wind divergence accumulate throughout the column, which in turn leads to larger errors in vertical velocity (e.g., Ray et al., 1980).

In contrast, simultaneous techniques treat mass continuity similar to other analysis constraints by inserting it directly into the cost function. This avoids accumulation of errors throughout the column since mass continuity is analyzed everywhere simultaneously. Moreover, these techniques are known to mitigate retrieval instabilities in poorly constrained regions like the dual-Doppler radar baseline (Bousquet and Chong, 1998; Dowell and Shapiro, 2003). Since simultaneous techniques are by definition 3DVAR techniques, we will refer to them as such throughout the remainder of this study. The 3DVAR approach has previously been shown to provide more accurate dual-Doppler retrievals than traditional techniques in observing system simulation experiments (OSSEs; Gao et al., 1999; Potvin et al., 2012a).

Several studies have investigated multi-Doppler wind retrieval uncertainties by identifying or utilizing (i) the importance of Doppler radar measurement errors and beam geometry (e.g., Doviak et al., 1976; Nelson and Brown, 1987; Matejka and Bartels, 1998; Bousquet et al., 2008), (ii) the influence of radar data objective analysis (e.g., Clark et al., 1980; Gal-Chen, 1982; Testud and Chong, 1983; Chong et al., 1983; Given and Ray, 1994; Majcen et al., 2008; Shapiro et al., 2010; Collis et al., 2010), and (iii) OSSEs (e.g., Fanyou and Jietai, 1994; Gao et al., 1999; Liou and Chang, 2009; Potvin et al., 2012b; Potvin and Wicker, 2012). However, few studies have compared practical retrieval performance to other independent air motion estimates from aircraft or ground-based profiling radars (e.g., Collis et al., 2013; Newsom et al., 2014). Collis et al. (2013) recently compared iterative dual-Doppler wind retrievals in tropical convection with those from a collocated dual-frequency wind profiler. However, because of the suboptimal location of the wind profiler near the dual-Doppler baseline, several assumptions had to be made before evaluating the two datasets.

While 3DVAR wind retrievals have been studied using OSSEs, an implementation, verification and sensitivity analysis using independent datasets from actual observations is noticeably missing. The US Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) program provides an excellent opportunity to investigate the benefits and relevant issues associated with multi-Doppler wind retrievals over the Southern Great Plains (SGP) in Oklahoma (Mather and Voyles, 2012). During the Midlatitude Continental Convective Clouds Experiment (MC3E), a joint field campaign between the DOE ARM program and the National Aeronautics and Space Administration (NASA) Global Precipitation Measurement (GPM) mission Ground Validation (GV) program (Jensen et al., 2016), the SGP site collected unique datasets from a distributed scanning Doppler radar network and radar wind profilers (RWPs) for several deep convective systems. These 3-D wind fields are strongly desired to analyze structures and characteristics of convective events (e.g., Liu et al., 2015; Donner et al., 2016). This study applies the 3DVAR radar wind retrieval to the MC3E deep convective events over the ARM SGP site and presents optimization of constraint weights used in the cost function. These 3DVAR retrievals are validated using data from the RWPs and compared with an iterative upwards integration technique.

Prominent convective events during MC3E, including a brief description and approximate time frame each event was sampled by UAZR-C1.

Event	Description	Time frame
		(UTC)
25 Apr 2011	Isolated, elevated convection	09:00–11:00
11 May 2011${}^{\mathrm{a}}$	Isolated convection, widespread stratiform precipitation	18:00–23:00
20 May 2011${}^{\mathrm{b}}$	Mesoscale convective system, squall line	06:00–16:00
23 May 2011${}^{\mathrm{b}}$	Isolated, severe convection	21:30–23:00
24 May 2011${}^{\mathrm{b}}$	Isolated, severe convection	21:00–22:30

${}^{\mathrm{a}}$ CSAPR-I7 nonoperational. ${}^{\mathrm{b}}$ XSAPR-I6 data recording issue.

This paper is organized as follows. A description of the dataset and radar data processing is presented in Sect. 2. Section 3 provides background information for 3DVAR wind retrievals from multiple scanning Doppler radars. Retrieval sensitivity and a method for producing physically sound wind fields is discussed in Sect. 4. Sections 5 and 6 presents 3DVAR wind retrieval results in the context of how they compare with those from independent collocated radar wind profilers as well as with an iterative upwards integration method. Section 7 is reserved for summary and concluding remarks.

Dataset and radar data processing

The MC3E took place during April–June 2011 in northern Oklahoma and surrounding states. A total of five events from MC3E were analyzed for this study and are listed in Table 1. These events represent a variety of warm-season convection over Oklahoma, including nocturnal elevated convection (25 April 2011), widespread stratiform precipitation with embedded convection (11 May 2011), mesoscale convective system (MCS) and associated squall line (20 May 2011), and isolated severe supercell thunderstorms (23–24 May 2011). The approximate time frame defining each event reflects profiling radar observations recorded at the SGP Central Facility (CF). The multi-Doppler radar wind retrieval in this study utilized plan position indicator (PPI) measurements operated by three X-band radars and one C-band radar from the ARM scanning precipitation radar network and one NEXRAD WSR-88D S-band radar. The ARM RWPs were used to evaluate these retrievals (Atmospheric Radiation Measurement (ARM) Climate Research Facility, 2011a, b, c). The environmental background wind fields are obtained from the ARM Merged Sounding value-added product that combines the observations from radiosonde soundings at the SGP CF (available every 3 h during MC3E), microwave radiometers, surface meteorological instruments, and European Centre for Medium Range Weather Forecasts (ECMWF) model output to produce a dataset at 1 min intervals and at 266 altitude levels (ARM, 1996).

ARM SGP site with locations of scanning and profiling radars surrounding the Central Facility (CF). Dashed black box in main panel corresponds to a 100 km $\times$ 100 km horizontal analysis domain (also shown in inset panel). Black circles are 40 km XSAPR maximum ranges (see Table 2). The separate analysis domain used in the sensitivity analysis (see Sect. 4) is shown as the dashed blue box surrounding the southeast radar wind profiler (UAZR-I9). The inset panel provides the large-scale view of the region including surface elevation in kilometers above mean sea level and the closest NEXRAD WSR-88Ds.

[Figure omitted. See PDF]

Operational parameters of the ARM XSAPR, CSAPR, UAZR, and NEXRAD WSR-88D (KVNX) radars during MC3E (convection mode).

Parameter	XSAPR	CSAPR	WSR-88D (KVNX)	UAZR (short-pulse/
				long-pulse modes)
Frequency (GHz)	9.4	6.3	2.85	0.915
Wavelength (cm)	3.2	5.4	10.5	33.0
PRF (kHz)	2.2	1.2	0.318–1.304 (short pulse)/	10.0/8.3
			0.318–0.452 (long pulse)
Pulse width (ns)	460	800	1570/4710	400/2833
Nyquist velocity (m s${}^{-\mathrm{1}})$	16.8	16.5	33.2	14.7/20.0
3 dB beamwidth (${}^{\circ }$)	1.2	1.0	0.9	9.0
Range resolution (m)	50	120	250	120/200
Temporal resolution${}^{\ast }$ (min)	6	7	5	0.1
Maximum range (km)	40	117	230	9.3/15
Number of elevations	22	17	14	–
Elevation range (${}^{\circ }$)	0.5–50.1	0.8–42.0	0.5–19.5	–

${}^{\ast }$ Temporal resolution of volume scan for scanning radars.

ARM scanning precipitation radar network

The ARM SGP site features a network of scanning Doppler and dual-polarization radars capable of providing coordinated coverage of cloud systems over a large domain. The locations of scanning and profiling radars around the SGP CF during MC3E are shown in Fig. 1. The radar facility includes a 6.3 GHz C-band scanning ARM precipitation radar (CSAPR-I7) and three 9.4 GHz X-band scanning ARM precipitation radars (XSAPRs, named I4, I5, and I6). During the MC3E, a specific deep convection volume coverage pattern (VCP) was implemented in an attempt to provide dense coverage throughout the depth of typical warm-season Oklahoma convection. The technical specifications of the ARM scanning radars are listed in Table 2.

Radar reflectivity observed by CSAPR-I7 was corrected for attenuation in rain using the CSAPR-I7 differential phase measurements as implemented using available open-source ARM python code utilities (e.g., Bringi and Chandrasekar, 2001; Giangrande et al., 2013b, 2014; Helmus and Collis, 2016). Because XSAPR reflectivity returns were significantly attenuated in rain, only XSAPR Doppler velocity measurements that are immune to partial attenuation were used in our retrievals. Aliased radial velocity measurements from all radars were corrected using the four-dimensional technique described in James and Houze (2001). Similar to Collis et al. (2013), this dealiasing technique was applied iteratively using horizontal wind profiles as obtained from the MC3E radiosonde network to produce robust results (e.g., Jensen et al., 2015). Each radar volume was manually inspected to check for conspicuous errors and artifacts.

NEXRAD WSR-88D radar

The NEXRAD WSR-88D S-band radar network surrounding the SGP site provides additional coverage and robust, unattenuated reflectivity measurement constraints for each event listed in Table 1 (e.g., Crum and Alberty, 1993). This was especially important for the 11 May 2011 event, as the CSAPR-I7 was nonoperational. The absence of the CSAPR implies these WSR-88D measurements carried additional weight in our retrievals, specifically since NEXRAD reflectivity factor measurements are the only ones not susceptible to attenuation in rain when paired with the remaining XSAPRs. As with the ARM radars, reflectivity factor and radial velocity measurements were corrected using similar methods.

Horizontal distributions of (a) nearest-neighbor distances and (b) nearest-neighbor weights within 20 km $\times$ 20 km surrounding the SGP CF for CSAPR-I7 (C), XSAPR-I4 (X), and KVNX (S) at 0 km (surface), 2 km, 6 km, and 10 km a.g.l.

[Figure omitted. See PDF]

The closest NEXRAD WSR-88D radar to the SGP CF, KVNX, located approximately 56 km west of the SGP CF, was used (Fig. 1). This relatively large distance, coupled with the 0.5${}^{\circ }$ base elevation scan of KVNX, ensures that its transmitted pulses are already 1 km above the surface directly over the CF. Figure 2a shows the nearest-neighbor distance between radar gates and grid points for the 20 km $\times$ 20 km region surrounding the CF assuming standard atmospheric refraction and Earth curvature (e.g., Doviak and Zrnić, 1993). The circular features seen in most cross sections are a result of discrete elevation scans. Between the surface and 2 km a.g.l., the ARM radars have enhanced coverage compared to KVNX. In particular, the ARM radars provide coverage that is ideal for characterizing the planetary boundary layer (PBL) since nearest neighbors are typically less than 150 m away from each other within this layer. At heights above approximately 2 km a.g.l., KVNX becomes increasingly valuable, especially for grid points close to and directly above the ARM radars. The dark red shades ($d$ $\ge$ 2 km) seen in the XSAPR panels of Fig. 2a highlight the radar cone of silence, a measurement gap due to no elevation scans past 50${}^{\circ }$ (see Table 2).

Mapping to the Cartesian coordinate grid

Radar reflectivity and Doppler velocity data from the ARM scanning and NEXRAD radars were mapped to a common Cartesian analysis domain. The domain for this study covers 100 km $\times$ 100 km $\times$ 10 km in meridional, zonal, and vertical extent, respectively, centered around the SGP CF, with 250 m grid spacing in each dimension. The horizontal area covered by the grid approximately encloses all available XSAPR coverage as shown in Fig. 1. Since the surface elevation of the analysis domain varies less than 30 m over its entire extent, this study neglects nuances associated with complex terrain (e.g., Chong and Cosma, 2000; Liou et al., 2011). The all-radar data are mapped using a single-pass isotropic Barnes distance-dependent weight (Barnes, 1964) with a constant smoothing parameter $\mathit{\kappa }=\mathrm{2}$ km${}^{\mathrm{2}}$ and convergence parameter $\mathit{\gamma }=\mathrm{0.5}$ (e.g., Trapp and Doswell, 2000; Majcen et al., 2008): $$\begin{array}{cc}{\displaystyle }{W}_{i,q}\left(d\right)=\exp \left({\displaystyle \frac{-d_{i,q}^{\mathrm{2}}}{\mathit{\kappa }{\mathit{\gamma }}^{p-\mathrm{1}}}}\right)\forall & {\displaystyle }i=\mathrm{1},\mathrm{\dots },n\\ & {\displaystyle }\text{and}& {\displaystyle }q=\mathrm{1},\mathrm{\dots },Q.\end{array}$$ Here $W_{i,q}$ is the weight for grid point $i$ and radar gate $q$ separated by distance $d$ for single pass ($p=\mathrm{1}$). The cutoff distance defining the $Q$ closest radar gates is the distance where the weight effectively vanishes, which is $d\approx \mathrm{4}$ km. The nearest 200 radar data points within the cutoff distance are used for interpolation at each grid point, but the number of data points would be less than 200 in sparse radar data regions (e.g., far distance from the radar). Several choices for weighting functions and their free parameters are found throughout the literature (e.g., Cressman, 1959; Barnes, 1964; Pauley and Wu, 1990; Askelson et al., 2000, 2005; Askelson and Straka, 2005; Trapp and Doswell, 2000); however, the weighting function used in this study is desirable for the preservation of phase and amplitude information of the input radar data, as well as its relative insensitivity to the spatial characteristics of the input data (Trapp and Doswell, 2000).

UAZR-C1 observations on 25 April 2011: (a) reflectivity, (b) Doppler velocity, (c) spectrum width, and (d) corresponding vertical air motion retrieval. Within the melting layer (2.5–3.5 km a.g.l.) no vertical air motion retrieval was attempted.

[Figure omitted. See PDF]

ARM radar wind profilers

Four 915 MHz UHF-band ARM zenith-pointing RWPs (UAZRs, named C1, I8, I9, and I10) were located within the scanning radar network (as shown in Fig. 1). Each UAZR was operated in a deep convective mode during the MC3E (Tridon et al., 2013). Technical details for the convective modes are also listed in Table 2. Vertical air motion retrievals from the wind profilers follow the method outlined by Giangrande et al. (2013a) that merges these modes to a single vertical velocity retrieval field ($\sim$ 6 s, 120 m) that is assumed accurate to within 1–2 m s${}^{-\mathrm{1}}$ in deep convective drafts. Our study uses datasets from UAZR-C1 and UAZR-I9 to evaluate the 3DVAR multi-Doppler wind retrieval as those operating and/or hit by convective cells for the events of this study. An example of UAZR-C1 convective cloud observations on 25 April 2011 and the corresponding vertical air motion retrieval are shown in Fig. 3.

Inherent differences between profiling and scanning radar observations restrict a direct comparison between the two datasets. The primary differences are related to temporal and spatial alignments of sampled radar volumes, beam broadening effects, and potential radar miscalibration, as seen in high-temporal and sampling resolutions in Fig. 3. The 6 s temporal resolution of the UAZR dataset is considerably higher than that of the scanning radars, which return to the same location every 6–7 min. Furthermore, the 120 m vertical resolution of the UAZR dataset is higher than the discrete elevation sampling of the scanning radars, especially for higher elevation scans where consecutive elevations are separated by more than 5${}^{\circ }$. Beam width and beam broadening effects must also be taken into account. Therefore, a two-dimensional time–height median filter is first applied to the UAZR dataset in order to remove high-frequency time-dependent phenomena and small-scale turbulent structures unresolvable by the scanning radars. We use a filter that is 61 time profiles in width (6 min) by 7 range gates (840 m) in height.

Using the time record and the fixed location of the UAZRs, the closest grid column in space and time is identified. Surrounding each identified grid column, we define a radius of influence $R_{\mathrm{s}}$ set to 750 m. This was implemented to account for the spatiotemporal sampling differences between the two datasets, as well as the advection of the cloud system. At each grid level, the median value within $R_{\mathrm{s}}$ is used as the best estimate, and the range of values within $R_{\mathrm{s}}$ are used to estimate the variability (e.g., spatial uncertainty) of the scanning radars. Similarly, at each UAZR range (height) gate, the median value within the time window of the corresponding scanning radar data (herein the valid time) is used as the best estimate for the UAZR, and the range of UAZR values within the valid time is used to characterize its variability. Several common error statistics are used to compare the two datasets: mean bias deviation (MBD), mean absolute deviation (MAD), root-mean-square deviation (RMSD), Spearman's rank correlation ($\mathit{\rho })$, and the Pearson product-moment correlation ($r)$. Difference statistics are computed subtracting UAZR from scanning radar, so a negative bias implies that the scanning radar dataset underestimated the corresponding UAZR dataset.

3DVAR wind retrieval methodology

Results presented in this study capitalize on the physical constraints of radial velocity observations ($J_{\mathrm{o}})$, anelastic mass continuity ($J_{\mathrm{c}})$, surface impermeability ($J_{\mathrm{p}})$, background wind field ($J_{\mathrm{b}})$, and spatial smoothness ($J_{\mathrm{s}})$. Assuming that $u=[u_{\mathrm{1}},\mathrm{\dots },u_n]$, $v=[v_{\mathrm{1}},\mathrm{\dots },v_n]$, and $w=[w_{\mathrm{1}},\mathrm{\dots },w_n]$ are the eastward, northward, and vertical wind components on the $n$-point analysis grid, respectively, we have the cost function $${\displaystyle }J(u,v,w)=J_{\mathrm{o}}+J_{\mathrm{c}}+J_{\mathrm{p}}+J_{\mathrm{b}}+J_{\mathrm{s}}.$$ The optimal wind field solution is at the (global) minimum of $J$ which implies that the gradient of $J$ with respect to $u$, $v$, and $w$ vanishes. For applications requiring large-scale (e.g., $\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}$ variables) nonlinear cost functional minimization, it is often necessary to use an iterative conjugate-gradient algorithm (Navon and Legler, 1987). In Gao et al. (1999), where a similar cost function and conjugate-gradient minimization algorithm were used, $u$ and $v$ were found to be well recovered within the first 50 minimization iterations; however, $w$ lacked both coherency and strength until 200$+$ iterations. We use these values as a reference point for the minimum number of iterations required to minimize Eq. (2).

Radial velocity observation constraint: $J_{\mathrm{o}}$

With radial velocity observations ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}$ defined on the same grid there is no need for an observation interpolation operator found in general 3DVAR schemes, and the observation constraint in Eq. (2) is instead given by $${\displaystyle }{J}_{\mathrm{o}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{l=\mathrm{1}}^m\left[{\left(v_{\mathrm{r}},-,{\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}\right)}^T,{\mathbf{\Lambda }}_{\mathrm{o}},\left(v_{\mathrm{r}},-,{\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}\right)\right].$$ The sum is over the $m$ radars used in the retrieval. ${\mathbf{\Lambda }}_{\mathrm{o}}$ is the set of constraint weights belonging to the radial velocity observations. It is an $n\times n$ matrix analogous to the inverse observation error covariance matrix in general 3DVAR schemes. We assume that observation errors are uncorrelated, meaning that ${\mathbf{\Lambda }}_{\mathrm{o}}$ is a diagonal matrix. The diagonal elements of radial velocity observation weights ${\mathbf{\Lambda }}_{\mathrm{o}}$ (${\mathit{\lambda }}_{\mathrm{o}})$ are estimated from the maximum value of Eq. (1) at each grid point (Fig. 2b) and observational data quality based on normalized coherent power for each radar. This naturally gives more weights to CSAPR and XSAPR observations within the PBL, and effectively ignores mapped observations propagated into sampling gaps such as the cone of silence. Elements of retrieved radial velocity $v_{\mathrm{r}}$ are

$$\begin{array}{cc}{\displaystyle }& {\displaystyle }{v}_{{\mathrm{r}}_i}=\left(u_i,\sin ,{\mathit{\varphi }}_i,+,v_i,\cos ,{\mathit{\varphi }}_i\right)\cos {\mathit{\theta }}_i+\left(w_i,-,w_{t_i}\right)\sin {\mathit{\theta }}_i\\ & {\displaystyle }& {\displaystyle }\forall i=\mathrm{1},\mathrm{\dots },n,\end{array}$$ where $\mathit{\phi }$ and $\mathit{\theta }$ are radar azimuth and elevation pointing directions, respectively, and $w_t$ is the bulk hydrometeor fall speed parameterized using radar reflectivity, temperature, and air density (Caya, 2001).

Radial velocity observations collected from two or more radars sampling the same convective cloud system are used in Eq. (3). The radial velocity observations are assumed to be closely matched in time. We required that (a) both KVNX and CSAPR-I7 be available (except 11 May 2011) and initiate a volume scan 2 min or less apart, and (b) any complementary XSAPR input initiate from a volume scan 2 min or less from either KVNX or CSAPR-I7. These criteria are designed to mitigate the errors associated with unaccounted advection and evolution of the cloud system (e.g., Gal-Chen, 1982; Shapiro et al., 2009).

Anelastic mass continuity constraint: $J_{\mathrm{c}}$

Anelastic mass continuity is known to be an adequate assumption in deep moist convection (e.g., Ogura and Phillips, 1962; Lipps, 1990). The general form of the mass continuity constraint is given by $$J_{\mathrm{c}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{L}^{\mathrm{2}}{D}^T{\mathbf{\Lambda }}_{\mathrm{c}}D,$$ where elements of $D$ are the anelastic mass continuity term, $$D_i={\displaystyle \frac{w_i}{{\mathit{\rho }}_i^{\mathrm{air}}}}{\displaystyle \frac{\partial {\mathit{\rho }}_i^{\mathrm{air}}}{\partial z}}+{\displaystyle \frac{\partial u_i}{\partial x}}+{\displaystyle \frac{\partial v_i}{\partial y}}+{\displaystyle \frac{\partial w_i}{\partial z}}\forall i=\mathrm{1},\mathrm{\dots },n,$$ and all vanish if anelastic mass continuity is perfectly satisfied. $L$ is a length scale inserted to unify the dimensions and magnitude of $J_{\mathrm{c}}$ with $J_{\mathrm{o}}$ (e.g., Legler and Navon, 1991; Bousquet and Chong, 1998; Shapiro et al., 2009). For this study we set $L=\mathrm{250}$ m, which is the grid spacing. In Eq. (6) ${\mathit{\rho }}^{\mathrm{air}}$ is air density derived from the MC3E radiosonde profiles. Although ${\mathbf{\Lambda }}_{\mathrm{c}}$ has $n\times n$ elements, this study sets them all to a constant value of ${\mathit{\lambda }}_{\mathrm{c}}$ for the diagonal elements (Table 3).

Note that for iterative upward/downward integration techniques, the vertical extent of the analysis domain controls the possible integration directions. If cloud tops are not adequately contained within the domain, a top boundary condition becomes impossible to define, making downwards integration impractical. For warm-season convective clouds in Oklahoma, a domain extending upwards of 15 km a.g.l. may be necessary in order to use downwards integration; however, these heights are poorly sampled by the scanning radar network and therefore poorly constrained by observations (see Fig. 2). Furthermore, Collis et al. (2010) showed that radar mapping artifacts aloft where radar coverage is poor leads to minimum vertical velocity errors of the order of 2 m s${}^{-\mathrm{1}}$ at these heights. This is the primary reason for capping our analysis domain at 10 km a.g.l. Although the upper level near the domain top can lack observation, the 3DVAR technique can produce better estimation compared with the iterative integration techniques (Potvin et al., 2012a).

Surface impermeability constraint: $J_{\mathrm{p}}$

This study imposes surface impermeability (Scialom and Lemaître, 1990) as a vertical velocity boundary condition at the ground level. Surface impermeability dictates that $w$ must vanish at the surface so we write $$J_{\mathrm{p}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}^T{\mathbf{\Lambda }}_{\mathrm{p}}w.$$ It is treated as a pseudo-strong constraint by heavily weighting its impact on surface grid points; non-surface grid points should not be influenced and their weights in ${\mathbf{\Lambda }}_{\mathrm{p}}$ are set to zero. This study uses constant values of ${\mathit{\lambda }}_{\mathrm{p}}$ as the diagonal elements of ${\mathbf{\Lambda }}_{\mathrm{p}}$ (Table 3).

Summary of 3DVAR constraint weights for stable solution derived from sensitivity analysis. Study values indicate the weight values used in this study. Nominal values reflect those used in previous OSSE studies.

Weight	Analysis	Study	Nominal
${\mathit{\lambda }}_{\mathrm{o}}$	–	(0, 1)	1
${\mathit{\lambda }}_{\mathrm{c}}$	(250, 1000)	500	1
${\mathit{\lambda }}_{\mathrm{b}}$	(0, 0.5)	0.01	0.01
${\mathit{\lambda }}_{\mathrm{p}}$	–	1000	–
${\mathit{\lambda }}_{\mathrm{s}u}$	(0, 100)	1	1
${\mathit{\lambda }}_{\mathrm{s}v}$	(0, 100)	1	1
${\mathit{\lambda }}_{\mathrm{s}w}$	(0, 100)	0.1	0.1

Background wind field constraint: $J_{\mathrm{b}}$

Including a background constraint helps improve the wind field solution in data-sparse regions based on additional observations. The background horizontal wind components $u_{\mathrm{b}}$ and $v_{\mathrm{b}}$ are typically those from a Merged Sounding profile nearest the analysis time. Since vertical velocity information is unavailable from these sensors, the background constraint is written as $$\begin{array}{cc}{\displaystyle }{J}_{\mathrm{b}}=& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\left(u,-,u_{\mathrm{b}}\right)}^T,{\mathbf{\Lambda }}_{\mathrm{b}},\left(u,-,u_{\mathrm{b}}\right),+,{\left(v,-,v_{\mathrm{b}}\right)}^T\right.\\ & {\displaystyle }& {\displaystyle }\left..{\mathbf{\Lambda }}_{\mathrm{b}},\left(v,-,v_{\mathrm{b}}\right)\right].\end{array}$$ Since $u_{\mathrm{b}}$ and $v_{\mathrm{b}}$ are assumed to be free of systemic errors, they are given the same (constant) weight (diagonal elements of ${\mathbf{\Lambda }}_{\mathrm{b}}\equiv {\mathit{\lambda }}_{\mathrm{b}}$, Table 3).

Spatial smoothness constraint: $J_{\mathrm{s}}$

The spatial smoothness constraint is essentially a low-pass filter designed to dampen high-frequency perturbations in the wind retrieval. Similar to Gao et al. (1999), we define this constraint as second-order spatial derivatives of $u$, $v$, and $w$: $$\begin{array}{cc}{\displaystyle }{J}_{\mathrm{s}}=& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{L}^{\mathrm{4}}[{\left({\mathrm{\nabla }}^{\mathrm{2}},u\right)}^T{\mathbf{\Lambda }}_{\mathrm{s}u}{\mathrm{\nabla }}^{\mathrm{2}}u+{\left({\mathrm{\nabla }}^{\mathrm{2}},v\right)}^T{\mathbf{\Lambda }}_{\mathrm{s}v}{\mathrm{\nabla }}^{\mathrm{2}}v\\ & {\displaystyle }& {\displaystyle }+{\left({\mathrm{\nabla }}^{\mathrm{2}},w\right)}^T{\mathbf{\Lambda }}_{\mathrm{s}w}{\mathrm{\nabla }}^{\mathrm{2}}w].\end{array}$$ In addition to reducing noise, Eq. (9) is able to extrapolate a wind field solution into data-sparse or poorly constrained regions. For instance, it may encourage usable solutions along the dual-Doppler baseline or add retrieval value to regions in close proximity to or directly above a radar (Bousquet and Chong, 1998). Although each of ${\mathbf{\Lambda }}_{\mathrm{s}u}$, ${\mathbf{\Lambda }}_{\mathrm{s}v}$, and ${\mathbf{\Lambda }}_{\mathrm{s}w}$ matrices has $n\times n$ elements, this study uses constant values of ${\mathit{\lambda }}_{\mathrm{s}u}$, ${\mathit{\lambda }}_{\mathrm{s}v}$, and ${\mathit{\lambda }}_{\mathrm{s}w}$, respectively, and ${\mathit{\lambda }}_{\mathrm{s}u}$ and ${\mathit{\lambda }}_{\mathrm{s}v}$ take the same value (Table 3).

Empirical wind retrieval sensitivity analysis

Typically, the constraint weight matrices ($\mathbf{\Lambda })$ found in Eqs. (5), (8), and (9) are treated as adjustable parameters ($\mathit{\lambda })$, controlling the degree to which each constraint influences the final solution. In essence, the values prescribed to each weight are often determined through trial and error (e.g., Gao et al., 1999). Fundamentally, there exists a range of values for each weight that produces a physically sound wind field. A thorough sensitivity analysis could be used to determine this parameter space, but this is often ignored because studies typically consider theoretical wind retrieval performance by comparing it to a known truth field (e.g., model output in an OSSE). The weights optimized to minimize the residual error between the retrieved and truth wind fields are then adopted (e.g., Gao et al., 1999; Potvin et al., 2012a). For applications involving real radar datasets where no truth field is available, one must consider (i) determining the parameter space which produces physically sound wind fields and (ii) characterizing the solution spread within the parameter space determined by (i).

This section addresses these two points through an extensive sensitivity analysis within the experimental domain indicated by the dashed blue box in Fig. 1. This domain has the same 250 m grid spacing and vertical extent as the larger domain, but covers a smaller horizontal area of 20 km $\times$ 20 km. Utilizing a smaller domain for the sensitivity analysis reduces processing time and allows for the isolation of specific cloud type regimes (e.g., convective versus stratiform). Since convective air motion retrievals are the primary interest of this study, the sensitivity analysis was done during a time when intense convection filled the experimental domain on 23 May 2011, using scanning radar observations valid between 22:36 and 22:43 UTC.

Point (i) is addressed by answering the following two questions. The first is, how well does the wind retrieval satisfy radial velocity observations? The second is, how well does the wind retrieval satisfy anelastic mass continuity? The second question is particularly important in the context of numerical modeling and convective parameterizations.

The wind retrieval is said to satisfy the radial velocity observations of one or more radars if the RMSD between the retrieval and observation is within the uncertainty estimate of the observations themselves. Since it is impractical to account for all sources of error inherent in mapped radial velocity observations, we establish a range of uncertainty and require the RMSD to be within this range. We employed radial velocity measurement error of approximately 0.5 m s${}^{-\mathrm{1}}$, which is a common value for regions of low signal-to-noise ratio ($\sim$ 20 dB; Fang et al., 2004) and larger Doppler spectrum width (Doviak and Zrnić, 1993; Bringi and Chandrasekar, 2001). The additional uncertainty introduced when mapping irregular radial velocity data to a regular grid is estimated to be of the order of 1 m s${}^{-\mathrm{1}}$. Therefore, we consider the wind field to satisfy radial velocity observations if it produces a RMSD with one or more radars within 0.5–1.5 m s${}^{-\mathrm{1}}$, computed over the entire analysis domain. To determine the degree to which the wind field satisfies anelastic mass continuity, following Shapiro et al. (2009), we define the normalized mass continuity residual (NMCR) as $$\begin{array}{cc}{\displaystyle }{\mathrm{NMCR}}_i=& {\displaystyle }{D}_i^{\mathrm{2}}\left[{\left({\displaystyle \frac{w_i}{{\mathit{\rho }}_i^{\mathrm{air}}}},{\displaystyle \frac{\partial {\mathit{\rho }}_i^{\mathrm{air}}}{\partial z}}\right)}^{\mathrm{2}},+,{\left({\displaystyle \frac{\partial u_i}{\partial x}}\right)}^{\mathrm{2}},+,{\left({\displaystyle \frac{\partial v_i}{\partial y}}\right)}^{\mathrm{2}}\right.\\ & {\displaystyle }& {\displaystyle }{\left..+,{\left({\displaystyle \frac{\partial w_i}{\partial z}}\right)}^{\mathrm{2}}\right]}^{-\mathrm{1}}\forall i=\mathrm{1},\mathrm{\dots },n,\end{array}$$ where $D_i$ is given by Eq. (6). As NMCR approaches zero, anelastic mass continuity becomes perfectly satisfied. However, this is not necessarily desirable since this condition is not exactly satisfied in nature, and even if it were, discretization errors would prevent precise satisfaction of Eq. (6). Therefore, we propose a range for NMCR, averaged over the entire analysis domain, between 1 and 10 %, whereby anelastic mass continuity is said to be adequately satisfied.

3DVAR constraint weight sensitivity analysis for two metrics: CSAPR-I7 radial velocity RMSD (left column) and NMCR (right column). Sensitivity analysis is performed by perturbing (a, b) ${\mathit{\lambda }}_{\mathrm{c}}$ versus ${\mathit{\lambda }}_{\mathrm{b}}$ and (c, f) ${\mathit{\lambda }}_{\mathrm{c}}$ versus ${\mathit{\lambda }}_{\mathrm{s}u}$, ${\mathit{\lambda }}_{\mathrm{s}v}$, and ${\mathit{\lambda }}_{\mathrm{s}w}$ constraint weights. The nominal values for weights not being tested in a given panel are set to those used in previous OSSE studies (see Table 3).

[Figure omitted. See PDF]

The response of CSAPR-I7 ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}$ RMSD and NMCR to perturbing multiple constraint weights is analyzed. The results are shown in Fig. 4. We first discuss the impact of the continuity and background weights (Fig. 4a–b). What is immediately evident in Fig. 4a is the strong dependence of ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}$ RMSD on ${\mathit{\lambda }}_{\mathrm{b}}$, with little to no dependence on ${\mathit{\lambda }}_{\mathrm{c}}$. Even with only a factor of 2 increase in ${\mathit{\lambda }}_{\mathrm{b}}$, the wind retrieval diverges substantially from the radial velocity observations and converges towards the background wind field. As ${\mathit{\lambda }}_{\mathrm{b}}\to \mathrm{0.5}$, CSAPR-I7 ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}$ RMSD approaches the specified upper limit of 1.5 m s${}^{-\mathrm{1}}$. This is important to note since 3DVAR retrievals have been found to be relatively insensitive to minor changes (e.g., not orders of magnitude) in other constraint weights (e.g., Gao et al., 1999; Potvin et al., 2012a). However, in Fig. 4b, ${\mathit{\lambda }}_{\mathrm{b}}$ has a decreased effect on the degree to which the wind retrieval satisfies mass continuity. As expected, this is primarily controlled by ${\mathit{\lambda }}_{\mathrm{c}}$, not only within the continuity-background parameter space but also in the continuity-smoothness parameter space shown in Fig. 4c–f. NMCR is particularly sensitive to ${\mathit{\lambda }}_{\mathrm{c}}$ when ${\mathit{\lambda }}_{\mathrm{c}}$ < 250. Outside of this range, NMCR is generally more stable with respect to ${\mathit{\lambda }}_{\mathrm{c}}$ and NMCR is typically less than 20 %. However, as seen in the right column of Fig. 4, in order to obtain NMCR $\le$ 5 %, ${\mathit{\lambda }}_{\mathrm{c}}$ must generally be 500 or larger.

Similar to Fig. 4a, b but for an iterative upwards integration technique.

[Figure omitted. See PDF]

Unlike the continuity-background sensitivity analysis, both RMSD and NMCR metrics appear highly unstable in certain regions of the continuity-smoothness parameter spaces investigated in Fig. 4c–f. For ${\mathit{\lambda }}_{\mathrm{s}u}$ and ${\mathit{\lambda }}_{\mathrm{s}v}$, which control the degree of smoothing of the horizontal wind components in Eq. (9), CSAPR-I7 ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}$ RMSD becomes unstable as these two weights approach values of 400 and larger. A similar phenomenon occurs for NMCR in Fig. 4d. These highly unstable regions of the parameter space are likely the result of nonlinear effects introduced by the squared second-order partial derivatives defined in $J_{\mathrm{s}}$ and should be avoided altogether. For values of ${\mathit{\lambda }}_{\mathrm{s}u}$ and ${\mathit{\lambda }}_{\mathrm{s}v}$ below approximately 100, CSAPR-I7 $\stackrel{\mathrm{̃}}{v}$ RMSD is within 1.5 m s${}^{-\mathrm{1}}$ and relatively stable. However, the parameter space in which this holds true gradually shrinks as ${\mathit{\lambda }}_{\mathrm{c}}$ increases towards 1000. Mass continuity is also adequately satisfied for ${\mathit{\lambda }}_{\mathrm{s}u}={\mathit{\lambda }}_{\mathrm{s}v}$ < 100 and ${\mathit{\lambda }}_{\mathrm{c}}$ > 250, with NMCR typically less than 10 %. Results for ${\mathit{\lambda }}_{\mathrm{s}w}$ are similar to those of ${\mathit{\lambda }}_{\mathrm{s}u}$ and ${\mathit{\lambda }}_{\mathrm{s}v}$ except for one aspect. Since ${\mathit{\lambda }}_{\mathrm{s}w}$ controls the degree of smoothing of the vertical wind component in $J_{\mathrm{s}}$, it has little influence on CSAPR-I7 ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}$ RMSD since the vertical wind component is generally not well sampled by scanning radars. This manifests itself in Fig. 4e, which shows CSAPR-I7 ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}$ RMSD to have much less dependence on ${\mathit{\lambda }}_{\mathrm{s}w}$ compared to ${\mathit{\lambda }}_{\mathrm{s}u}$ and ${\mathit{\lambda }}_{\mathrm{s}v}$.

Each panel in Fig. 4 contains over 2000 wind field realizations, each of which was concurrently saved. Therefore, we compute the 3DVAR vertical velocity solution spread from these thousands of realizations, allowing us to address point (ii) above. The ranges of the optimized $\mathit{\lambda }$ values derived from this sensitivity analysis and $\mathit{\lambda }$ values used for the retrieval case studies are recorded in Table 3. It is found that within the range of constraint weights defined in the analysis column of Table 3, the vertical velocity solution spread is relatively narrow at 1.5 m s${}^{-\mathrm{1}}$. This provides a form of uncertainty estimate for the 3DVAR wind retrievals presented in this study. It follows that we expect the 3DVAR vertical velocity retrievals to be relatively stable over a large range of constraint weights, with an uncertainty estimate of the order of 1–2 m s${}^{-\mathrm{1}}$.

C-SAPR radial velocity RMSD and NMCR from the 3DVAR and iterative upward integration techniques.

Event (UTC)	3DVAR		Iterative upward integration
	Radial velocity	NMCR	Radial velocity	NMCR
	RMSD		RMSD
25 Apr 2011 (09:16–09:23)	0.97	9.13	1.62	38.56
11 May 2011${}^{\ast }$ (18:12–18:18)	1.10	5.63	1.62	40.17
20 May 2011 (10:37–10:48)	2.03	8.36	3.48	38.74
23 May 2011 (22:36–22:44)	1.31	7.25	2.39	38.16
24 May 2011 (22:12–22:20)	1.95	10.21	4.38	42.74

${}^{\ast }$ KVNX radial velocity RMSD.

Comparison with iterative upwards integration technique

This section investigates the benefits of the 3DVAR approach for convective events as compared to an iterative upwards integration technique. First, a similar sensitivity analysis was performed for the iterative upwards integration technique, the results of which are shown in Fig. 5. Similar to the 3DVAR results in Fig. 4a, CSAPR-I7 ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}$ RMSD is highly dependent on ${\mathit{\lambda }}_{\mathrm{b}}$ and less so on ${\mathit{\lambda }}_{\mathrm{c}}$. However, for continuity-background parameter spaces where ${\mathit{\lambda }}_{\mathrm{c}}>\mathrm{5}$ or ${\mathit{\lambda }}_{\mathrm{o}}>\mathrm{0.5}$, there are sharp increases in ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}$ RMSD exceeding approximately 2.5 m s${}^{-\mathrm{1}}$, well outside the 1.5 m s${}^{-\mathrm{1}}$ upper limit. In fact, the parameter space in which ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}$ RMSD is below 1.5 m s${}^{-\mathrm{1}}$ is very small (approximately ${\mathit{\lambda }}_{\mathrm{c}}<\mathrm{1}$ and ${\mathit{\lambda }}_{\mathrm{o}}<\mathrm{0.1}$), and when looked at together with the normalized mass continuity residual (NMCR), no continuity-background parameter space exists in which both metrics are reasonably satisfied for an iterative upwards integration technique. It is worth noting that as ${\mathit{\lambda }}_{\mathrm{c}}$ is increased, NMCR appears to approach an asymptote around a value between 10 and 15 %. This indicates that even in the parameter space where radial velocity observations are effectively ignored (e.g., ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}$ RMSD greater than 3 m s${}^{-\mathrm{1}})$, iterative upwards integration techniques still have difficulty properly satisfying mass continuity. Results are also poor for the continuity-smoothness sensitivity analysis, in particular they were more unstable, and therefore they are not shown.

Comparison of 20 May 2011 squall line wind retrieval between 3DVAR and iterative upwards integration technique, showing (a) radar reflectivity, (b, c) 3DVAR and iterative vertical air motion, respectively, (d, e) 3DVAR and iterative horizontal wind divergence, respectively, (f) CSAPR-I7 radial velocity RMSD profile, (g) NMCR profile, and (h) vertical air motion bias (3DVAR minus iterative; 2 m s${}^{-\mathrm{1}}$ bin width). Select heights for panels (a–e) are 1 km, 2 km, and 8 km a.g.l. Origin in panels (a–e) corresponds to CF.

[Figure omitted. See PDF]

Next, we compare the actual wind fields retrieved by both techniques to determine if the difference found in the sensitivity analysis impacts the wind retrievals. The RMSD of radial velocity and the NMCR are estimated from the 3DVAR and an iterative upward integration method for the five cases at times when a strong convective region passed over the SGP CF for a 20 km $\times$ 20 km domain centered around the SGP CF, and the result is listed in Table 4. The 3DVAR technique provides lower NMCR and radial velocity RMSD values than the upward integration technique for the five cases. We note that the radial velocity RMSD values for the upward integration technique from the 25 April and 11 May cases are very low, close to the values found with 3DVAR techniques. These two events are nocturnal elevated convection (25 April) and widespread stratiform precipitation with embedded convection (11 May), respectively. Both cases included narrow or weaker convective regions, and exhibit propagation speeds that are slower than in the remaining MC3E cases featuring isolated severe convective cells and organized MCS events. The results suggest that the upward integration technique is comparable to the 3DVAR approach for the two cases where the mass continuity would be satisfied, whereas the 3DVAR technique demonstrates an advantage for the severe convective events.

Detailed comparisons of the retrieved wind fields from the two techniques are performed using the squall line event on 20 May 2011 at 10:40 UTC. This event featured the largest areas of strong convection for this MC3E dataset and substantial surface wind convergence ahead of the convective line that was well sampled by the scanning radar network. Strong wind convergence at or near the surface was indirectly observed around 10:40 UTC by UAZR-C1 as strong upwards motion lasting close to 5 min. The upwards integration technique used here is ideal for surface-driven events since the horizontal wind divergence profile should be well defined, particularly near the lower boundary. Figure 6 presents the wind retrievals from these two retrieval techniques within the 20 km $\times$ 20 km area surrounding the CF. Both techniques retrieve similar wind convergence patterns and upwards motion at 1 km a.g.l., with the iterative upwards integral technique retrieving a slightly stronger convergence line near the surface and therefore enhanced upwards air motion near the surface. Both methods also satisfy radial velocity observations below 2 km a.g.l., as shown in Fig. 6f. As expected, there is a large discrepancy between the two retrieval techniques when it comes to satisfying mass continuity. At each analysis height in Fig. 6g, the 3DVAR retrieval is adequately satisfying mass continuity, with NMCR < 10 % at each height and NMCR $=$ 5 % over the entire domain. For the iterative upwards integration retrieval, NMCR never gets below 30 % at any given height, and over the entire domain NMCR $=$ 52 %. Differences in the vertical velocity field between these techniques becomes more pronounced with increasing altitude due in part to the iterative upwards integration retrieval not adequately satisfying mass continuity throughout the column. At 8 km a.g.l., the vertical velocity fields no longer exhibit similar spatial patterns or intensities.

Radar reflectivity comparisons between CSAPR-I7 and UAZR-C1 at three characteristic heights for all events. Error statistics have units of dBZ.

Event	Height	Sample	MBD	MAD	RMSD	$\mathit{\rho }$	$r$
	(km a.g.l.)	size
25 Apr 2011	2	28	$-$1.00	2.22	3.79	0.95	0.80
	6	23	$-$1.62	2.73	3.51	0.93	0.64
	8	24	$-$1.09	2.50	3.81	0.92	0.67
	All	864	$-$0.94	2.97	3.96	0.94	0.55
11 May 2011${}^{\ast }$	2	70	$-$2.11	2.45	3.57	0.86	0.47
	6	85	$-$2.04	2.23	2.98	0.62	0.44
	8	69	$-$0.41	1.62	2.11	0.87	0.38
	All	2462	$-$1.16	1.98	2.79	0.92	0.40
20 May 2011	2	52	$-$2.11	2.33	3.24	0.97	0.77
	6	68	$-$1.90	2.28	3.79	0.79	0.76
	8	54	$-$1.29	1.79	3.60	0.85	0.70
	All	1983	$-$1.66	2.32	3.71	0.90	0.62
23 May 2011	2	6	$-$1.10	3.29	3.68	0.73	0.53
	6	19	$-$3.86	4.02	4.83	0.86	0.59
	8	11	$-$0.94	1.97	2.34	0.73	0.39
	All	398	$-$1.75	3.18	4.03	0.82	0.34
24 May 2011	2	6	$-$0.98	2.45	3.03	0.66	0.48
	6	16	$-$1.51	3.41	4.85	0.90	0.69
	8	29	$-$1.31	3.77	4.45	0.72	0.58
	All	557	$-$0.93	3.50	4.58	0.86	0.44

${}^{\ast }$ Comparison between KVNX and UAZR-C1.

The accumulation of differences with height is most evident in Fig. 6f. For the iterative upwards integration retrieval, CSAPR-I7 ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{r}}$ RMSD quickly grows larger than 2 m s${}^{-\mathrm{1}}$ at heights above 5 km a.g.l., whereas the 3DVAR method is able to satisfy CSAPR-I7 radial velocity observations at almost all analysis levels. As a result, the large spread in vertical velocity differences between the two retrieval techniques shown in Fig. 6h primarily comes from analysis levels above 5 km a.g.l. Over the entire analysis domain, vertical velocity MBD, MAD, and RMSD are all large at 1.3, 5.8, and 7.7 m s${}^{-\mathrm{1}}$, respectively. Below 5 km a.g.l., the MBD, MAD, and RMSD decrease substantially to 0.5 m s${}^{-\mathrm{1}}$ (38 %), 2.6 m s${}^{-\mathrm{1}}$ (45 %), and 3.9 m s${}^{-\mathrm{1}}$ (51 %), respectively. Nonetheless, these still represent large differences in the wind fields retrieved by each technique. That the radar coverage near ground level is not sufficient for these retrievals (particularly, as viewed from the KVNX radar, Fig. 2b) could be an explanation for the unsatisfactory mass continuity behavior in the iterative upwards integral technique.

Radar reflectivity observations on 20 May 2011 at three characteristic heights: (a) 2 km, (b) 6 km, and (c) 8 km a.g.l. CSAPR-I7 and KVNX error bars indicate the full range of radar reflectivities within $R_{\mathrm{s}}$ and the volume scan time (valid time). UAZR-C1 observations have been filtered using a 61 $\times$ 7 time–height median filter. KVNX observations are shown between 10:00 and 11:00 UTC, exclusively.

[Figure omitted. See PDF]

Evaluation with collocated profiling radars

Radar reflectivity comparisons

The comparison method between scanning radar and wind profiler measurements described in Sect. 2.4 is evaluated by comparing radar reflectivity measurements observed by UAZR-C1 with those from CSAPR-I7. The RWP receiver is known to saturate below 1 km range in heavier precipitation (e.g., $R$ > 10 mm h${}^{-\mathrm{1}}$, or $Z$ $\sim$ 40 dBZ at 1 km range); therefore, comparisons are performed above this level. Below the melting layer, the two reflectivity time series are moderately correlated, with Spearman's rank correlation $\mathit{\rho }$, and the Pearson product-moment correlation $r$ typically above 0.8 and 0.6, respectively, and MAD below 3 dBZ. Table 5 provides a summary of each comparative statistic at three characteristic heights for all events. The characteristic heights represent those below the melting layer (2 km a.g.l.) and above the melting layer (6–8 km a.g.l.). We note that the comparisons typically deteriorate within the melting layer (2.5–3.5 km a.g.l.) due to additional wavelength-dependent scattering phenomena.

Figure 7 shows the reflectivity time series of both CSAPR-I7 and UAZR-C1 at the three characteristic heights on 20 May 2011. The two datasets are visually highly correlated, with correlation coefficients typically above 0.8, MAD below 2.5 dBZ, and RMSD below 3.8 dBZ. This particular event included the formation and subsequent passage of a squall line directly over CSAPR-I7 around 10:40 UTC. The large difference in reflectivity between CSAPR-I7 and UAZR-C1 during the period 10:30–11:00 UTC is a result of rain and additional radome-induced attenuation effects on the measured CSAPR-I7 reflectivity (that typically cannot be resolved even using dual-polarization corrections). This highlights the advantage of incorporating observations from a longer-wavelength radar such as KVNX that is less susceptible to attenuation in heavy rain. To demonstrate this, KVNX reflectivity is superimposed in Fig. 7 exclusively between 10:00 and 11:00 UTC to illustrate the usefulness it offers in terms of reflectivity observations and associated fall speed corrections within the squall line. Overall, the relative reflectivity time-series comparisons shown here still indicate that these two independent dataset records are reasonably well matched for vertical velocity comparisons. The observed reflectivity discrepancies between the scanning radars and UAZR (see MAD and RMSD values in Table 5) are arguably comparable to the calibration limits achievable using natural media under Oklahoma conditions (Ryzhkov et al., 2005; Giangrande and Ryzhkov, 2005).

Two independent vertical air motion retrievals on 25 April 2011 from (a) 915 MHz radar wind profiler (UAZR-C1) and (b–e) 3DVAR. Radar reflectivity background and vertical velocity contours of 4 (light), 6 (medium), and 8 m s${}^{-\mathrm{1}}$ (thick) are shown in all panels. Wind vectors are shown in every 3DVAR panel. The horizontal cross sections shown in panels (b–c) are at 4 and 8 km a.g.l., respectively, with dashed black lines indicating the corresponding vertical cross sections in panels (d–e). The 3DVAR retrieval was derived from scanning radar observations recorded between 09:16 and 09:24 UTC. The origin in panels (b–e) corresponds to the location of UAZR-C1.

[Figure omitted. See PDF]

Vertical air motion comparisons between 3DVAR and UAZR-C1 at three characteristic heights for all events. Error statistics have units of m s${}^{-\mathrm{1}}$.

Event	Height	Sample	MBD	MAD	RMSD	$\mathit{\rho }$	$r$
	(km a.g.l.)	size
25 Apr 2011	2	20	$-$0.46	0.65	0.92	0.54	0.53
	6	20	1.18	1.58	1.86	0.43	0.51
	8	20	0.43	2.17	3.09	0.23	$-$0.03
	All	676	0.50	1.63	2.22	0.14	0.11
11 May 2011	2	27	0.02	0.79	1.05	0.34	0.35
	6	27	0.66	0.93	1.08	0.51	0.55
	8	27	0.23	0.61	0.74	0.68	0.63
	All	897	0.45	0.87	1.07	0.48	0.48
20 May 2011	2	55	0.38	0.99	1.24	0.22	0.69
	6	58	0.35	0.86	1.22	0.44	0.61
	8	51	0.14	0.94	1.49	0.38	0.15
	All	1871	0.29	0.92	1.27	0.32	0.50
23 May 2011	2	6	0.54	1.01	1.45	0.33	0.33
	6	19	0.42	0.99	2.31	0.32	0.41
	8	11	0.50	1.88	1.89	0.21	0.40
	All	398	0.52	1.74	2.11	0.29	0.40
24 May 2011	2	6	0.36	0.98	1.21	0.34	0.31
	6	16	0.39	1.32	1.75	0.51	0.29
	8	29	0.49	1.87	2.02	0.42	0.35
	All	557	0.50	1.41	2.01	0.49	0.47

Vertical air motion comparisons

The 25 April 2011 event was the first coordinated aircraft–ground precipitation mission during MC3E. Convective cells developed during the nighttime across northern parts of Oklahoma and along an elevated front, aided by mid- to upper-level ascent associated with the passage of an upper-level trough. The convective cells were relatively shallow in depth. Figure 8a shows the time–height cross section of vertical air motion retrievals and corresponding reflectivity from UAZR-C1. The closest available 3DVAR retrieval and its reflectivity field are shown in Fig. 8b–e. The time axis in Fig. 8a is reversed to better represent what the updraft retrieved by UAZR-C1 would look like in the north–south vertical cross section in Fig. 8e. The 3DVAR retrieval used scanning radar observations recorded between 09:16 and 09:24 UTC, spanning a total of 8 min. The most prominent features in both retrievals are a deep updraft region above 3 km altitude and strong updraft values greater than 8 m s${}^{-\mathrm{1}}$. The local system advection speed was estimated to be 18 m s${}^{-\mathrm{1}}$ by Giangrande et al. (2013a), with a north-northeast direction inferred by comparing successive CSAPR-I7 reflectivity displays. The base of the updraft retrieved by UAZR-C1 is first seen at approximately 09:23 UTC, which is near the end of the 3DVAR valid time window. The base of the updraft in the 3DVAR retrieval is approximately 2 km south of UAZR-C1, and with the prescribed system motion would pass over UAZR-C1 2 min later. This 2 min (2 km) offset is consistent with the UAZR-C1 retrieval if we assume that the 3DVAR retrieval is valid at 09:21 UTC, which is within the valid time. The two independent vertical air motion retrievals are qualitatively consistent with one another in terms of the relative location of the main updraft, its base height and depth, and its overall intensity.

A more direct time–height comparison between these two methods for the same event covering 09:00–10:45 UTC is provided in Fig. 9. In this case, the UAZR-C1 data have been filtered using the 61 $\times$ 7 time–height median filter. The elevated updraft seen in Fig. 8 is easily identifiable in the 6 km and 8 km a.g.l. panels between 09:15 and 09:30 UTC for both retrievals, with each showing the updraft strength to be stronger at 8 km rather than 6 km a.g.l. Visually the two retrievals are reasonably correlated at each height, with $r=\mathrm{0.51}$ and MAD less than 2 m s${}^{-\mathrm{1}}$ at 6 km a.g.l. (Table 6). The updraft retrieved by the 3DVAR method as seen in Fig. 8 that was offset by approximately 2 min (2 km) from the UAZR-C1 location is well captured by the 3DVAR error bars between 09:15 and 09:30 UTC in Fig. 9b–c. Table 6 reports the remaining errors and correlations between the two methods for this event. At most characteristic heights, vertical velocity bias is near 0.5 m s${}^{-\mathrm{1}}$ and absolute error is less than 2 m s${}^{-\mathrm{1}}$. These are considered to be small errors in comparison to the intensity of convective draft regions (e.g., 10–20 m s${}^{-\mathrm{1}})$ and negligible when considering the inherent differences between the two retrieval methods and spatiotemporal mismatch. The correlation coefficients are moderate at heights 6 km a.g.l. and below (not all shown), with values of $\mathit{\rho }$ and $r$ between 0.4 and 0.6. At higher altitudes such as 8 km a.g.l., correlations are weaker and errors are larger, likely the result of the stronger dynamics aloft for this elevated convective event (errors attributed to larger gradients of wind found aloft and unaccounted storm motion). That said, the 3DVAR vertical air motion time series at 8 km a.g.l. still appears to show some skill, with a bias less than 0.5 m s${}^{-\mathrm{1}}$ and error bars indicating a better correlation than is otherwise shown at this height in Table 6.

Vertical air motion time series on 25 April 2011 at three characteristic heights: (a) 2 km, (b) 6 km, and (c) 8 km a.g.l. 3DVAR error bars indicate the full range of vertical velocities within $R_{\mathrm{s}}$ and the 3DVAR valid time. UAZR-C1 retrievals have been filtered using a 61 $\times$ 7 time–height median filter.

[Figure omitted. See PDF]

The 20 May 2011 event was the longest-lived propagating MCS sampled by the scanning radar network during MC3E. UAZR-C1 observed leading stratiform precipitation and shallow convection throughout 06:00–10:00 UTC, followed by deep convection between 10:00 and 11:00 UTC which ultimately produced a large region of trailing stratiform precipitation that existed over UAZR-C1 for another 4–5 h. The most interesting feature of this event from a wind retrieval standpoint was the development of a well-organized squall line, passing over UAZR-C1 around 10:40 UTC. Similar to Fig. 9, the temporal comparisons at three characteristic heights between UAZR-C1 and 3DVAR wind retrievals for this event are shown in Fig. 10, covering the 6 h between 07:00 and 13:00 UTC. The 3DVAR vertical velocities near this time at 2 km a.g.l. reach upwards of 13 m s${}^{-\mathrm{1}}$, similar to the instantaneous, unfiltered values retrieved by UAZR-C1 (not shown). The large range of vertical velocities in $R_{\mathrm{s}}$ (e.g., $\sim$ 8 m s${}^{-\mathrm{1}}$ at 2 km a.g.l.) between 10:15 and 10:45 UTC are an indication of the strong dynamics associated with the squall line. Overall, there is good agreement between the two methods surrounding the squall line as well as throughout the rest of the 6 h period. Vertical velocity correlations as high as 0.7–0.8 were found at select heights between 2 km and 8 km a.g.l. (not shown), with the entire event producing a moderate correlation of $r=\mathrm{0.5}$ (see Table 6). Vertical velocity errors were also relatively small, with biases approaching 0.3 m s${}^{-\mathrm{1}}$, absolute errors generally smaller than 1 m s${}^{-\mathrm{1}}$, and root mean square errors below 1.5 m s${}^{-\mathrm{1}}$.

Similar to Fig. 9 but for 20 May 2011.

[Figure omitted. See PDF]

Two independent vertical air motion retrievals on 23 May 2011 from (a) 915 MHz radar wind profiler (UAZR-I9) and (b–e) 3DVAR. Radar reflectivity background and vertical velocity contours of $-$4 (light), $-$6 (medium), and $-$8 m s${}^{-\mathrm{1}}$ (thick) are shown in all panels. Wind vectors are shown in all 3DVAR panels. Horizontal cross sections shown in panels (b–c) are at 2 and 6 km a.g.l., respectively, with dashed black lines indicating the corresponding vertical cross sections in panels (d–e). The 3DVAR retrieval was derived from scanning radar observations recorded between 22:36 and 22:43 UTC. The origin in panels (b–e) corresponds to the location of UAZR-I9.

[Figure omitted. See PDF]

The 23 May 2011 event was part of an active sequence of severe convective outbreak days over the central plains. A surface low-pressure system located over the Texas panhandle and the associated surface boundaries were focal points for late afternoon convection. The environmental forcing coupled with strong daytime heating led to significant instability in addition to deep layer shear consistent with the eventual development of strong, discrete supercells. Convection captured within the analysis domain developed ahead of a surface dry line in western Oklahoma, coinciding with the passage of a shortwave trough. Supercells propagated eastward into the analysis domain by 21:00 UTC, with UAZR-I9 observing intense, deep convection between 22:00 and 23:00 UTC. Near 22:35 UTC, UAZR-I9 retrieved strong downdrafts reaching the surface with magnitudes larger than 8 m s${}^{-\mathrm{1}}$, with the core of the downdrafts increasing in height in the 8–10 min that followed. These results are shown in Fig. 11a. The closest available 3DVAR retrieval, valid between 22:36 and 22:43 UTC, is shown in Fig. 11b–e. Due to the east-northeast propagation of the local cloud system, the time axis in Fig. 11a was reversed to better reflect the east–west cross section through the 3DVAR retrieval in Fig. 11d. Cloud advection speed was estimated to be 17 m s${}^{-\mathrm{1}}$, indicating that the total downdraft feature retrieved by UAZR-I9 between roughly 22:35 and 22:45 UTC (10 min) covered approximately 10 km in length. This length is consistent with the east–west length of the total downdraft feature retrieved by the 3DVAR method in Fig. 11d, which covers the zonal length roughly between $x=-\mathrm{6}$ km and $x=\mathrm{4}$ km. Furthermore, the behavior of the retrieved 3DVAR downdraft is consistent with that of UAZR-I9, namely a surface-bound downdraft which appears to elevate as time (displacement) increases.

Discussion and summary

This study addresses the utilization of the ARM SGP scanning radar network in retrieving robust air motion estimates in five deep convective cloud events observed during MC3E, with the inclusion of surrounding NEXRAD WSR-88D assets. The ARM SGP site promises to provide a unique convective air motion dataset moving forward due to its continuous operation. Plans are currently being made for a second C-band scanning dual-polarization Doppler radar to be placed due south of the SGP CF, which would provide an additional constraint for these 3DVAR retrievals in the future.

First, this study optimized constraint weights of the cost function used in the 3DVAR multi-Doppler wind retrieval. The sensitivity analysis of the constraint weights suggests that the 3DVAR vertical velocity retrievals are relatively stable over a large range of constraint weights, with an uncertainty estimate of the order of 1–2 m s${}^{-\mathrm{1}}$. The retrieved updrafts from 3DVAR retrievals are compared with an iterative upward integration technique and independent ARM wind profiler observations. Our comparisons in terms of the constraint weight sensitivity testing indicated that there is a large constraint weight parameter space in which 3DVAR retrievals are able to simultaneously satisfy radial velocity observations and mass continuity and give stable solutions for vertical velocity retrieval. This is in contrast to the iterative upwards integration technique that has difficulty properly satisfying mass continuity. These comparisons for the five MC3E cases suggested that the 3DVAR technique can produce smaller errors in updraft retrievals than the iterative upward integration technique. Particular emphasis for this improvement was on the severe convection events that included large areas of strong convection.

Additional focus on the squall line case on 20 May 2011 revealed that the two techniques retrieved similar vertical velocity spatial patterns, including a large region of upwards motion associated with the surface convergence zone. However, the magnitudes of the vertical velocities between the two methods were considerably different, with MAD and RMSD of the order of 3 and 4 m s${}^{-\mathrm{1}}$, respectively. These large differences were likely caused by the iterative upwards integration technique inadequately satisfying mass continuity. In particular, the mean NMCR was 52 % for the iterative upwards integration retrieval compared to 5 % for the 3DVAR retrieval, and the radial velocity RMSD quickly grew larger than 2 m s${}^{-\mathrm{1}}$ above 5 km a.g.l. for the iterative upwards integration retrieval.

The 3DVAR retrieval was also evaluated in terms of how well it behaved when compared to collocated column measurements of radar reflectivity (for alignment) and retrieved vertical velocity from wind profilers. Time–height comparisons showed good visual agreement between reflectivity measurements, which was reinforced by correlations greater than 0.8 at most heights, MBD near $-$1.5 dB, and MAD typically less than 3 dB. The spatial and temporal characteristics of 3DVAR vertical velocity retrievals were also generally in good agreement with those from the wind profilers. Prominent updraft and downdraft features retrieved by the UAZRs were repeatedly observed in the 3DVAR dataset. Time–height comparisons showed reasonable agreement for most events analyzed, with moderate correlations of the order of 0.5, MBD less than 0.5 m s${}^{-\mathrm{1}}$, MAD within 1 m s${}^{-\mathrm{1}}$, and RMSD generally less than 1.5 m s${}^{-\mathrm{1}}$. In the context of deep convective drafts, where velocities can easily exceed 15 m s${}^{-\mathrm{1}}$, these differences are arguably negligible.

One omission for this sensitivity analysis was that the advection and time evolution of convective clouds during the typical PPI volume scan window was not considered (e.g., Gamache et al., 1995; Protat and Zawadzki 1999; Collis et al., 2013). A radar PPI scan from the NEXRAD and ARM radars during MC3E generally took 6–7 min to complete, allowing for substantial cloud movement and evolution in faster-moving deep convective events. In this regard, the co-gridding of the network radars cannot represent an actual snapshot of the 3-D convective structure, ultimately limiting the ability for this 3DVAR approach to satisfy the mass continuity equation (e.g., Clark et al., 1980; Gal-Chen, 1982). A further analysis using radar forward simulator and high-resolution ($\sim$ 0.5 km) and frequent (every 20 s) model output will be needed to better address these issues and other potential source of errors (e.g., PPI strategy, radar beam width, sensitivity) in future retrievals.

All ARM datasets used for this study may be downloaded at http://www.arm.gov (ARM, 1996, 2011a, b, c).

The authors declare that they have no conflict of interest.

Acknowledgements

This paper was authored by employees (Pavlos Kollias, Scott Giangrande) of Brookhaven Science Associates, LLC under contract no. DE-SC0012704 with the US Department of Energy (DOE). The contribution of Scott Collis through Argonne National Laboratory was supported by the US Department of Energy, Office of Science, Office of Biological and Environmental Research, under contract DE-AC02-06CH11357. Data were obtained from the Atmospheric Radiation Measurement (ARM) Program, sponsored by the US Department of Energy, Office of Science, Office of Biological and Environmental Research, Climate and Environmental Sciences Division. Edited by: Gianfranco Vulpiani Reviewed by: four anonymous referees