The Effects of Spatial Interpolation on a Novel, Dual-Doppler 3D Wind Retrieval Technique

In this study, we examine the effects of spatial interpolation on 3D wind retrievals using observing system simulation experiments (OSSEs, e.g., Potvin et al., 2012a, b; Dahl et al., 2019). Under the OSSE framework, observations are simulated using a known source, meaning the accuracy of the resulting experiments may be quantified exactly (in contrast to experiments with real data, where the “truth” is unknown). The source of our OSSE data is an Advanced Regional Prediction System (ARPS, Xue et al., 2000) simulation of the 8 May 2003 Oklahoma City tornado, first described in Xue et al. (2014). The simulation contains a mature, nearly steady-state tornadic supercell shown in Fig. 1. The reader is referred to Xue et al. (2014) and Schenkman et al. (2014) for a detailed description of the supercell. The ARPS model contains a terrain-following grid with 53 vertical levels, increasing from $\sim$20 m spacing at ground level to $\sim$750 m spacing at the top of the domain, and a uniform 50 m spacing in the horizontal dimensions. The simulation begins at 2200 UTC and extends $20\times 60\times 80$ km in the $z$, $y$ and $x$ dimensions, respectively, having been successively downscaled through a series of three nested domains with 9 km, 1 km and 100 m horizontal grid spacings. Mesonet, rawinsonde and Weather Surveillance Radar-1988 Doppler (WSR-88D) data are all assimilated into the analysis, and the reader is referred to Xue (2001) and Xue et al. (2003) for additional details on the model physics.

Modelled wind fields from the 20th minute of the 50 m resolution ARPS simulation (2220 UTC) serve as the “true” atmospheric state²²2This analysis time is identical to that in Dahl et al. (2019), meaning retrieved winds may be qualitatively compared to their findings despite the smaller analysis domain used in their study., upon which radar data is emulated at various positions around the domain. More specifically, our OSSE experiments contain an ensemble of fifteen different dual-Doppler observation locations (Fig. 1). The azimuthal spacing between radar pairs is set to 51.7^∘, so that the model grid lies at the center of the respective dual-Doppler lobe in each experiment (defined where the cosine of the cross-beam angle $<$ 0.9, shown as black dotted lines in Fig. 1), eliminating errors due to poor cross-beam angles. We implement the experimental ensemble to reflect the various sampling geometries possible in operational domains, thereby ensuring the generality and repeatability of the findings. We choose to extend our analysis grid from the lowest grid point at 500 m above sea level (ASL, note the terrain varies gradually between $\sim$270–390 m ASL across the domain), up to the tropopause height ($\sim$15 km altitude) and over the entire ARPS domain (60 km meridionally and 80 km zonally). The chosen domain is significantly larger than previous rapid-scan OSSE experiments (e.g., the $5\times 20\times 20$ km analysis grid in Dahl et al., 2019). Radars are positioned at an altitude of 350 m, either 60, 70 or 80 km from the center of the domain (illustrated by the varying line styles in Fig. 1), resulting in baseline distances for dual-Doppler pairs of $\sim$50, $\sim$60 and $\sim$70 km, respectively. Furthermore, the farthest range for the most distant radars in our experiment is roughly 130 km. These settings reflect our focus on retrieval methodologies applicable to operational radar networks (as opposed to small-scale, rapid-scan analyses with mobile radars), where dual-Doppler baselines are commonly large³³3For example, the largest dual-Doppler baseline distance tested in this study corresponds to the distance between the two operational radars (70 km, Terrey Hills and Wollongong) in Sydney, Australia’s most populous city., and analysts wish to retrieve winds over large domains.

Radar data simulation is implemented identically to previous studies (e.g., Dahl et al., 2019; B22), and is described here for completeness. The radar beam propagates assuming standard atmospheric refraction, and accounts for the curvature of the earth using the 4/3 effective earth radius model (Doviak and Zrnic, 1993). Model points within each radar voxel are weighted according to their proximity to the center of the beam, effectively functioning as individual scatterers. Radar observations are calculated as a weighted average of each scatterer as follows,

Radar reflectivity is calculated based on the modeled hydrometeor mixing ratios according to Tong and Xue (2005), and all Doppler measurements coincident with $<$5 dBZ reflectivity are masked to mimic data gaps in non-precipitating regions (e.g., Potvin et al., 2012c). Finally, we assume all radar information is collected instantaneously to eliminate temporal discretization effects due to flow evolution and translation, and do not include any consideration on how hydrometeor terminal velocities affect the measured Doppler velocities (i.e., we calculate radial wind velocities directly from the modeled $u$, $v$ and $w$ wind components assuming particle velocity is nil). While both error sources must be accounted for in retrievals from observations [e.g., using a vertical vorticity constraint and advection correction procedure (Shapiro et al., 2009; Potvin et al., 2012b; Shapiro et al., 2010), along with a reflectivity-based terminal velocity correction (Atlas et al., 1973)], we control for these effects in our experiments to isolate the errors due to spatial interpolation alone. Wind retrieval outputs from all methods are masked using the 5 dBZ threshold from the true model reflectivity to permit the direct comparison of error statistics. These statistics include root mean squared errors of various quantities, including the wind components themselves, the total wind magnitude, and derived products such as vertical vorticity $\zeta\equiv\partial v/\partial x-\partial u/\partial y$ and horizontal divergence $\delta\equiv\partial u/\partial x+\partial v/\partial y$. Consistent with previous studies (e.g., Potvin et al., 2012c), centered finite derivatives with Neumann boundary conditions were used for these derived fields. Fractions Skill Scores (FSS) are also calculated on column-maximum, or horizontal slices of vertical velocity fields, to quantify how well each method retrieves updraft/downdraft regions (Roberts and Lean, 2008). FSS scores are well-suited for assessing the ability to retrieve large vertical motions when slight spatial offsets are acceptable, owing to its original use as a neighborhood verification method for NWP rainfall accumulations.⁴⁴4We average all FSS calculations over five spatial scales in our study (namely 1, 2, 3, 4 and 5 km), and found that our results were insensitive to shortening, lengthening or adding additional length scales (not shown). Such cases may involve operational nowcasting or informing convection parameterizations, where identifying the size and strength of updrafts/downdrafts may be the priority. Finally, backward trajectory calculations for retrieved wind fields in Section 33.1 were calculated similarly to Potvin et al. (2012c) by using a second-order Runga-Kutta scheme with an adaptive time step, along with trilinear interpolation for the surrounding velocity values.

The variational 3D wind retrieval method proceeds with a cost function similar to those in previous literature, albeit with two important distinctions: an altered observational constraint ($J_{o}$) to facilitate the various interpolation methodologies, and a novel denoising constraint ($J_{d}$). The overall cost function is given below,

Most recent 3D wind retrieval studies (e.g., North et al., 2017; Dahl et al., 2019; Gebauer et al., 2022), along with open-source wind retrieval software (PyDDA, Jackson et al., 2019), require that dual-Doppler wind data must be pre-gridded prior to analysis. In this study, we aim to investigate how this pre-processing step may degrade the resulting wind retrievals relative to the direct assimilation approach described above, and provide guidance on how common gridding techniques perform relative to each other. We also include “control” dual-Doppler observations, which are upscaled directly from the original, high-resolution (50 m horizontal) model grid using Cressman (1959) weighted average interpolation with a radius of influence equal to the isotropic, analysis grid spacing (500 m). These control experiments will illustrate the baseline performance of our wind retrieval code for “perfect” observations, before they are degraded by the spatial interpolation errors in the “radar sampled” experiments.⁷⁷7Note that retrievals with “perfect” control observations do not result in wind fields with zero errors due to the effects of discretization errors, ill-posed boundary conditions and model assumptions (such as the anelastic mass continuity approximation) in the retrieval methodology. Furthermore, we test two common radar gridding techniques involving Cressman weighted averages in two and three dimensions, both of which are detailed below.

In this section, we analyze the performance of the first of fifteen ensemble members pictured in Fig. 1 (i.e., experiment 1 of 15). We aim to derive a qualitative understanding of the performance of each retrieval method through the analysis of this single experiment, before providing more quantitative ensemble statistics in Section 33.2. Firstly, Fig. 2 illustrates the various radial velocity inputs for the northwesternmost radar in experiment 1. The control field in Figs. 2a and 2c has not undergone any radar sampling, rather, it is a direct projection of the true Cartesian wind components in the radial direction. This setup renders it a ground truth description of the underlying radial velocities, limited spatially to regions with radar reflectivity $\geq$ 5 dBZ.

Perhaps the most noteworthy signature in the control field is the pronounced velocity couplet in region 1 of Fig. 2a, indicating the presence of a strong low-level mesocyclone. The shape and azimuthal shear values across the couplet are largely smoothed in the pre-gridded radial velocities in Figs. 2b and 2e. These effects are due to spatial aliasing in the radar gridding methods, which reliably filters observational noise, at the expense of high-frequency information required to accurately resolve such small-scale features (B22). We hypothesize that this type of smoothing will hamper the ability of pre-gridded methods to retrieve important dynamical features such as mesocyclone rotation in the resulting 3D wind retrievals. In contrast, the simulated radar data used in the RA method (Fig. 2f) accurately resolves the velocity couplet at $\sim$60 km range in experiment one, albeit in the presence of notable speckle noise added during the radar simulation process. This level of observational noise is common throughout the simulated radar data (e.g., in region 2), and it is not clear a priori whether the benefits of finer data resolution will outweigh the added observational noise in the RA method.

Another important consideration for each wind retrieval method is the spatial continuity; or equivalently, the size and shape of data voids, in the radial velocity data. The “true” extent of the simulated radar echoes is shown in the control experiment (with the minimum detectable signal defined at 5 dBZ). The 3D Cressman method artificially extends valid data past its true extent (e.g., between regions 3 and 4 in Fig.  2), effectively “filling in” many true data voids. It is not obvious whether this data extension is advantageous in 3D wind retrievals, as the added observations may actually assist in retrieving valid winds in data voids, contingent on the quality of the extrapolated radial wind values. Velocity observations in the RA method contain large gaps aloft due to the spacing between constant elevation radar sweeps (e.g., region 3 in Fig. 2h), and this may have serious implications for the accuracy of retrieved winds if the interpolation operator ($\mathcal{C}$) and the smoothing constraints ($J_{s}$) cannot effectively propagate radar information into these data voids. Similarly, the radial velocities in the 2D Cressman method are severely limited above the highest tilt, and below the lowest tilt (e.g., region 2 in Fig. 2e and region 4 in Fig. 2g). A notable consequence of these spatial continuity considerations may be observed in region 4, where a signature indicating easterly inflow in the low levels is evident in the control observations (Fig. 2c). This signature is largely absent in the pre-gridded methods, due to over-smoothing in the 3D Cressman observations, and data gaps below the lowest tilt in the 2D Cressman observations. Once again, we hypothesize that the pre-gridded retrievals will be adversely affected by the omission of this important dynamical information in the underlying radial velocity data.

Horizontal cross sections of retrieved vertical velocity for each method are given in Fig. 3. Visual comparisons of both shading colors and streamlines between Figs. 3a and 3c indicate that the control experiment is able to reproduce the true 3D wind field very accurately. One notable exception in the accuracy of the control experiment are small patches of spurious vertical velocity retrievals, such as those in region 2 of Fig. 3a. These are due to boundary effects in data-sparse regions, and are an example of how errors are still present in 3D wind retrievals with “perfect” radial velocities. For these reasons, the accuracy of any 3D wind retrieval method simulated using practical radar scanning geometries should also be judged against such a control experiment, and not merely the “true” model winds. The reader is referred to Appendix B for further discussion on boundary effects in 3D wind retrievals.

All three retrievals that ingest radar simulated radial velocities (either pre-gridded in 3D/2D Cressman methods, or directly in the RA method) are able to reproduce the broad-scale southwesterly horizontal flow, along with the approximate size and shape of the updraft within the inset region in Fig. 3 (hereafter referred to as the analysis region). The largest errors arise in the vertical velocity retrievals for the pre-gridded methods, which largely underestimate the true magnitude of the updraft velocities (by more than a factor of 2 in the 3D Cressman method), and are unable to retrieve the small downdraft features around the edges of the main updraft (blue shading, northwest and southeast of region 1 in Fig. 3c). We attribute these poor vertical velocity retrievals to the errors introduced by pre-gridding the underlying radial velocity fields. More specifically, filtering and over-smoothing of the highest radial velocity magnitudes, along with the poorly resolved low-level convergence discussed in Fig. 2, both contribute to the large vertical velocity errors noted in Fig. 3.

Unlike pre-gridded methods, the RA method does accurately resolve the amplitude and extent of the main updraft region (with maximum updraft velocities in Fig. 3e approaching the true values evident in Fig. 3c). Furthermore, the position of the aforementioned downdraft features are retrieved well in this method, while slightly underestimating maximum downdraft speeds. We attribute this success to the higher resolution of the input radial velocity data, which allows the technique to resolve the small-scale and highly dynamic circulation within the analysis region. Another example of the retrieval of important small-scale dynamic features is the cyclonic rotation within region 1 in Fig. 3c. Streamlines in Fig. 3e indicate this feature is retrieved well in the RA method, where the easterly component of the mesocyclone rotation is clearly visible against the prevailing westerly background flow. Both pre-gridded methods were unable to retrieve any easterly wind component within the mesocyclone region, which has important implications for identifying supercell thunderstorms, which cause a disproportionate amount of large hail and significant tornadoes (e.g., Duda and Gallus, 2010; Smith et al., 2012).

Vertical cross sections of vertical velocities from experiment 1 in Fig. 4 further illustrate the deficiencies of the pre-gridded retrieval techniques. We will highlight these deficiencies in the context of two important dynamic features present in the analysis inset in Fig. 4c: a rear flank downdraft (RFD) signature close to the surface within region 1, and an overturning circulation resulting from an updraft/downdraft couplet within region 2. Both pre-gridded methods were unable to retrieve any downdraft areas within the inset region, meaning both the RFD and overturning circulation are absent from the resulting wind fields. However, the RA was again able to retrieve these features accurately in terms of their shape and position, whilst slightly underestimating the true vertical velocity maximum values. The retrieval of these features has important implications for nowcasting meteorological hazards, such as identifying potentially strong straight line winds at the surface due to the RFD, and the existence of favorable hail growth trajectories due to the overturning circulation (e.g., Kumjian et al., 2021).

While the wind retrieval clearly benefits from the added resolution of radial velocity observations in the RA method, this method still displays considerable limitations when compared to the control retrieval. Firstly, consider the 3rd highlighted region in Fig. 4c, which contains two spatially proximal, but distinct regions of very high vertical velocities ($>$ 30 m s^-1). The control method in Fig. 4c accurately retrieves these separated updraft regions, meaning this scale of motion (these features, and most convective towers, are approximately 1–2 km in diameter) is practically resolvable on the 500 m isotropic analysis grid used in this study. The RA method is once again able to retrieve the magnitude of this updraft feature in region 3 much more accurately than the pre-gridded methods, however, the two updraft features are aliased into one large updraft region in Fig. 4e. This spatial aliasing is present throughout the upper levels ($>$10 km), and is a result of the sparsity of observed radar information at these altitudes (refer to Fig. 2h). Clearly, the accuracy of the RA method is still fundamentally limited by the scanning geometry of the underlying radar observations, albeit to a lesser extent than the pre-gridded methods. Finally, we also note the poor performance of all three radar-derived retrieval methods in the forward flank of the storm (e.g., region 4 in Fig. 4). Here, spurious vertical velocity signatures are present throughout the depth of the storm, and are coincident with significant data gaps due to clear air in the low levels. We posit these errors are caused by a well-known problem in the wind retrieval literature, whereby the impermeability boundary condition ($w=0$ at $Z=0$) in the mass continuity constraint is poorly imposed in regions with data voids near the surface. Studies have shown these issues may be alleviated through the inclusion of a vertical vorticity constraint (Shapiro et al., 2009; Potvin et al., 2012b), which is not implemented here due to our sole focus on the effects of spatial interpolation.

The final qualitative assessment of the wind fields from experiment 1 is a comparison of simulated air parcel trajectories in Fig. 5 (the control experiment is omitted due to its similarity with the model truth). Trajectories are initiated within the inflow of the storm at 600 m, 700 m and 800 m altitudes, at 1 km spacing within the horizontal bounds of the analysis box shown in the inset in Fig. 3. Broadly speaking, trajectories calculated using the true wind fields in Fig. 5a follow two main paths: strong vertical ascent to the top of the domain ($>$8 km) within the updraft [mostly those initiated in the western (blue) and northern (orange) quadrants], or initial ascent followed by ejection into the forward flank of the storm in the mid-levels [4–8 km, mostly the eastern (green) and southern (red) quadrants]. The RA winds in Fig. 5d produce trajectories that are qualitatively similar to the true trajectories for all four quadrants. Most notably, the RA method reproduces the strong split between trajectories from the southern (blue) and western (red) trajectories, which is also observed in Fig. 5a. Many trajectories in the pre-gridded fields also exhibit very similar trajectories to those in Fig. 5a, however, there are also some notable differences. For example, the aforementioned split between the western and southern inflow regions is not observed in Figs. 5b and 5c. Rather, there exists a smooth transition between these pathways, resulting in a roughly continuous spatial distribution between the trajectories at the top of the storm and those ejected into the mid-levels.

Another notable difference in the pre-gridded retrievals (vs. the radar assimilation retrievals) is the comparative smoothness, or uniformity, of the trajectories, especially in the forward flank. This is once again a result of filtering of important information in the underlying gridded radial velocities, resulting in over-smoothed retrieved wind fields. Lastly, perhaps the most notable difference in the pre-gridded fields is the spurious collection of low-level trajectories which initiate from the northern quadrant of the initiation region. The presence of these trajectories indicates a lack of convergence within the main updraft, which ultimately results in the aforementioned underestimation of updraft strength for these methods. The lack of convergence in the low-levels was hypothesized through the analysis of Fig. 2, where the inflow signature was poorly resolved in the pre-gridded radial velocities. Fig. 5 qualitatively confirms this assumption, and emphasizes the importance of ingesting accurate radial velocity information, especially in the low levels.

Table 1 presents a summary of the experimental results for each retrieval method, including ensemble mean and standard deviation values. The control experiment achieves excellent results across all radar geometries tested in this study. Considering the average wind magnitude for the experiment is $\sim$51 m s^-1, the ensemble average RMSE of 2.00 m s^-1 results in an average percentage error of just $\sim$4% in wind magnitudes across the domain. A standard deviation of 0.05 m s^-1 for the 15 experiments also indicates the technique is not particularly sensitive to different radar positions or cross beam angles. As expected, the vast majority of error comes from the poorly observed vertical wind component, which has an average RMSE of 1.96 m s^-1, compared to the horizontal component RMSE’s of 0.23 and 0.22 for $u$ and $v$, respectively. We also note the control experiment is able to retrieve winds with an average RMSE of 2.68 m s^-1 in the analysis region, which indicates that the “perfect” radial velocities on a 500 m analysis grid are capable of resolving the scales of motion very well in this highly dynamic region of the storm. Finally, we observe that Fractions Skill Scores (FSS) for updrafts outperform those for downdrafts in the control experiment (0.95 vs. 0.80, respectively), and note this finding is consistent across the three radar observed methods. This is likely a result of the proximity of downdraft regions to the top of the analysis domain ($z$ = 15 km) at this time step, where boundary errors are prevalent (compare downdrafts in Figs. 4a and 4c).

As expected from the qualitative analysis in Section 33.1, the quality of retrievals is considerably degraded when winds are retrieved with simulated radar radial velocities, resulting in a 300% – 350% increase in error magnitudes. Interestingly, we did not observe reliable accuracy reductions for ensemble members with more distant radars (not shown), indicating the storm was adequately sampled by the volume coverage pattern at each of the three ranges used in this study (60, 70 and 80 km). Firstly, both pre-gridded methods score similarly in terms of overall RMSE scores (7.17 m s^-1 for 3D Cressman and 7.43 m s^-1 for 2D Cressman), however, ensemble standard deviations indicate that error scores are much more consistent in the former. This is easily explained by the superior error filtering qualities in the 3D Cressman analysis, which limits the amount of observational error propagating into the retrieved wind fields. While the smoothing properties of the 3D Cressman analysis are beneficial in the consistency of error scores, they also hinder the retrieval of strong winds in the analysis region, where the RMSE score is considerably higher at 11.61 m s^-1. The inability of the 3D Cressman method to retrieve winds in highly dynamic regions is further underscored by the poor updraft and downdraft FSS values (defined by the 5th and 95th percentile of column-maximum vertical velocities, $w\approx 8.5$ and $w\approx-3.4$, respectively) relative to the control experiment (0.70 vs. 0.95 for updrafts, and 0.26 vs. 0.80 for downdrafts). These results ultimately illustrate the trade-off between the 3D Cressman and 2D Cressman gridding methods: the former results in more consistent retrievals with slightly lower total error measures, but degrades the retrieval in the most dynamic regions of the storm, which are perhaps the most important in terms of understanding storm dynamics and nowcasting potential hazards.

The average RMSE score of 6.22 m s^-1 in the RA method supports the aforementioned qualitative improvements over the pre-gridded methods. The use of this technique results in an average reduction in errors of 18% and 22% compared with the 3D Cressman and 2D Cressman methods, respectively, relative to the control experiment. These domain-wide improvements relative to the pre-gridded methods are also reflected in RMSEs in the analysis region (8.01 m s^-1), and for individual wind components over the entire domain ($\sim$3.5 m s^-1 errors for $u$, $v$ and $w$). Interestingly, while the average vertical velocity error is slightly lower in the RA method, the majority of RMSE improvements come from reducing errors in the horizontal components, which are much lower than in the pre-gridded retrievals. This belies the qualitative observations made in Figs. 3 and 4, which showed considerable improvements in RA vertical winds, especially within the analysis region. An explanation for this seeming discrepancy is that high-amplitude fields are penalized more heavily by pixel-to-pixel accuracy measures such as RMSE, as slight spatial offsets in wind fields result in higher error scores compared to smoother fields such as the 3D Cressman method. This effect also accounts for the lower RMSE values observed in the 3D Cressman method compared to the 2D Cressman method, despite the considerable over-smoothing noted in the former in Section 33.1. We further this hypothesis by noting that the FSSs for updraft and downdraft regions are considerably higher for the RA method. The updraft and downdraft FSSs in the RA method are only 7% and 21% lower, respectively, than the control experiment, while the 3D Cressman and 2D Cressman methods show larger a reduction of 20% and 46%, and 26% and 68%, respectively. These statistics indicate that whilst the RA method may not achieve a significantly lower RMSE for vertical velocities, it does retrieve regions of high vertical velocity much more reliably than the pre-gridded methods.

After presenting the mean RMSE and FSS values for the entire domain in Table 1, we proceed to examine the variation of these statistics with the magnitude of the underlying variables in Fig. 6. Firstly, Fig. 6a confirms that horizontal velocities are retrieved very well with the “perfect” observations in the control experiment ($<$0.5 m s^-1 errors across the entire range of magnitudes). For the other retrieval methods, the sampling effects of the simulated radar observations considerably degrade the accuracy of horizontal winds, with comparatively large RMSEs ($>$3 m s^-1). The lowest horizontal velocity errors actually arise in regions with large velocity magnitudes ($\sim$45 m s^-1), which correspond to the mean environmental flow in the upper levels of the storm (8–12 km, refer to Fig. 4). Horizontal velocity errors are lower here due to the uniformity of the flow, with limited storm-induced perturbations and smaller vertical velocities. In contrast, regions with low horizontal velocities ($<$10 m s^-1) occur only very close to the surface, or in strongly dynamic regions of the storm, which oppose the broad-scale environmental flow (such as within the mesocyclone, refer to Fig. 3). Both of these regions are poorly resolved in the pre-gridded retrievals, leading to large errors at low velocity magnitudes (especially for the 2D Cressman method, see Fig. 6a). The 3D Cressman method also predictably shows larger errors for the strongest horizontal velocities, as the over-smoothed radial velocities cannot reproduce the highest magnitude winds. Generally speaking, the RA method produces the lowest horizontal velocity errors and has the most consistent accuracy across the range of horizontal velocity magnitudes.

The difference in error magnitudes between the control experiment and the other retrieval methods is less pronounced for vertical velocities (compare the red lines in Figs. 6a and 6b). This reflects the comparative difficulty of retrieving accurate vertical velocities with low-elevation radial velocities, even with “perfect” radar observations. Errors are lowest for all methods in areas with small vertical velocities ($\sim$0 m s^-1), and generally increase at higher velocity magnitudes. As in the horizontal velocities, 3D Cressman retrievals most poorly resolve the highest wind magnitudes due to over-smoothing of important high-frequency information, and the RA method produces the lowest errors across most vertical velocity magnitudes. The RA method is also able to retrieve considerably larger vertical velocity magnitudes (both updrafts and downdrafts) relative to the pre-gridded techniques, however, these retrieved maximum values are still considerably underestimated. We also examine how FSS values vary with respect to the thresholds used to delineate the updraft and downdraft regions. We test thresholds between the 60th and 98th percentile of the column-maximum vertical velocity field to compute updraft FSSs, and between the 2nd and 40th percentile of the equivalent column-minimum field for downdrafts. Figs. 6c and 6d show the qualitative FSS results presented in Table 1 (which used the 5th and 95th percentile values) are not sensitive to this threshold choice, as the RA method consistently outperforms the pre-gridded methods over the range of thresholds. This is particularly evident for the updraft and downdraft thresholds at the 98th and 2nd percentiles, confirming the RA method is better suited to retrieving regions with the most intense vertical velocities. These ensemble findings concur with the qualitative observations of vertical velocities made based on Figs. 3 and 4 in Section 33.1.

We also analyze how each of the diagnostics presented in Table 1 vary with altitude in Fig. 7. Consistent with Fig. 6, the control experiment exhibits very small horizontal velocity errors throughout the depth of the domain, but most notably at the lowest grid point (500 m altitude). This indicates that the impermeability condition imposed at the lowest grid point, as opposed to the ground surface, does not adversely affect the retrieval in this case (a more sophisticated implementation considering the surface elevation should be used in areas with complex terrain; Chong and Cosma, 2000; Liou et al., 2012). All three radar sampled methods exhibit large horizontal velocity errors ($>$5 m s^-1) at the surface due to the missing data below the lowest radar tilt elevation (0.5^∘), especially the 2D Cressman method (refer to Fig. 2g). Recall that Table 1 showed the majority of RMSE reductions in the RA method were attained through improvements in the horizontal velocity components. Fig. 7a illustrates the origin of these improvements. The RA method produces roughly consistent RMSE scores throughout the depth of the storm ($\sim$4 m s^-1), however, both the pre-gridded methods show considerable deficiencies at different height levels. Firstly, the 3D Cressman method produces large horizontal velocity errors in the mid-levels of the storm ($\sim$4 km) due to a combination of not resolving the small-scale dynamics such as the mesocyclone within the core of the storm (e.g., region 1 in Fig. 3b), and boundary effects due to missing low level data within the forward flank of the storm (e.g., region 4 in Fig. 4b, also refer to Appendix B). Secondly, the 2D Cressman method exhibits large errors in the upper-levels ($>$ 10 km), where the linear interpolation scheme used to fill large data gaps between PPI scans propagates observational noise and creates first-order discontinuities in the underlying radial velocity data (Trapp and Doswell, 2000; Askelson et al., 2000; B22). The ability of the RA method to mitigate these error sources results in overall lower RMSE scores.

Errors in vertical velocities for all four methods (including the control) increase steadily to a height of roughly 4 km, coinciding with the lower portion of the strong ($>$30 m s^-1) updraft in the mid-levels (refer to Fig. 4c). Above 4 km, the errors in the control experiment stay roughly constant at around 2 m s^-1, whereas the radar sampled methods are significantly degraded above 8 km. We attribute this poor performance to the operational radar scanning pattern used in our experiments, which results in large elevation gaps that fundamentally limit the accuracy achievable in the upper levels. As noted in Table 1, the average improvement in vertical velocities gained by the RA method is modest throughout the depth of the storm, with 2D Cressman retrievals actually producing the lowest mean RMSE between 5–9 km in altitude. Again, we note that this result contrasts with the substantial qualitative improvements in vertical velocity observed for the RA method in 3 and 4, and speculate that the discrepancy may be due to the RMSE scores being biased in favor of pixel-to-pixel comparisons of smoother analysis fields. To verify this assertion, we recompute updraft and downdraft FSSs for the 5th and 95th percentile thresholds at each horizontal slice of the analysis grid. Figs. 7c and 7d show that the RA method exhibits modest improvements in FSS compared to the 2D Cressman method between 5–9 km in altitude, despite having slightly higher RMSE scores over the same altitude range. This may indicate a slight spatial displacement in vertical velocities, or that errors from vertical velocities below the FSS thresholds contribute significantly to the RMSE in this region. Furthermore, Fig. 7d illustrates particularly low FSS values and high ensemble spread for the 3D Cressman method, further emphasizing its inability to accurately capture highly dynamic regions of the storm. Overall, Fig. 7 shows that the RA method produces lower error scores, and more reliably retrieves updraft and downdraft regions throughout the depth of the storm.

While the error statistics presented thus far are useful for understanding the error distribution over the entire domain (Table 1), for various magnitudes (Fig. 6), or over a range of altitudes (Fig. 7), some analysts may prioritize accurately identifying the maximum values of certain quantities, such as the maximum horizontal velocity in the low-levels for forecasting severe wind hazards, or the maximum updraft velocity for diagnosing storm severity or verifying/informing convection parameterizations in numerical weather prediction models. These maximum value statistics are provided in Fig. 8. As expected, the pre-gridded retrieval methods considerably underestimate the true horizontal and vertical velocity maxima, due to over-smoothing of important radial velocity information. Encouragingly, the RA method is able to accurately reproduce maximum vertical velocity values below 8 km, before considerably underestimating those farther aloft (at least partly due to the aforementioned observational gaps at high altitude). However, the poor performance of the control experiment at the domain top and bottom in Fig 8b also indicates some of these errors may be partly attributed to the ill-posed vertical boundary conditions in the mass continuity equation.

Finally, the over-smoothing in pre-gridded methods is especially apparent in the derived maximum vorticity and divergence fields in Figs. 8c and 8d, and the qualitative effects of underestimating these dynamic quantities were observed in Figs. 5b and 5c, where simulated trajectories appeared unrealistically uniform or laminar. The RA method considerably outperforms these methods, while also underestimating the true vorticity and divergence maxima throughout the depth of the storm. The most notable improvements gained by the RA method occur in the lower levels in Fig. 4f, where the maximum convergence near the surface is retrieved accurately. This confirms the qualitative findings from the simulated radial velocity data in Fig. 2 and simulated trajectories in Fig. 5, which both suggested the superior characterization of low-level convergence and inflow in the RA method. The accurate characterization of low-level storm inflow is crucial for resolving storm dynamics (Coffer et al., 2022), and this factor likely contributes to the qualitative and quantitative improvements observed from the RA method.

In order to investigate the applicability of the OSSE findings to real cases, we now introduce a supercell case study from November 27, 2014 in Brisbane, Australia. This case was selected as an ideal candidate for verifying the ability to retrieve strong vertical motions, which are known to have occurred based on prior studies (e.g., Parackal et al., 2015; Soderholm et al., 2017). The storm tracked northward through the Brisbane metropolitan area between 0200–0700 UTC, resulting in over AUD 1.5 billion in insured losses due to giant hail ($\sim$70 mm), severe straight-line winds, and localized flooding (Insurance Council of Australia, 2017). The storm was observed by two S-band, Doppler radars: the CP2 research radar and the operational Mt. Stapylton radar (hereafter MS), located to the southwest and southeast of Brisbane, respectively. Fig. 9 shows the analysis domain for our dual-Doppler 3D wind retrieval, centered on the leading edge of the storm at (-27.4^∘S, 153.05^∘E), situated almost entirely within the dual-Doppler lobes ($\sim$40 km between radars). Quality controlled radar volumes are sourced from the Australian Unified Radar Archive (Soderholm et al., 2022), containing nine and fourteen constant elevation sweeps for CP2 and MS, respectively⁹⁹9Exact elevations are $\theta$ = (0.9, 1.7, 2.4, 3.2, 4.7, 6.5, 9.1, 12.8, 17.8)^∘ for CP2, and $\theta$ = (0.5, 0.9, 1.3, 1.8, 2.4, 3.1, 4.2, 5.6, 7.4, 10.0, 13.3, 17.9, 23.9, 32.0)^∘ for MS. Half-power beamwidths/range gates are 0.93^∘/150 m and 1.0^∘/250 m for CP2 and MS, respectively..

In order to process real data, two minor methodological adjustments are made to the partially idealized OSSE retrieval methodology. First, we apply a reflectivity-based correction to account for the effects of hydrometeor terminal velocities in radial velocity measurements. We utilize the terminal velocity correction from the PyDDA package (Jackson et al., 2020), which varies from pure rain (Joss and Waldvogel, 1970) to hail (Conway and Zrnić, 1993) based on reflectivity thresholds. Second, we employ a simple horizontal advection correction to account for storm motion during the finite sampling period of radar volumes (e.g., Shapiro et al., 2009). The average horizontal advection velocity was calculated by taking the mean optical flow velocity (PySTEPS Lucas–Kanade implementation, Pulkkinen et al., 2019) for constant-altitude reflectivity slices between 1–5 km. Reflectivity grids were interpolated using the variational method outlined recently by Brook et al. (2022). Observations from both radars were then shifted to the mid-point of the radar scanning time (06:44 UTC) according to this average advection velocity prior to gridding or ingestion into the RA retrieval method. The analysis grid spans 60 and 50 km in the $x$ and $y$ axes, respectively, with a 500 m grid spacing. The grid is also spaced at 500 m increments in the $z$ axis, but extends from 200–15200 m ASL.

Fig. 10 presents retrieved wind fields for the 2014 Brisbane hailstorm for each of the retrieval strategies discussed in this study. Horizontal streamlines in Figs. 10a–c are broadly similar across each method, resolving the southeasterly change to the rear of the storm, and the strong mesocyclone circulation around highlight region 1 (as in Soderholm et al., 2017). Similar to the OSSE experiments in Section 33.1, horizontal and vertical velocities appear visually smoother in the pre-gridded retrievals (Figs. 10b and 10c), whereas small-scale structures in the vertical velocity fields, especially in the southwest of the domain, are resolved in the RA method (Fig. 10a). The largest vertical velocity differences between methods occur to the north of region 1 in this case. All three methods capture an updraft at the leading edge of the storm, likely a result of dynamic lifting ahead of the advancing gust front. However, the leading-edge updraft appears spuriously large and intense for 2D Cressman method when compared to the main rotating updraft in region 1. This is likely a result of boundary effects caused by data scarcity at the storm’s leading edge in the 2D Cressman method (refer to Appendix B), which under-represents the spatial extent of radar data, particularly in the low-levels (e.g., region 4 in Fig. 2g).

The overall position and southwest–northeast orientation of the updraft at $Z=4$ km in highlight region 1 is retrieved by all retrieval methods. The shape of the inflow notch (indicated by the reflectivity contour in region 1) aligns well with the mesocyclone position and updraft orientation. Maximum updraft velocities for the RA method in Fig. 10a are much larger than that in the pre-gridded retrievals. Quantitatively, the maximum updraft speed retrieved by the RA method was 52.5 m s^-1, compared to 30.8 and 30.3 m s^-1 for the 3D and 2D Cressman methods, respectively. While the veracity of the increased updraft velocities in the RA method cannot be directly confirmed in this real data experiment, the retrieved maximum value of $\sim$50 m s^-1 is closer to the mean values found in supercell simulations (Peters et al., 2020), and the OSSE experiments suggest it is likely a more accurate estimate of the true updraft speed. Further qualitative evidence of considerable updraft strength for this case is given by the large, elevated hail core (¿65 dBZ, $\sim$0 Z${}_{\text{DR}}$), and bounded weak echo region well-above the freezing level (Soderholm et al., 2017; Brook et al., 2021).

Figs. 10d–f corroborate these findings for a vertical cross section through the main updraft. The pre-gridded retrievals adequately estimate the position and shape of the updraft (as evidenced by the position of the bounded weak echo region in region 2 for Figs. 10d–f), but underestimate vertical velocities relative to the RA method. Notably, the strong rear flank downdraft (RFD) signature shown in region 2 for the RA method is also entirely absent in the pre-gridded retrievals. This outflow feature is well documented for this event (Parackal et al., 2015; Soderholm et al., 2017; Brook et al., 2021), resulting in severe wind gusts and strong (¿10 km) hail advection toward the west. The RA method’s ability to resolve the RFD also leads to stronger divergence at lowest grid level, resulting in 40.2 m s^-1 maximum horizontal winds in the outflow region (west of region 2 in Figs. 10d–f), compared with 24.0 and 30.3 m s^-1 for the 3D and 2D Cressman methods, respectively. By comparison, a 39.2 m s^-1 wind gust was measured roughly 10 minutes prior to the analysis time at Archerfield Airport (refer to Fig. 9), along with widespread $\sim$20 m s^-1 gusts throughout Brisbane’s western suburbs (Parackal et al., 2015; Soderholm et al., 2017). While the difficulties associated with comparing retrieved winds aloft to wind measurements at the surface preclude a more quantitative comparison¹⁰¹⁰10Retrieved winds at the lowest grid level ($\sim$170 m AGL) are calculated using a free slip, impermeable, constant elevation boundary condition, and do not resolve the small-scale interactions (i.e., due to topography and surface roughness) that strongly influence wind gust measurements at the surface., the strength and position of these surface measurements provide qualitative support for the RFD signature retrieved by the RA method.

Overall, the results for the Brisbane supercell case study support the findings of the OSSEs, suggesting that the advances made in this study may be realized in practice. However, the practical implementation of the RA method also warrants a brief discussion on the relative computational efficiency of each retrieval method. The introduction of an additional interpolation step in the observational operator for the RA method (i.e., $\mathcal{P}\mathcal{C}$ rather than $\mathcal{P}_{g}$, refer to Section 22.2) roughly doubles the computational cost of the operator relative to the pre-gridded methods. The observational constraint is only one of four otherwise identical constraints in the cost function (Equation 4), however, it must be performed twice at each iteration of the optimization procedure (forward and adjoint). In our implementation, the overall execution time of the RA method is thus $\sim$2.3 times slower than the pre-gridded methods for this case study (15.5 vs. 6.8 seconds, respectively). However, when considering the additional costs involved in pre-gridding the observations (4.4 and 8.0 seconds for 3D and 2D Cressman, respectively), the RA method introduces a relatively modest computational overhead.