Detecting Wave Features in Doppler Radial Velocity Radar Observations

. Mesoscale, wave-like perturbations in horizontal air motions in the troposphere (velocity waves) are associated with vertical velocity, temperature, and pressure perturbations that can initiate or enhance precipitation within clouds. The ability to detect velocity waves from horizontal wind information is an important tool for atmospheric research and weather forecasting. This paper presents a method to routinely detect velocity waves using Doppler radial velocity data from a scanning weather radar. The method utilizes the difference ﬁeld between consecutive PPI scans at a given elevation angle. Using the difference 5 between ﬁelds a few minutes apart highlights small scale perturbations associated with waves because the larger scale wind ﬁeld changes more slowly. Image ﬁltering retains larger contiguous velocity bands and discards noise. Wave detection scales are limited by the size of the temporal difference relative to the wave motion and the radar resolution volume size.

not and Section 6 contains our conclusions. There is a video supplement which includes animated versions of all figures in this manuscript.

60
The technique described in this paper is designed for use in the lower troposphere and works well even if the background wind field is complex. By subtracting two consecutively scanned Doppler radial velocity fields, we can remove the more slowly varying background wind and isolate relatively faster moving velocity wave features that have noticeable movement between consecutive scans.
The WSR-88D Level II Doppler radial velocity data are first processed to unfold (dealias) velocities (Helmus and Collis, binary threshold value balances removing noise and detecting the features of interest. The binary threshold value appropriate for any specific data set will vary based on characteristics such as the phenomena of interest, the sampling frequency of the radar, and the noise level of the radar. For the next step, the polar coordinate binary field was interpolated to a 0.5 km Cartesian grid using a nearest neighbor algorithm. A grid resolution of 0.5 km is sufficient for Level II WSR-88D data based on the 0.93°beamwidth and the unambiguous 95 range (137 km in this example). The goal is to oversample the data to preserve the detail native to the original polar field.
As a last step, we filtered out objects smaller than 16 km 2 . We then used the skimage morphology Python package (van der Walt et al., 2014) to remove eight-connected, continuously connected areas of positive wave detection smaller than 16 km 2 . In implementing a minimum size filter, the users are making a choice about the size scale of waves they are interested in. Waves in storms are common across a variety of spatial scales and orientations. Myriad waves often overlap. A given minimum or 100 maximum value for a size filter will emphasize the detection of a certain scale of waves. The appropriate Cartesian grid resolution and area size to filter out is a function of the radar characteristics and intended application. In this example, we focused on the larger velocity waves and filtered out the smaller scale features including most of what appears to be spatial noise in the western half of the domain.
The resulting field (Fig. 2e and associated animation in the Video Supplement) contains several sets of parallel velocity wave 105 features. The most prominent are the set with their long axis aligned from SSW to NNE. The horizontal wavelength of these waves is on the order of 12 to 18 km. Examination of sequences of detected waves shows the wave train propagating toward the NW and spanning the radar domain. Embedded among the longer SSW to NNE waves are collections of shorter wavelength horizontal waves oriented in a variety of directions. In the southeastern quadrant of the radar domain, the detected velocity features are less spatially coherent in part because the radar echo is more sporadic in that region.

110
Within the limitations described below in Section 2.1, we can use this method of taking the difference between successive Doppler radial velocity fields and converting to a binary field to identify velocity wave structures. Examining the binary difference fields as they evolve gives information on the speed and direction wave trains propagate. Measuring the distance between adjacent wave features gives the horizontal wavelength. Examining several different elevation angles and/or computing a vertical cross section based on multiple elevation angles provides information on the spatial coherency of the waves with 115 height and the depth of the waves. The velocity difference technique described here can also be applied to sequential Doppler velocity data from RHIs scanned along the same azimuth to detect possible waves (see Section 4).
Tracking coherent features in the wave detection data across successive volumes will allow the user to calculate the wave propagation direction and speed by measuring the direction and distance of the displacement of these features and accounting for the difference in sampling time. There are several different approaches to performing such an analysis using manual and 120 automated methods. The details of any successful approach will be largely dependent on the specifics of the radar data, the features of interest, and a user's applications. In some situations, the application of non-rigid image registration tools as well as optical flow and related methods for feature tracking have shown promise. The ability to resolve the velocity waves in Doppler radial velocity data is constrained by the spatial resolution of the radar, the temporal sampling interval between consecutive scans, and the temporal scale of changes in the background wind field.
Additionally, there is the intrinsic limitation of radial velocity data in that the radar only measures the component of the wind along azimuths originating at the radar. The WSR-88D radar can detect velocity changes ≥ 0.5 m s -1 (NOAA, 2017). The size of radar sample volumes is a function of the radar antenna beamwidth, the size of the range gates (pulse length), and the 130 range from the radar (Battan, 1960). At minimum, the velocity waves must propagate horizontally at least the equivalent of one radar resolution volume in distance along the direction of wave motion to show as a difference in the PPI radial velocity data. In practice, we have found that the waves need to propagate at least 3 times the horizontal radar resolution size to reduce uncertainty related to beam filling and velocity aliasing and to have a signal clearly distinguishable from the background noise.
For the US WSR-88D network, the radars have a typical half-power beamwidth of 0.93°and the length of the range gates for 135 Doppler velocity are 250 m (NOAA, 2006). At a range of 100 km, the WSR-88D radar resolution volume size is an approximate cylinder 1.6 km wide by 250 m long. The temporal sampling interval varies among different volume coverage patterns (VCPs).
For the US WSR-88D network VCPs, the time between consecutive scans at the same elevation angle is typically between a few minutes for scans designed for severe weather and up to approximately 10 minutes for scans designed for clear air and light rain or snow precipitation (Rauber and Nesbitt, 2018). For a radar resolution volume horizontal distance of 2000 m and a time 140 interval between consecutive scans of 300 seconds, the velocity wave needs to propagate at a minimum of speed of 6.67 m s -1 to traverse a sample volume.
For this wave detection method to work best, the background wind field must change as little as possible between the radial velocity observations used to calculate the difference field. Turbulence must also be low since this creates noise in the radial velocity difference fields. Time scales for storm evolution usually increase as storm spatial scales increase. Large, synoptically 145 driven systems typically evolve more slowly than rapidly changing, convective-scale thunderstorms. If the spatial scale of the background wind change is small, it creates the equivalent of small scale noise that can be mitigated by techniques such as removing contiguous areas in the binary difference field that are small in area and not likely to be waves. If the changes in the background wind field are large in spatial scale, they will mask the presence of velocity waves.
Not every linear velocity band feature detected by this method will correspond to a wave. Wind shifts along lines, including 150 convergence lines, gust fronts, and fronts, will be detected as single velocity bands and will usually move with the prevailing wind. In contrast, propagating wave features will appear as sets of roughly parallel linear velocity bands that move in concert and may or may not move with the same speed or direction as the prevailing wind. Users should keep in mind the meteorological context of the observations so they can make informed judgements about which velocity band features are waves and which are not.

Wave motion
As compared to simple wave detection, there are more stringent conditions to obtain an accurate speed and direction for propagating waves. To obtain the correct wave speed, the wave train must move less than half its wavelength between sample times. A wave train that moves more than half its wavelength can appear to propagate in the wrong direction. A wave train that moves exactly its wavelength in a sample period will appear stationary. The sample time, in seconds, must be less than where t is the sample time in seconds, λ is the wavelength in meters, and V is the wave propagation speed in meters per second.
An NWS radar in VCP-12 completes a volume approximately every 240 s. A wave with a wavelength of 10 km would therefore have to have a propagation speed of less than 20.8 m s -1 to avoid errors in estimating propagation speed and direction. Further discussion of this issue is in Section 3.1.

Application to Idealized Waves and Storm Examples
In order to test and illustrate the velocity wave detection method in a controlled environment, we wrote software to create idealized plane waves and sample them as a radar would (Miller, 2021). This allowed us to verify that we are able to correctly estimate wavelengths, depths, speeds, and directions without the uncertainty and noise found in real-world data. In the idealized examples below, radial velocity data were generated for plane waves with various defined wavelengths, depths, orientations, 170 and amplitudes of the horizontal velocity perturbations. The code was used to compute the wave characteristics in a Cartesian coordinate system as well as to sample the waves along elevation angle PPIs for a set of user-specified radar characteristics, including maximum unambiguous range, number and size of range gates, and beamwidth. Standard refraction and earth curvature were accounted for using the 4/3rds Earth approximation (Doviak and Zrnić, 1993;Rinehart, 2010). All simulated radar points in a PPI were sampled all at the same time as compared to ray-by-ray sequential sampling as a real radar would scan.

Single Wave Example
We constructed an idealized plane wave by prescribing horizontal velocity perturbations as a sinusoidal plane wave. Figure 3 and When the waves are tilted (Fig. 3bd), they appear bowed in PPI radar images. This effect is caused by the height of the radar beam over the surface increasing as a function of range.

185
When the time between samples is too long (Eq. 1), the apparent wave speed will be wrong. Video Supplement movies Animation-MotionStudy1 and Animation-MotionStudy2 illustrate mischaracterizations of wave speed when the temporal sampling is insufficient. As described in Section 2.1.2, if the temporal sampling is equal to the time a wave takes to travel an integer multiple of its wavelength, the wave will appear to not be moving. A wave could also appear to be propagating the incorrect direction if the temporal sampling frequency is too low. The animations feature a single wave with a wavelength of Animation-MotionStudy2, the wave propagates at 40 m s -1 and moves more than half its wavelength between samples but less than its full wavelength. This results is the appearance that the wave is moving in the opposite direction at a slow speed. These effects can also be illustrated by using a fixed wave propagation speed and a variable radar sampling frequency.

Two Waves Example
In real storms, several waves with different wavelengths, orientations, and depths often coexist. Figure 4 and Video Supplement In the 0.5°elevation PPI (Fig. 4), the beam tops the 1000 m wave at about 70 km range. In cases where one is certain the beam is overshooting a wave train, the range at which this occurs can be used to infer the depth of the wave by referencing how the radar beam height increases as a function of range. When two waves with different wavelengths overlap, the longer wavelength wave typically dominates and can obscure the signal from the smaller wavelength plane wave. The obscuring of the 210 shorter wavelength waves is most evident when the orientation of the two wave sets are similar. As the number of superimposed waves increase, the more the wave detection information degrades to noise.
The contrasting wave characteristics in vertical cross section are illustrated in cross sections of the idealized wind field and in the RHIs constructed from the PPI scans in Fig. 4 and Animation-Figure4. The RHIs were constructed as one would from operational radar data, by finding the subset of data along a given azimuth in each elevation angle's PPI in the radar 215 volume scan, defining the radar resolution volume center locations corresponding to slant range coordinates, and determining the range-dependent radar resolution volume sizes. The constructed RHIs are along an azimuth of 90°, which is between the two idealized wave orientations of 290°and 250°azimuth and coincides with the cross sections also shown. The set of elevation angles, the beamwidth, and the gate spacing were chosen to match the NWS VCP-12 (Rauber and Nesbitt, 2018) VCP. The tops of the two prescribed waves are denoted in the RHI panels by dashed, gray lines.

220
The wave detection technique is capable of accurately estimating the tops of wave trains where they are sufficiently resolved.
When using RHIs constructed from sets of PPIs, the wave top estimates are most accurate close to the radar where there is greater effective vertical resolution than at farther ranges. The 1000 m deep wave is more difficult to discern because of the overlap between the two waves in the 0 to 1000 m layer. Examination of sequences of cross sections can help mitigate the visualization problem if the waves do not move in phase with each other. Scanned RHIs, when the radar antenna varies 225 elevation angle while holding a given azimuth, typically have higher vertical resolution than RHIs constructed from PPIs and better enable the resolving of wave features.

Winter Storm Examples
We applied the wave detection method to PPI and RHI radar data from a winter storm on 1 February 2021 that impacted the northeast US. Figure 5 shows radar data from 09:03 UTC. The data are from the NWS KOKX radar located on eastern The KASPR radar has a narrow beamwidth (0.32°) and, for this case, a range resolution of 24 m for 1154 range gates yielding a maximum unambiguous range of 29 km., which yields higher spatial resolution than a WSR-88D radar. Prior to wave detection, the KASPR fields were linearly interpolated to a grid with 10 m grid spacing. For this storm, KASPR RHI 245 scans alternated from north-to-south to south-to-north. Because of this, the waves were calculated using the difference from every other scan so that the temporal spacing between elevations in the scans remain consistent. We used a binary threshold value of -0.4 m s -1 and a minimum area filter of 0.03 km 2 threshold. The wave detection product derived from the high spatial resolution RHI data shows the presence of vertically oriented, coherent, linear wave patterns between 8 and 16 km south from the radar. Examination of the radial velocity difference field suggests the presence of two different modes of waves.

250
structures revealed by the wave detection method propagate mainly along their long axis which is not consistent with motion of the synthetic wave data examined in 3.2.

6 Conclusions
Our method identifies translating bands of wind perturbations (velocity waves) by subtracting successive radial velocity fields.
This method can be used to detect the presence of several types of velocity waves and serve as the foundation for analysis on the role of gravity and Kelvin-Helmholtz waves in storms. Manual or automated measurement of the flagged wave features can be used to determine their wavelength. Tracking the waves in time allows the estimation of the wave propagation direction 290 and speed. Such estimates must be done with care to ensure the spatial and temporal resolution is sufficient to resolve the features of interest and to capture their motion with fidelity (Section 2.1). This technique should be extendable to Doppler radial velocities from other instruments such as cloud radar and lidar. Specific filtering thresholds for such instruments should be selected consistent with instrument specifications and intended application. Processing of radar observations with different dynamic ranges and noise floors than a WSR-88D radar will likely use different threshold values.

295
The wave detection method was tested using sets of idealized planar waves and accurately captures the idealized waves' characteristics within spatial and temporal resolution limits. Examination of winter storms impacting the northeast US reveals the presence of overlapping packets of quasi-planer, monochromatic waves with wavelengths spanning 1 to 40 km. Radar PPIderived cross sections of the waves suggest they have depths of several kilometers. Within the same storm, some parallel sets of waves move with the prevailing wind and others move independently of the prevailing wind. Some sets of waves coincide 300 with banded features of higher reflectivity while other sets move independently of mesoscale reflectivity bands. The nature of interactions among radar-detected velocity waves of several types and reflectivity bands is a topic of ongoing work.
The wave detection output requires informed interpretation because waves are not the only source of linear velocity perturbations. Outflow boundaries, sea-breeze fronts, frontal boundaries, etc., can produce a propagating change in velocity that will be flagged by our method. The user must assess the flagged locations to ensure that they possess characteristics consistent with 305 waves such as presenting as a set of parallel banded features. A creative user can apply various filtering techniques to the data to remove areas that are not likely to be waves. A full discussion of potential approaches is beyond the scope of this paper because it is highly application dependent.
Code and data availability. MATLAB Code for Idealized Wave Examples The matlab-RadarSim library (Miller, 2021) used to create the idealized wave figures and animations is available on GitHub.