Objective identification of pressure wave events from networks of 1-Hz, high-precision sensors

. Mesoscale pressure waves including atmospheric gravity waves, outflow and frontal passages, and wake lows are outputs of and can potentially modify clouds and precipitation. A wavelet-based method for identifying and tracking these types of wave signals in time series data from networks of low-cost, high-precision (0.8-Pa noise floor, 1-Hz recording frequency) pressure sensors is demonstrated. Strong wavelet signals are identified using a wave period-dependent (i.e., frequency-dependent) threshold, then those signals are extracted by inverting the wavelet transform. Wave periods between 1 minute and 5 120 minutes were analyzed, a range which would include several types of mesoscale disturbances in the troposphere. After extracting the signals from a network of pressure sensors, the cross-correlation function is used to estimate the time difference between the wave passage at each pressure sensor. From those time differences, the wave phase velocity vector is calculated using a least-squares fit. If the fitting error is sufficiently small (thresholds of RMSE < 90 s and NRMSE < 0.1 were used), then a wave event is considered robust and trackable. 10

Section 2 of this paper describes the pressure sensors used in this study and the data they provide.Section 3 describes the methodology for objectively identifying pressure waves from the pressure time traces.Section 4 provides five examples of events captured by the wave identification method.Finally, a summary and avenues for future work are discussed in Sect. 5.
To examine the properties of gravity waves which are detectable by these pressure sensors, we consider an internal gravity wave occurring in an environment with constant background wind u 0 The relationship between the pressure perturbation p ′ and the horizontal velocity perturbation u ′ associated with the wave is described by Nappo (2013): where ρ 0 is the environmental air density and c is the phase speed of the gravity wave.From Eq. 1, the maximum pressure perturbation p ′ max can therefore be related to the maximum horizontal velocity perturbation u ′ max by: Figure 3 shows p ′ max values according to Eq. 2 at an air density of 1.225 kg m −3 (standard air density at sea level; American Meteorological Society, 2022) for u ′ max and |c − u 0 | values up to 15 m s −1 .In order to :::::: reliably : detect a wave event :::: signal, the amplitude likely needs to be at least an order of magnitude greater :::::::::: substantially ::::: larger : than the noise floor, meaning gravity waves with p ′ max of at least 0.08 should be detectable using the pressure sensors in this study.

Operational Weather Observations
For context, we compare extracted wave signals with available operational weather observations.
We  Miller et al. (2022), by calculating the difference in radial velocity from successive scans, converting those differences to a binary (positive/negative) field, and filtering out small objects in that binary field.

Methods
The methods outlined here for identifying wave events in the pressure time traces are adapted from the techniques used by Grivet-Talocia and Einaudi (1998) and Grivet-Talocia et al. (1999).The method uses wavelet transforms to identify wave events in time-wave period (or, equivalently, time-frequency) space.Wavelet transforms are preferable to Fourier transforms for the purpose of identifying transient waves which are localized in time (Torrence and Compo, 1998).To illustrate the stepby-step procedure, an example corresponding to an gravity wave event on 23 February 2023 in the Toronto pressure network is described in detail.

Identifying wave events in a single sensor
The full pressure time series for a gravity wave event on 23 February 2023 captured by sensor 25 in Toronto is shown in Fig. 4a.As an initial pre-processing step, 10-second samples of pressure (i.e., averages of 10 pressure measurements) are used to smooth out noise and pressure perturbations due to high frequency turbulent eddies in the data (Fig. 4b).Hereafter, time series labeled as total pressure are the 10-second sub samples of the original pressure measurements.
A wavelet transform W of a finite energy signal f (t) (a pressure time series in this study) can be defined as (Grivet-Talocia and Einaudi, 1998, their Eq. 1): where a is the scale (related to the wave period), and b shifts the wavelet in time (t).ψ * represents the mother wavelet.An analytic Morse wavelet was used (e.g.Olhede and Walden, 2002;Lilly and Olhede, 2012) via the cwt function within Matlab (Lilly, 2021).In this study, W (b, a) will always refer to the wavelet transform of a pressure time series.The resulting wavelet transform is an array of complex values in time-scale space.The absolute value of the wavelet transform |W (b, a)| can be considered the wavelet power at a given time and scale.Figure 4c shows the wavelet power associated with the wave event at sensor 25 on 23 February 2023.In this study, wave periods between 1 minute and 120 minutes were analyzed corresponding to expected periods for mesoscale disturbances.
To objectively identify wave event centers according to wavelet power, a scale-dependent (i.e., wave period-dependent) threshold function A(a) is defined as the mean wavelet power across all available data for the sensor network by scale, multiplied by a constant K: A scale-dependent threshold K is necessary because the 'background' wavelet power for a pressure time series generally increases with scale (e.g., Canavero and Einaudi, 1987).Grivet-Talocia and Einaudi (1998) and Grivet-Talocia et al. (1999) used 2 as an appropriate value for K. Lower values of K lead to more wave events being detected, which can include potential artifacts in the pressure time trace.For the present study, a K value of 10 was used to ensure that only the strongest wave signals were identified (solid contour in Fig. 4d).::: This ::::::::: threshold ::: can ::: be :::::::: adjusted, :::: and ::::::: different ::::::::::: applications :::: may ::::::: warrant ::::::: different ::::: values :: of ::: K. : The mean wavelet power as a function of wave period is shown for each regional sensor network and for all networks combined in Fig. 5. Event centers were identified as local maxima in wavelet amplitude which exceed A(a), which are located at (b max , a max ).In Fig. 4, an event center is located within the solid contour.From the identified event centers, the first iterations of event regions (Ω ′ ) were identified in time-scale space as connected regions where the wavelet power exceeds K 2 ⟨|W (b, a)|⟩ b , i.e., half of the event center threshold.In Fig. 4, Ω ′ is represented by the region within a dashed contour which contains a solid contour.
The watershed transform (Meyer, 1994) was used to refine Ω ′ .Watersheds (i.e., catchment basins) were identified in the negative wavelet power array −|W (b, a)|.Any watersheds within Ω ′ whose period range was entirely outside the period range of the watershed containing the event center were removed from the event region Ω ′ .This step was included to correct cases where multiple "peaks" in wavelet power were present within Ω ′ at different wave periods, with a "valley" in wavelet power in between where the wavelet power still exceeded K 2 ⟨|W (b, a)|⟩ b , which likely represented distinct wave modes and should be considered separate wave events.
Ω ′ was extended to define the final event region Ω for each event, first by taking the bounding box of Ω ′ , then by extending the bounding box along the time axis in both directions until it reaches a local minimum in the wavelet amplitude to ensure that the entire signal of interest is contained in the event region.This could result in overlapping event regions.Figure 4d shows the wavelet power normalized by the mean wavelet power by scale (|W | / ⟨|W |⟩ b ) for the 23 February 2023 example in sensor 25, with the outline of the event region overlaid with the magenta box.

Matching corresponding wave events between multiple sensors
Once wave events were identified for each sensor individually, the following steps were taken to identify coherent wave events across multiple sensors.For this purpose, the terms primary sensor (denoted by i) and secondary sensor (denoted by j) will be used to describe a pair of sensors for which events are identified and paired together.
The process of matching events between sensors described above was repeated for each possible combination (within a sensor network) of primary and secondary sensors in order to obtain the full set of lag times between each pair of sensors which captured each event.In other words, N 2 pairs of sensors, order-dependent, were analyzed, where N is the number of sensors in the network with data at a given time.Then, each event in each sensor was assigned an ID based on which other sensors had a matching event in order to track events across 3 or more sensors.This process required iterating through each sensor in a network.Each event in the first sensor was assigned a new (i.e., arbitrary) ID.For each subsequent sensor s c , events with no two-way matches in any prior sensor were also given new IDs.If there were two-way matches with an event in one or more prior sensor(s), the event in the sensor s c would share the ID assigned to the matched event in the prior sensor.If there were multiple prior sensors with matched events, and those events had different IDs, the ID associated with the higher maximized cross-correlation between the event traces was assigned to the event in sensor s c .If this process results in multiple events in sensor s c sharing the same event ID D, the event in sensor s c associated with the highest maximized cross-correlation with any one prior sensor for an event with event ID D is assigned event ID D, and the previously outlined steps are repeated for the other event(s) in sensor s c , with event ID D and associated sensors excluded.
The result of this process is a set of events with associated ID numbers for every sensor in the network.For a single sensor, each event has a unique ID.For each ID number that appeared in at least three sensors, the wave phase velocity vector was calculated using the set of lag times between each pair of sensors which captured the event.

Estimating wave phase velocity vector
Once sets of matched events were identified, the wave propagation velocities (two-dimensional vectors) could be estimated for events which occurred in three or more sensors.It is hypothesized for each wave event that a plane wave crosses the sensor network with slowness vector s = (s x , s y ), where s x and s y are the inverses of the x-and y-components of the wave propagation vector (in s m −1 ), respectively.s can be solved for from the following equation (Del Pezzo and Giudicepietro, 2002): where t is the column vector of the ∆t opt values for each possible pair of sensors which captured the event, and ∆x is the two-column matrix of the x-and y-components of the distance vector between each pair of sensors which captured the event.t and ∆x each have N s (N s − 1)/2 rows, where N s is the number of sensors which captured the event.Equation 6 can be considered an overdetermined system of N s (N s − 1)/2 linear equations, as long as N s ≥ 3, and is solved for s by a least squares approach represented by: where superscript T indicates the transpose of a matrix (Del Pezzo and Giudicepietro, 2002).
Once s x and s y are solved for, they can be inverted to obtain the wave phase velocity components, c x and c y , respectively.
Additionally, we require that wave events be captured by at least 4 sensors to be considered trackable.While the slowness vector and corresponding error metrics can be calculated for events captured by only 3 sensors, the calculation is less constrained because there are only 3 delay times δt opt in the calculation (compared to 6 delay times for events captured by 4 sensors, 10 delay times for events captured by 5 sensors, etc.).The result is that events captured by only 3 sensors can have small RMSE and NRMSE by chance much more easily than events captured by 4 or more sensors.For each robust and trackable wave event, the mean amplitude was calculated by averaging the difference between the maximum and minimum values in the extracted event trace for each sensor which captured the event.The center period for the event was calculated as the mean of the wave period corresponding to b max for each sensor which captured the event.

4 February 2022: Outflow pressure jump and subsequent oscillations over Long Island
Between 1730 and 1900 UTC on 4 February 2022, an event with amplitude of roughly 1.8 hPa was detected by the pressure sensor network ::: five ::::::: pressure ::::::: sensors : in the New York City metro area and Long Island.This event was a positive jump in pressure followed, to varying degrees in each sensor, by weak oscillations in the pressure trace (Fig. 11).Prior to the jump in pressure, there had been a decreasing trend in the pressure traces for several hours.At the same time, WSR-88D weather radar data from Upton, NY, showed widespread precipitation echo over Long Island.A wave feature was apparent in the Doppler radial velocity data, which could also be identified following the methods of Miller et al. (2022) (Fig. 12 and Video Supplement Animation-Figure -S02).This wave event had a phase speed of 21.1 m s −1 and direction of 118.2 • clockwise from north (i.e., southeastward).The values are consistent with the radar-detected Doppler velocity wave feature (Video Supplement Animation-Figure -S02).
Operational one-minute Automated Surface Observing System (ASOS) data (Fig. 13a) also recorded a jump in the surface pressure of nearly 2 hPahPa.Near the time of this jump, there was also a peak in the wind speed and gusts, along with a brief shift in the wind direction from north-northeasterly to north-northwesterly (Fig. 13b).These features, along with the modest decrease in the temperature (Fig. 13a), are consistent with a convective outflow boundary (i.e., gust front).A "fine line" can be seen in WSR-88D reflectivity data at roughly the same location as the wave, which further suggests that a convective outflow was responsible for the pressure rise (Fig. 12 and Video Supplement Animation-Figure -S02).

15 November 2020: Cold front passage over Toronto
A robust and trackable wave event was detected by 5 ::: five : pressure sensors in Toronto coincident with a cold front passage at roughly 2000 UTC on 15 November 2020.The pressure steadily dropped in the hours leading up to the frontal passage before abruptly rising 1-2 hPa as the cold air mass arrived.The pressure then dropped roughly 1 hPa about 30 min min later.Some sensors recorded oscillations in the pressure trace embedded within the gradual pressure rise in the proceeding hours (Fig. 14).
One-minute ASOS data from Buffalo, NY, also captured the pressure jump at roughly the same time as the temperature and dew point drop indicating the cold front passage (Fig. 15).
The pressure wave event for the cold front passage had an estimated phase speed of 27.5 m s −1 at 65 • (i.e., to the eastnortheast), a mean amplitude of 1.8 hPa, and a center wave period of 00:02:08.WSR-88D radar data from Buffalo, NY, show a narrow band of high reflectivity and a Doppler velocity wave (identified following the methods in Miller et al., 2022) associated with the cold front advancing over Toronto at roughly 2000 UTC at a speed and direction consistent with the pressure wave (Video Supplement Animation-Figure -S03).
4.5 14 September 2021: Wake low associated with a mesocale convective system Between 0300 and 0400 UTC on 14 September 2021, a pressure drop of roughly 5 hPa and subsequent recovery occurred at 4 ::: four : of the pressure sensors in New York City metro area and Long Island network (Fig. 16).This was detected as a wave event with an estimated propagation speed of 20.6 m s −1 and propagation direction of 67.5 • (i.e., to the east-northeast).ASOS data from Islip, NY (KISP; Fig. 17), and other stations in the area (not shown) also recorded the pressure minimum.This wave event occurred near the time of a mesoscale convective system (MCS) passage over Long Island as indicated by reflectivity data from the WSR-88D radar in Upton, NY (Fig. 18 and Video Supplement Animation-Figure -S04).In addition to the precipitation echo associated with the MCS translating from northwest to southeast, there was also a stationary region of weak echo with low dual-polarization correlation coefficient (shown in greyscale following Tomkins et al., 2022) in the vicinity of the radar.The stationary weak echo was likely non-meteorological and due to either birds or insects.The precipitation echo appears to be entirely past KISP by 0342 UTC (Fig. 18c and Video Supplement Animation-Figure -S04), which is roughly the same time as the minimum pressure at KISP (Fig. 17a).
This pressure minimum appears to be consistent with a wake low, associated with subsidence heating in the rear inflow jet (Markowski and Richardson, 2010;Johnson and Hamilton, 1988).The subsidence heating does not necessarily lead to warming at the surface, which was not observed in the ASOS data (Fig. 17a), but decreased air density aloft due to warming will still lead to a surface pressure decrease.Markowski and Richardson (2010) also note that a property of wake lows associated with a translating squall line is that the center of convergence due to the wake low does not perfectly align with the center of the wake low.Rather, the convergence center slightly lags behind the wake low center.In the 14 September 2021 example, the ASOS time series data have a wind speed minimum co-occurring with a shift in the wind direction from near 100 • (eastsoutheasterly) to near 280 • (west-northwesterly) which can be interpreted as the convergence maximum (Fig. 17b).This convergence maximum occurs slightly after the pressure minimum associated with the wake low (Fig. 17a), consistent with the Markowski and Richardson (2010) description.
Code and data availability.Data: The pressure time series data used throughout this publication can be found at https://doi.org/10.5281/zenodo.8136536(Miller and Allen, 2023) Lindzen and Tung (1976) and Koch and O'Handley (1997).

:
:::::: represent : a case where where the origin of the waves is known and the phase speed is constrained to the speed of sound :::::: known ::as : a ::::::: function ::: of :: air ::::::::::: temperature(Amores et al., 2022).A gravity wave train which passed over Toronto on 25 Feb ::::::: February 2022 occurred coincident with a surface cyclone 100 km distant but in local conditions of sparse radar echo.A wave event on 4 Feb ::::::: February : 2022 is a case associated with an outflow boundary clearly captured by Doppler radar data from the WSR-88D radar located in Upton, NY.In another example, waves coincided with a cold front which passed over Toronto on 15 Nov ::::::::

Figure 2 .
Figure 2. Locations of pressure sensor networks.(a) US northeast regional map with the locations of the Toronto, New York City and Long Island, and Raleigh networks indicated.Detailed maps of sensor locations in (b) Toronto, (c) New York and Long Island, and (d) Raleigh.

Figure 3 .
Figure 3.The maximum pressure perturbation p ′ max (hPahPa) contoured and colored as a function of the maximum velocity perturbation u ′ max (m s −1 ) and absolute value of the difference between the wave phase speed and background wind speed |c − u0| (m s −1 ) at a density ρ0 of 1.225 kg m −3 , according to Eq. 2. In (a), u ′ max and |c − u0| up to 15 m s −1 are shown.In (b), u ′ max and |c − u0| up to 3 m s −1 are shown.The color scales differ in (a) and (b).

Figure 4 .
Figure 4. Steps in process of identifying an event corresponding to a gravity wave train passage on 23 Feb 2023 in sensor 25.(a) Original 1-second pressure time series.(b) 10-second moving average of the pressure time series, with every 10th point kept.(c) Wavelet power corresponding to the pressure time series in (b).(d) Wavelet power normalized by the mean for each wave period corresponding to the pressure trace, with contours for values of 5 (dashed) and 10 (solid).The event region corresponding to the wave described in the text is outlined in magenta in (d).(e) Time series of the extracted wave event.Wave periods in (c) and (d) are shown on a logarithmic scale.All times are UTC.

Figure 5 .
Figure 5. Mean wavelet power as a function of wave period across the entire data set (purple curve), and for the individual sensor networks around Toronto, ON, Canada (blue), New York City and Long Island, NY (red), and Raleigh, NC (yellow).

Figure 7 .
Figure 7. Extracted waveforms corresponding to the gravity wave train on 23 Feb, 2023 for sensors 25 (a), 04 (b), 23 (c), 24 (d), and 34 (e).In (b), (c), (d), and (e), brown lines show the extracted wave event with no time shift, and black lines show the extracted wave event shifted in time according to the peak cross-correlation with the extracted wave event time series for the event in sensor 25 (time shift and correlation coefficients are shown in subplot titles).All times are UTC.

Figure 8 .
Figure 8. Pressure traces and extracted waveforms for two events associated with the shockwave :::: Lamb ::::: waves caused by the Hunga Tonga-Hunga Ta'apai eruption on 15 Jan 2022.Extracted waveforms (black lines) are overlaid on the total pressure time series (blue lines).(a) the Outbound (ii) wave event propagating from Tonga toward the antipode location in Algeria in Sensor 02 (Table 1, column 3).(b) the Rebound

Figure 9 .
Figure 9. Pressure traces and extracted waveforms for the 25 Feb 2022 wave event in sensors (a) 04, (b) 24, (c) 26, and (d) 34.Extracted waveforms (black lines) are overlaid on the total pressure time series (blue lines).All times are UTC.Cross-correlations and lag times are indicated relative to sensor 04.Cross-correlations are computed for each variation of pairs of sensors (not shown).

Figure 10 .
Figure 10.Upper air sounding data from Buffalo, NY, valid for 1200 UTC 25 Feb 2022.(a) Dry bulb temperature (blue), dew point (orange), wet bulb temperature (green), and frost point (red) profiles (all in • C).(b) Equivalent saturation potential temperature (θ * e ) profile (black line, K) overlaid on the vertical gradient in θ * e (K km −1 ).Positive values (blue) of the vertical gradient in θ * e indicate absolute stability, while negative values (red) indicate conditional or absolute instability.(c) Horizontal wind profile (barbs, kts; colored according to wind speed).Annotation indicates the vertical extents of a ducting layer and a trapping layer according to the gravity wave duct criteria described byLindzen and Tung (1976) andKoch and O'Handley (1997).

Figure 11 .
Figure 11.As in Fig. 9, but for the 04 Feb 2022 wave event in sensors (a) 21, (b) 22, (c) 11, (d) 14, and (e) 20, with cross-correlations and lag times indicated relative to sensors 21.Extracted waveforms (black lines) are overlaid on the total pressure time series (blue lines).

Figure 12 .
Figure 12.(a) Reflectivity and (b) Doppler velocity wave detection for NWS WSR-88D radar data from Upton, NY, at 0.5 • tilt, at 1920 UTC on 4 Feb 2022.In (a), reflectivity values are shown in greyscale where there is likely enhancement due to melting (Tomkins et al., 2022).Filled blue circles indicate locations of pressure sensors which captured the wave event described in Sect.4.3, and unfilled blue circles indicate locations of pressure sensors which did not capture the wave event.Filled green circle indicates location of Islip, NY, ASOS station (KISP).An animation of this figure showing the time sequence from 1541 to 2130 UTC is in Video Supplement Animation-Figure-S02.

Figure 13 .
Figure 13.Time series of one-minute ASOS data from Islip, NY (KISP), on 4 February 2022.(a) Temperature (purple) dew point (red), and pressure (blue).(b) Wind speed (orange), wind gust speed (red), and wind direction in degrees clockwise from northerly (yellow).Wind direction is not plotted when it changes by more than 180 • in consecutive observations (e.g., when crossing 0 • or 360 • ) or when the wind speed is below 1.5 m s −1 .All times are UTC.

Figure 14 .
Figure 14.As in Fig. 9, but for the 15 Nov 2020 wave event in sensors (a) 34, (b) 04, (c) 24, (d) 25, and (e) 25, with cross-correlations and lag times indicated relative to sensor 34.Extracted waveforms (black lines) are overlaid on the total pressure time series (blue lines).

Figure 15 .
Figure 15.Time series of one-minute ASOS data from Buffalo, NY (KBUF), on 15 November 2020.(a) Temperature (purple) dew point (red), and pressure (blue).(b) Wind speed (orange), wind gust speed (red), and wind direction in degrees clockwise from northerly (yellow).Wind direction is not plotted when it changes by more than 180 • in consecutive observations (e.g., when crossing 0 • or 360 • ) or when the wind speed is below 1.5 m s −1 .All times are UTC.

Figure 16 .
Figure 16.As in Fig. 9, but for the 14 Sep 2021 wave event in sensors (a) 21, (b) 27, (c) 14, and (d) 18, with cross-correlations and lag times indicated relative to sensor 21.Extracted waveforms (black lines) are overlaid on the total pressure time series (blue lines).

Figure 17 .
Figure 17.Time series of one-minute ASOS data from Islip, NY (KISP), on 14 September 2021.(a) Temperature (purple) dew point (red), and pressure (blue).(b) Wind speed (orange), wind gust speed (red), and wind direction in degrees clockwise from northerly (yellow).Wind direction is not plotted when it changes by more than 180 • in consecutive observations (e.g., when crossing 0 • or 360 • ) or when the wind speed is below 1.5 m s −1 .All times are UTC.

Figure 18 .
Figure 18.Maps of radar reflectivity at 0.5 • tilt from the NWS WSR-88D radar in Upton, NY, on 14 September 2021.Reflectivity values in color show meteorological echo and those are shown in greyscale are likely non-meteorological echo such as insects and birds in this case (Tomkins et al., 2022).Filled blue circles indicate the locations of pressure sensors which captured the wave event described in Sect.4.5, and unfilled blue circles indicate the locations of pressure sensors which did not capture the wave event.Filled green circle indicates the location of the Islip, NY, ASOS station (KISP).The sequence of images from (a) 0303 UTC, (b) 0323 UTC, (c) 0342 UTC, and (d) 0402 UTC shows the southeastward movement of region of convective cells > 40 dBZ from closer to further off the southern coast of Long Island.The wake low is inferred to be near the trailing edge of the weaker stratiform precipitation region behind (west of) the convective cells.The minimum pressure at KISP associated with the wake low occurred near the time of the scan shown in (c).An animated version of this figure, with Doppler velocity wave detection, is shown in Video Supplement Animation-Figure-S03.
use Automated Surface Observing Systems (ASOS; NOAA National Centers for Environmental Information, 2021a) data including surface temperature, dew point, wind speed and direction, and additional pressure measurements coincident with wave events.ASOS data are recorded each minute.For wave events detected in New York and Long Island, we examined ASOS . The NWS NEXRAD Level-II data used in Figs. 12 and 18 can be accessed from the National Centers for Environmental Information (NCEI) at https://www.ncei.noaa.gov/products/radar/next-generation-weather-radar(NOAA National Weather Service Radar Operations Center, 1991).The NWS ASOS surface station data used to create Figs. 13 and 17 can be accessed from NCEI at https://www.ncei.noaa.gov/products/land-based-station/automated-surface-weather-observing-systems(NOAA National Centers for Environmental Information, 2021a).The radiosonde data used to create Fig. 10 can be accessed from NCEI at