Utility of ocean wave parameters in ambient noise prediction^a)

Ambient noise (AN) in the ocean is produced by shipping, biologics, wind, precipitation, sea ice, and other sources. The goal of this study is to improve prediction of the component of ocean AN commonly referred to as “wind noise.” Important early work on this topic was performed by Knudsen (1948) and Wenz (1962). In these works, it was recognized that the sea state (“surface agitation”) has a central role in the creation of this noise, while the underlying mechanisms were imprecisely ascribed to some combination of the action of bubbles in the water and/or the impact of sea spray on the surface. Wilson (1980) argued that the latter was a dominant mechanism, but by the end of the decade, there was general recognition that bubbles play the more important role (e.g., Prosperetti, 1988; Medwin and Daniel, 1990). Prosperetti (1988) outlined three frequency regimes: (1) low frequencies, 100–200 Hz, where bubbles in the water act as amplifiers of “water turbulence noise,” (2) intermediate frequencies, 200 Hz to 1 or 2 kHz, where there is typically a broad maximum in the wind/wave noise, caused by “collective oscillation of bubble clouds,” (3) high frequencies, 1 or 2 kHz to around 10 kHz, where individual bubble oscillation dominates. The last regime is the most well understood; this noise is emitted immediately (e.g., 10 s of ms) (Manasseh , 2006) after the bubbles are formed, and the frequency of the emission is directly related to the bubble radius (Minnaert, 1933; Manasseh , 2006).

It is generally understood that wave breaking is the link between wind speed and bubble generation. Further, Farmer and Vagle (1988) established that the third noise mechanism above is a direct result of wave breaking. However, in the real ocean, wave breaking is not a simple function of wind speed. While wind speed is of primary importance, the frequency and severity of wave breaking is affected by other environmental factors, e.g., see review by Babanin (2011, Chaps. 2 and 6). Consider, as illustration, an idealized case of 15 m/s wind speed over 100 m basin vs the same wind over 100 km basin. The oceanographic difference (specifically, fetch) will result in waves of different height and length, different strength of breaking, different bubble generation, and therefore, different acoustic emissions for these two cases. Current AN models use wind speed, plus water depth [usually binary, i.e., deep vs shallow, as in Wenz (1962)] or receiver depth [as in Hildebrand (2021)] and predict AN as a function of frequency. Our hypothesis is that since the AN is more directly connected to wave breaking than to wind speed, then it should be possible to improve a prediction of wind noise by incorporating oceanographic information, particularly that associated with ocean waves.

The above hypothesis is not new. It was central to the landmark study of Felizardo and Melville (1995) (henceforth denoted “FM95”), who quantified the correlation between AN and several wind/wave parameters for a relatively short (10 day) field experiment, looking at three frequencies, from 4.3 to 14 kHz. They found that the AN correlation with predicted wave dissipation rate is equal to, or slightly higher than, the correlation with wind speed. In the present study, we perform a similar analysis, though here we apply multiple input parameters simultaneously to predict AN, rather than one parameter at a time, as in FM95. Further, we use a larger dataset; the cumulative duration of our dataset is 1291 days, making it more than 100 times larger than that of FM95, and a much broader range of frequencies is studied (127 Hz to 18 kHz vs 4.3–14 kHz). We employ wind and wave information at seven different hydrophone locations. At one location, wave characteristics are taken from observations, while at the other six locations, wave parameters are derived from numerical ocean wave models.

In Sec. II, we develop further our thought experiment of idealized wave growth, to illustrate the uncertainty around the behavior of emitted AN as a type of “wave parameter.” In Sec. III, we introduce the datasets used in the study. In Sec. IV, we present an analysis in which we study the correlation of AN with wave parameters, after the correlation with wind speed is constrained. In Sec. V, we perform multi-linear regression (MLR) to predict AN using combined wind and wave parameters and quantify the benefit of including wave parameters in terms of Root Mean Square (RMS) error in decibels, relative to the baseline performance, using wind information alone. In Sec. VI, we discuss limitations and future work, and in Sec. VII, we summarize our conclusions from this study.

Since Wenz (1962), empirical predictions of AN using wind speed have demonstrated good skill (we use the word “skill” here to indicate the relative ability of different models to predict the observed AN). For example, analysis by Yang (2023), in their Fig. 2(c), finds that 80% of the data are within 3 dB of the mean for 6–8 m/s wind speed at 5 kHz. The present study seeks to provide quantitative information on the potential for inclusion of wave information to improve this skill. We do not present a new predictive formula as was done, for example, by Wenz (1962), Hildebrand (2021), and Pflug (2021). Also, bubbles are not explicitly part of the predictive models developed herein. The problem of noise from breaking waves in the surf zone, e.g., see Deane (1997) and Fabre and Wilson (1997), is not addressed here.

In Sec. I, we introduced a thought experiment with 15 m/s wind speed blowing over two basins, one with 100 m fetch, and the other 100 km, and stated that the latter case will logically result in more bubbles and AN. We elaborate further here. According to prevailing theory, the wave spectrum is controlled by the net effect of three “source terms” that dominate in deep, open water: wind input, wave breaking dissipation, and weakly non-linear four-wave interactions (Cavaleri , 2007). The third is a “direct and inverse cascade,” meaning a transfer of energy from near the spectral peak to the higher and lower frequencies, respectively (e.g., Young and van Vledder, 1993). At higher frequencies, say 1.5 to 3 times the peak frequency, these three source terms roughly balance each other, so that the spectral level is in rough equilibrium with the wind (e.g., Phillips, 1985). At frequencies nearer the peak, in growing seas, the local (in spectral space) imbalance of the three terms produces a spectral downshifting with increasing fetch or duration of a wind event, producing larger wave height (thus, greater wave energy) and longer dominant wave lengths, e.g., Kahma (1981). It might be understood intuitively that breaking in more mature or “fully arisen” seas produces more turbulence, bubbles, and AN than breaking in less energetic windsea, but it is also supported in the literature. Using observational data, Zhao and Toba (2001), in their Fig. 1., show whitecap coverage increasing with wave age. Using fitted curves based on energy balance equations and observational data, Hwang and Sletten (2008) show in their Fig. 3(c), parameter $α$ increases with the duration of a wind event, where $α$ is the coefficient of proportionality between total energy dissipation and the wind speed cubed.

The above suggests that, for any particular wind speed, AN should grow with fetch. However, this is merely a qualitative relation, and so it does not assist us in predicting AN. A numerical ocean wave model such as Simulating WAves Nearshore (SWAN) (Booij , 1999), is based on the energy balance described above and can be used to predict many different wave parameters, each of which is associated with wave breaking to a greater or lesser degree. Table I shows SWAN results for an idealized case with a 10 m wind speed of 15 m/s.

Parameters given in Table I include moments of the wave energy density spectrum $Ef.$ They are computed as $mn=∫f1f2Effndf$⁠, where $f1$ and $f2$ are the upper and lower wave spectra frequency bounds, with 0.038 and 1.0 Hz used here, respectively.¹ Spectral moments typically have a physical relevance. For example, the significant waveheight $Hsig$ is defined as $Hsig=4m0$⁠; $m2$ is proportional to surface orbital velocity, $m3$ is proportional to surface Stokes drift, and $m4$ is proportional to mean square slope of the sea surface. Parameter $Dtot$ is the integrated rate of dissipation of wave energy by whitecapping, and $Tp$ is peak period. The calculation and physical relevance of these parameters are described in greater detail in  Appendix A.

The wave parameters grow at a markedly different rate. For example, $Dtot$⁠, $m4,$ and $m5$ grow very quickly with fetch, while total energy (⁠ $m0$⁠) grows more slowly. For the wave parameter AN generated by breaking waves, the growth rate is unknown, though intuition would suggest a growth rate similar to $Dtot$⁠, $m4$⁠, and $m5$⁠. Of course, it must be recognized that even this unknown idealized growth rate of AN is likely dependent on the acoustic frequency.

This study relies on three datasets: (1) hydrophone observations of AN, (2) observed wind/wave parameters, and (3) modeled wind/wave parameters. Each is described below.

The hydrophones are described in this section. Below, we refer to critical depth (CD), which is the depth at which the sound velocity matches the sound velocity at the surface (Urick, 1984). We briefly review its relevance here. For cases when the ocean depth is below the CD, acoustic energy does not interact significantly with the bottom, and long-distance propagation is supported. In cases where the ocean bottom is above the CD, sound is more significantly attenuated by the bottom. The amount of ambient acoustic energy reaching a hydrophone is maximized when a receiver is shallower than the CD and the ocean bottom is deeper than the CD. When both receiver and ocean bottom are shallower than CD, then most of the acoustic energy that reaches the receiver has interacted with the ocean bottom and is more attenuated. When the receiver and ocean bottom are both below the CD, the receiver therefore lies in the “shadow zone.” Here, most ambient sound waves are refracted away from the receiver, and most received AN is that which is generated locally. An assumption of “direct path only” propagation is appropriate in this case (Wilson, 1983).

The first hydrophone dataset is from a deployment by the NOAA/PMEL at Ocean Station Papa (OSP), south of the Gulf of Alaska. The location is very remote, around 1000 km from the Alaska and Canada mainlands, so fetch here can be limited only by the scale of the meteorological forcing. It is in deep water (4250 m depth), with the hydrophone at 500 m depth, which is above CD. Data were processed for October to December 2018. The processing of these data include eight analysis frequency bands, with center frequencies from 127 Hz to 1.7 kHz, with 5 kHz sample rate.

The second and third hydrophones are offshore of southwest Canada, deployed and managed by Ocean Networks Canada (ONC). Both are near the seafloor, and data were processed for October 2020 to December 2021. One hydrophone is the Cascadia Basin deployment, 202 km west of Vancouver Island. The more nearshore hydrophone is the Clayoquot Slope deployment, 78 km west of Vancouver Island. The wind climate for this region is primarily alongshore, so fetch-limited conditions are not common. These data include 18 analysis frequency bands, with center frequencies from 127 Hz to 18 kHz, with 64 kHz sample rate. At most frequencies, Cascadia is quieter than Clayoquot, implying that at the latter, the impact of being within the sound channel is relatively more important than any suppression of sound by bottom interaction.

Some of the results herein from the OSP and ONC datasets were previously presented in the gray literature (Rogers , 2023).

The remaining four datasets are from hydrophone data provided directly by Dr. John Hildebrand (UCSD, University of California, San Diego). These are the HARP instruments described in Hildebrand (2021).

The instrument locations are shown in Fig. 1 and features of the hydrophone dataset are summarized in Table II. This table includes a simple metric for roughly estimating the fraction of each dataset that is fetch- or duration-limited based on the concept of separating between “forced conditions” (i.e., waves actively receiving energy from the wind) and “unforced conditions.”² This is denoted as “% $Cr/U<1$⁠,” and here, $U$ is the 10 m wind speed, $Cr$ is an ocean wave phase velocity representative of the wave spectrum, calculated using the deep-water relation, $Cr=(gTm,−1,0/2π),$ where $g$ is gravitational acceleration, and $Tm,−1,0$ is a definition of mean wave period, $Tm,−1,0=m−1/m0$ (other definitions of mean period exist; see  Appendix A). The smaller percentage—implying larger $Cr/U$—at the OSP and ONC sites indicates a more mature sea state, implying that a sea state prediction will be less heavily penalized for using only local wind speed.

The procedure for wind noise estimation is similar for all datasets, with modifications to account for factors such as data bandwidth, hydrophone depth, and acoustic propagation characteristics. Raw hydrophone time series data were first processed using Welch's method (Welch, 1967) to produce spectral estimates over stationary time intervals. A Hann window was used to minimize off-frequency-bin energy leakage via spectral sidelobes (Harris, 1978). Received sound levels were then estimated over analysis frequency bands spanning 100 Hz to 20 kHz with center frequencies and bin widths chosen to avoid redundancy. Analysis bin center frequencies for the 18 (eight in the case of OSP, due to lower sampling rate) AN analysis bands are shown in Table III. For all datasets, lower frequencies (<100 Hz) were not processed due to elevated ambient sound levels from shipping activity. High-frequency processing cutoffs were determined based on identified low pass filter roll-off and data acquisition system noise floor. Additional description of the quality control process used on the hydrophone data is given in  Appendix B.

The OSP hydrophone is co-located with a PMEL buoy equipped with an anemometer at 4 m elevation [Ocean Climate Station (OCS) Project Office], and an Applied Physics Laboratory/University of Washington wave buoy (Thomson , 2015), which is Coastal Data Information Program (CDIP) Station 166.³

The ONC hydrophones are not co-located with wind/wave observations, but a number of U.S. and Canadian buoys are in the region, eight of which are used to validate input/output to/from the southwest (SW) Canada wave model described in Sec. III C. The HARP locations similarly have two to five buoys in each area that were used for wave model validation, with the exception of the “Antarc” hydrophone, for which no contemporaneous/nearby in situ wave observations exist.

Wind and wave information from numerical models are generally less accurate than information taken from direct observations. However, in an operational environment, observations are not available at all locations, and one must rely on output from numerical models. This is especially true when a forecast is needed.

For all locations, other than the Gulf of Mexico (GofMX), we use a regional grid of the SWAN wave model (Booij , 1999) nested within a global 1/4° hindcast with the wave model WAVEWATCH III^® (WW3) (Tolman, 1991; ^®). The grid design of Rogers and Linzell (2018) is used for the global grid and resolutions from 1.8 to 2.8 km are used in the SWAN grids. As an example, in the SW Canada (ONC) case, the resolution of the regional model is between 2.5 and 2.8 km and is illustrated in Fig. 2. For the GofMX case, a global model was not needed; instead, a 2.1 × 2.2 km coastal SWAN grid was nested within a 11.9 × 13.9 km regional SWAN grid. The model settings and SWAN grid information are described in the Supplementary Material.

In the case of the two ONC locations, 10 m wind vector forcing fields for the models comes from the Navy's operational global atmospheric model navy global environmental model (NAVGEM) (Hogan , 2014), given at 0.18° resolution. The global wave model used to provide boundary conditions also includes ice fraction information taken from NAVGEM. In the case of the HARP locations, wind and ice forcing is taken from the European Centre for Medium–Range Weather Forecast (ECMWF) Reanalysis v5 (ERA5) (Hersbach , 2020). These fields are given hourly at 1/4° resolution, and the 10-m neutral winds are used.⁴ Model results for the OSP location are not included in this manuscript;⁵ only the observed wave parameters are used for OSP herein.

As mentioned in Sec. III B, with the exception of the Antarc case, the SWAN implementations are validated using wave observations from moored buoys. These buoys are listed in the Supplemental Material.

For further information on the wave model governing equations, source terms, and methods, the reader is referred to Tolman (1991), Tolman and Booij (1998), Booij (1999), WW3DG (2019), and SWAN team (2019).

In the case of OSP observed wind/wave parameters, the hydrophone data are already geographically co-located. For other cases, i.e., those that employ model wind/wave parameters, bilinear interpolation is used for geographic co-location.

Temporal co-locations are made by taking the wind/wave parameters at 30 min intervals, and for each temporal data point, averaging all valid AN data within ±15 min of the wind/wave time. From this, there are between 4000 and 14 000 30 min co-locations for each hydrophone, as listed in Table II. Across all locations, this results in a total of 1291 days of co-locations. Of the co-locations in this total, some include frequencies for which the AN data are marked as invalid by the quality control. The average fraction of valid AN data was 68% for the OSP case, 86%–93% for the two ONC cases, and 95%–98% for the four HARP cases.

Waves that are not fully developed are either fetch- or duration-limited, which may be caused by non-uniformity and non-stationarity (respectively) of the winds. In the classic scenario, fetch is determined by the distance from a nearby coastline. Wind climatology can help us to understand the likelihood that waves are not fully developed, for each of our hydrophone locations, so here, we briefly review this, using information taken from the GLObal Spectral WAve Climate (GLOSWAC) database of Portilla-Yandún (2018). All else being equal, fetch- or duration-constrained (i.e., “young”) wave conditions are more likely to occur when wind speed is high; the “wave age” is typically defined as $C/(Ucosθ)$⁠, where $C$ is wave phase velocity (proportional to wave period if in deep water), $U$ is the wind speed, and $θ$ is the angle between wind and wave directions. Below, for brevity, we use shorthand notation for cardinal directions, e.g., WSW for winds from the west-southwest.

Ocean Station Papa (OSP, NOAA/PMEL): Westerly (especially WSW) winds are most common. Winds are often strong (⁠ $U10$>12 m/s), and occasionally reach 20 m/s.

Cascadia and Clayoquot (ONC): Winds are primarily along the coast, with winds from NW being most common but with the highest winds occurring more often from the southeast (SE).

Gulf of Alaska (GofAK): Winds from east-southeast (ESE) and west (W) are most common, and winds from the north are least common. In the case of westerly winds, there is potential for limiting of wave conditions by fetch from the Alaska coastline, but since the distance is quite large (225–300 km), waves are expected to be fetch-constrained only for the strongest winds (e.g., $U10$>15 m/s).

Offshore of North Carolina (WAT): Similar to the ONC location, winds are often along the coast, with SW winds being most common.

Drake Passage (Antarc): Westerly winds are by far the most common. Winds are strong, occasionally reaching 20 m/s. Clarence Island is only 30 km west of the hydrophone, which strongly implies fetch-limitations, but the island is small, so it would have reduced effect on winds from the northwest (NW) and SW.

GofMX: Winds from the east (E) and SE are most common, though northerly winds are also common. Winds typically do not exceed 16 m/s. The coast is 120–150 km from this hydrophone. Waves are expected to be fetch-constrained only weakly, and only for the highest wind speeds.

Based on the above, only the Antarc location is likely to be fetch-constrained by topography for a significant fraction of the time series used. At other locations, the wave conditions may still be constrained (i.e., not fully developed), but this would occur primarily as a result of non-uniformity and/or non-stationarity of the wind speed and direction.

AN is created by bubbles that are created by breaking waves, and the wave energy is generated by wind action. Thus, the AN is strongly correlated with wind speed, and so models such as Wenz (1962) and Hildebrand (2021), have had success. The direct correlation between AN and wave parameters (e.g., FM95) is difficult to interpret, because much of that AN-to-wave correlation is via mutual correlation with wind speed, which is already used in the AN predictive model. However, it is possible to anticipate the potential for wave parameters to improve the prediction by performing analysis while controlling for this mutual correlation.⁶ We have done this using methods described in Secs. IV A and IV B.

Figure 3 shows examples from the OSP observations at a middle frequency, 841 Hz. The color scale here is the wave parameter $m3$⁠, normalized by the mean $m3$ within a wind speed bin of width 1 m/s. Horizontal black lines indicate these bins, and the mean AN in the bin, and thus, these black lines can be taken as an example of an AN prediction that is based on wind speed alone. Below we discuss a “wind speed only” model, which refers to a hypothetical step-function empirical model that is composed of these black lines. Figure 4 shows results for the ONC model output at two frequencies, 841 and 4700 Hz. Since this case has the larger population (10 900 vs 4400), plotting individual points would result in many concealed points, so instead, the normalized $m3$ values are plotted as the averages on a 1 m/s × 2 dB grid.

The $m3$ parameter was selected for this example based on analysis that we will subsequently present in Sec. IV B and Fig. 6. It is not the only wave parameter that shows the trend shown in Fig. 3. The parameters $m2$⁠, $m4$⁠, $m5$⁠, and $Dtot$–and to a lesser extent, $m0$ and $m1$—also show this. Similar example plots, for $m2$⁠, $m4$⁠, $m5$⁠, and $Dtot$ are provided in the Supplementary Material. Examples for a higher frequency, 7.2 kHz, are also provided in the supplementary material.

Figure 3 shows that at higher wind speeds, there is little variation of AN from the black horizontal lines. This indicates that a wind speed only prediction of AN is already quite accurate, and little can be done to improve upon it. At lower wind speeds, there is much more scatter above and below the black lines. The scatter is presumably due to sources other than wind. The colors indicate that where observed AN is higher than the wind speed only model, normalized $m3$ tends to be higher (red and orange), and below the line, the opposite is true (blue colors). This suggests that some improvement can be made on the wind speed only model by incorporating $m3$ into the prediction model.

Figures 3 and 4 indicate a dependence of AN on wind/wave parameter even at low wind speeds. This may seem surprising given that whitecaps typically start to form only at wind speeds of 4–6 m/s (e.g., Monahan, 1971); 5 m/s is sometimes used as a rule of thumb. A study by Updegraff and Anderson (1991) sheds light on this, and here, we rely on a summary given in Callaghan (2018), in the context of so-called microscale breaking at low wind speeds (emphasis added): “There is acoustic evidence that small oceanic breaking waves occurring at wind speeds as low as 1.5 m/s may entrain small amounts of air but do not form distinctive whitecaps (Updegraff and Anderson, 1991).”

A limitation of the diagrams of Figs. 3 and 4 is that these only present results for a single wave parameter and single acoustic frequency at a time. Here, we introduce a method that permits averaging over frequencies. First, we group our dataset into 1 m/s bins, essentially looking at vertical slices of plots such as Fig. 4, and then if the population of a bin is not small—we use a threshold of 200 here—we compute the Pearson correlation between AN and the wind/wave parameter. Since correlation with $U10$ is constrained via the binning, we refer to this computed correlation as “residual correlation.” Examples for $U10$ and $m3$ are shown in Fig. 5(a), for the wind speed bin of 4–5 m/s. Importantly, there is still some correlation between AN and wind speed within this 1 m/s bin: $RP$ = 0.10. This establishes a baseline against which to judge the correlation of AN with the wave parameter. We find a correlation for normalized $m3$ of $RP$ = 0.46. The substantial exceedance of the wave parameter's $RP$ over the wind's $RP$ can be taken as an indication that there is correlation of AN with the wave parameter independent of their mutual correlation with wind speed. However, this result is not seen across all wind speeds. In Fig. 5(b), we plot the correlation values as a function of wind speed and, consistent with what we saw in Figs. 3 and 4, the $RP$ of the wave parameter exceeds that of the wind speed only at lower wind speeds (⁠ $U10$ from 2 to 7 m/s).

The correlation functions $RP(U10)$ at individual frequencies, such as Fig. 5(b), can then be averaged to represent all frequencies in a single plot. Figure 6 shows the mean magnitude of $RP(U10)$ for several wind/wave parameters. The result for wind speed is included as a thick blue line, to indicate the range in which wave parameter correlation can be considered insignificant (equal to or smaller than thick blue line). This comparison includes two parameters that have not been previously introduced, summarized here: “STEEP” is the mean steepness $Hsig/L$⁠, where $L$ is mean wavelength, and “precip” is mean precipitation rate. The parameters are described in detail in  Appendix A. A comparison similar to Fig. 6, for five additional wave parameters, is included with the supplementary material.

This figure indicates that parameters such as $m2$⁠, $m3$⁠, $m4$⁠, mean steepness, and integrated dissipation (⁠ $Dtot)$ can help to improve the AN prediction, by a modest amount (correlation from 0.25 to 0.42), at lower wind speeds (2–7 m/s). We find that despite the effort to remove precipitation during the acoustic data processing stage (see  Appendix B), there is significant correlation with precipitation, even at higher wind speeds.

In the moderate to high wind speeds, 10–15 m/s, $Dtot$ is the only wave parameter in Fig. 6 that has correlation $RP,$ elevated relative to that of wind speed, and the difference is slight. Taking example values that are used in the calculations of the average of Fig. 6; for the 10–11 m/s bin, at 841 Hz, the $RP$ for $Dtot$ and $U10$ are 0.23 and 0.14, respectively. At 1.7 kHz, it is 0.22 and 0.14, respectively. This suggests that $Dtot$ has only slight potential for improving AN prediction in this wind speed range.

For each location, we perform linear regression to fit individual wind/wave parameters to AN and MLR to fit multiple parameters (e.g., wind speed plus one to five wave parameters) to AN. All fits were to sound levels (AN) as measured in decibels. For each regression, using the notation that $C1$ is the first in a series of fitting coefficients, and $p1$ is the first in a series of $n$ input parameters (e.g., wind speed), we create both a linear fit, $ANlin=∑i=1nCnpn=C1p1$ +  $C2p2$ +⋯, and a log fit, $ANlog=∑i=1nCnlog10pn=C1log10p1$+ $C2log10p2$+⋯, and select the one with higher skill. We also experimented with power fits, but those results are not presented here. The parameter sets included here have $n$ as low as $n=1$ (the case with $p1=U10)$⁠, and as high as $n=8$ (the case with $p=U10,m−1, m0, …,m5$⁠). For example, in Sec. V A, we have seven parameters sets, and with 18 frequency bands, this means that for every hydrophone location, 126 fits are recorded, each being of type $ANlin$ or $ANlog$⁠.

MLR was performed for all seven locations, but the OSP results are excluded here. This is for two reasons. First, the frequency range (127 Hz to 1.7 kHz) is smaller than that from the other cases (127 Hz to 18 kHz), and second, AN levels at the lowest frequencies at the OSP location are elevated relative to levels anticipated by models such as Wenz (1962) and Pflug (2021), perhaps from ubiquitous, distant shipping. Both factors complicate the presentation of MLR results from OSP alongside the other cases, in particular the “cross-application” section below. Because OSP is excluded, all MLR analysis herein is performed using modeled wind and wave parameters.

Below, we present results of MLR using two broad approaches. In the first approach, that of “local application,” we calibrate using a dataset from a location and then apply the MLR model to independent data from the same location. Each AN time series is segmented by “day of month,” with a separator at day of month equal to 21.0, so that roughly 2/3 of data are used for creating the regression, with 1/3 withheld and subsequently used for prediction and evaluation.

In the second approach, that of “cross-application,” we calibrate MLR models for each location individually and apply the models at all six locations, including the location used for the calibration, resulting in 36 calibration/application pairs for each model type. No data are withheld.

Hildebrand (2021) use an explicit formula for sound absorption, and Wenz (1962) divides data into two water depth categories: deep and shallow. Since our prediction models do not control for sensor depth, we should anticipate that they suffer a penalty to overall accuracy. However, since the penalty applies to both the baseline model (using $U10$ alone) and the experimental models (using $U10$ plus wave—and, in some cases, precipitation—parameters), general conclusions about relative skill are still valid. Unlike those earlier works, our goal here is not to create a predictive model that can be applied to other locations. In both cases (local application and cross-application), the goal is the same: to quantify the impact of inclusion of wave parameters on model predictive skill.

Figure 7 shows evaluations of the multi-linear predictions, comparing the wind speed only model against six examples of experimental models. The vertical axis is the RMS error of the AN prediction, in dB re: μPa²/Hz. The plot includes two wave parameters that have not yet been introduced: “FSPR” is a metric for frequency narrowness of the wave spectrum and “DSPR” is a metric for directional spread of the wave spectrum. The parameters are described in detail in  Appendix A.

The wave parameter set of spectral moments $m−1, m0, …,m5$ was chosen as a list that we expected to adequately describe spectral shape. The parameters sets with $U10,$ plus three other parameters, were selected based on early experiments with the MLR, in which the best performing first parameter (which was $m5$⁠) provided a starting point, and then the “most complementary” second parameter was chosen (⁠ $m−2$⁠), and then several of the most complementary third parameters (FSPR, DSPR, STEEP, $Dtot$⁠, precip) were added. Here, most complementary refers to root mean square error (RMSE) reduction for that preliminary dataset, and we expect that the outcome was determined, at least in part, by the orthogonality of the parameters.

We find that there is modest but consistent improvement in the predictive skill, relative to the baseline model. The degree of improvement is mostly insensitive to the selection of wave parameters.

Some conclusions and conjecture can be made from the results shown in Fig. 7:

• The improvements vs the baseline model are typically less than 0.3 dB RMSE.

• Improvements are typically small at the lowest and highest frequencies, and largest at the middle frequencies. This implies that our dataset is well-bounded for our problem, i.e., extending the experiment to lower or higher frequencies is not justified.

• The start and end of this mid-frequency range varies from one location to another.

• The selection of additional parameters has only minor impact, i.e., it appears that one set of wave parameters is as good as another.

• The Antarc case is somewhat of an outlier: the baseline error is relatively high (over 4 dB at 2.5–6 kHz), as is the degree of improvement (more than 0.7 dB at 2–4 kHz). This result is not surprising given that the region of the hydrophone has sea ice present during approximately 50% of the time period in question. While the baseline model does not account for the impact of sea ice on wind noise, the wave model does.

Though not shown here, we also compared the accuracy of models with just two input parameters, e.g., $U10$ and $Dtot$⁠. As a general rule, these were slightly less accurate than the models shown in Fig. 7, which use four or more input parameters.

With the cross-application experiments, there are 36 combinations per model, making it impractical to present in the format of Fig. 7. Instead, we use tabulations of RMSE and RMSE averages. Table IV presents the RMSE for each combination for a single frequency bin near the middle of the frequencies used here, 2.85 kHz, and showing a single experimental model, the one using the parameter set $U10,Dtot.$⁷ The table indicates that the experimental model is an improvement over the baseline model using parameter set $U10$ for 29 of the 36 combinations. In the case of the baseline model, the least successful calibration is that which was made using the Antarc case, with a mean RMSE of 3.71 dB. This is unsurprising, since, as mentioned already, the baseline model does not take into account the effect of sea ice on wind noise. The experimental model $U10,Dtot$ calibrated to Antarc outperforms the baseline model $U10$ calibrated to Antarc for all application cases (mean RMSE of 3.29 dB, reduced by 11.3% from the RMSE of the baseline). Of the experimental model calibrations, the calibration to the WAT case is most successful (mean RMSE of 3.14 dB). For cases other than Antarc, improvement is marginal: 0.15–0.19 dB from a baseline of 3.32–3.71, or a 4%–5% improvement from the baseline.

Table V shows results for all frequencies, for the same two models (⁠ $U10$ baseline and $U10,Dtot$ experimental). We see that the experimental model is more accurate than the baseline for 16 out of 18 frequencies, and for all frequencies $≥$ 408 Hz. To show results from several experimental models, we must further aggregate the outcomes. This is done in Table VI, by indicating the number of frequency bins for which the experimental model outperforms the baseline model, using the mean across all 36 calibration-application combinations. Here, we find that the model with only two inputs (⁠ $U10,Dtot)$ outperforms the same model with a third input parameter (⁠ $m3, m4,m5,$ or DSPR), which differs from our conclusions from the local application study. This may imply that the models with several input parameters are “over-tuned” to their specific locations.

Non-local AN generation is not addressed here. All past studies have, to our knowledge, also omitted this. In principle, this can be addressed with an acoustic propagation model, but the calibration process would become prohibitively complex.

The degree of reduction of RMSE in AN predictions by including wave parameters is found to be marginal herein, but may be larger for fetch-constrained cases. With the possible exception of the Antarc case, datasets that we have used here are not severely fetch-constrained by proximity to land. We recommend that additional experiments be performed for cases where the hydrophone is less than 100 km from shore, and off-shore winds are climatologically predominant.

We find here that the degree of improvement to model predictive skill by including wave parameters is increased in specific cases: when predicting AN for frequencies mostly falling in the middle of our frequency range, and for regions that are sometimes ice-covered. Further, our correlation study suggests that more improvement can be anticipated for cases of lower winds. AN prediction tools using wave information that targets these specific cases, may be most fruitful. To this end, it would be worthwhile to experiment with frequency-space fittings at discrete wind speed bands, in lieu of the MLR fittings at discrete frequency bands used herein.

Our MLR models treat the wave parameters as the independent variables, along with $U10$⁠. It may be useful to experiment with instead using the wind-normalized wave parameters, as was done in our correlation study (e.g., Fig. 3).

Where both frequency and wind speed are high, it has been shown that there is a relative reduction in the sound levels. Thus, AN is non-monotonic in wind speed dependence at high frequencies. The earliest reference to this phenomenon was, to our knowledge, from Farmer and Lemon (1984). They propose that the effect, observed at 14.5 and 25 kHz, is caused by the scattering and absorption of sound by near-surface bubbles from wave breaking, and we find no studies since then that contradict this explanation. In the parametric AN model of Hildebrand (2021), this starts to appear at frequencies above 10 kHz and $U10$>15 m/s. Yang (2023) find a “fast drop-off of [ambient noise level] for frequency higher than 4 kHz when wind speed exceeds 15 m/s.” Inspecting our ONC Clayoquot and HARP GofMX datasets, we find that AN is unambiguously non-monotonic with wind speed at the frequency bands 14 and 18 kHz, for wind speed greater than 15 m/s; these plots are included in the Supplemental Material.

The non-monotonic behavior is a concern in the context of our MLR fitting, which, like Wenz (1962), is monotonic with wind speed. The natural outcome is reduced accuracy of the fit for these data points (high-frequency, high wind speed). Hildebrand (2021) use more advanced fitting methods, with fitting parameters that are a function of wind speed and frequency, but this is not practical for a regression to multiple wave parameters so it was never considered for this study. Machine learning (ML) methods may be able to address this problem.

We make the following conclusions, which are based on comparisons that are omitted from this manuscript for sake of brevity:

• In our evaluations of RMSE for the OSP case (analogous to Fig. 7), we find significantly higher skill (up to around 0.5 dB) using input from observations (wind only or wind and waves) relative to using model values of the same. This indicates that accuracy (or lack thereof) of the winds is a primary determiner of AN predictive skill for all cases. This comparison is included in the Supplementary Material.

• We found that prediction of AN from wind/wave parameters using ML exhibits skill that is not significantly better than that using the simpler method of MLR. However, since there is great variety in ML methods, this is an approach that should be explored further using our dataset.

Herein, wave age is used for general description of climatology (Table II), and our thought experiment (Sec. II) centers on the issue of wave parameter development relative to a fully developed condition. However, wave age is not included as an input parameter to our MLR models, which may seem counter-intuitive. There are two reasons for this. First, the actual wave age depends on both wave frequency and direction, i.e., it is a “wave age spectrum,” $C(f)/Ucos(θ)$⁠, where $C$ is wave phase velocity and $θ$ is the angle between the wind and the wave component. One can imagine a careful calculation as an energy-weighted integration of this spectrum; we have not done this. Second, the dependence of AN on wave age should be non-monotonic, making any linear fitting problematic. This was first shown by Dragan-Górska (2024) who studied the trend for two datasets: one shows no trend (⁠ $R=−0.04$⁠) and another shows negative trend (⁠ $R=−0.65$⁠). Our assertion for non-monotonic dependence is based on the following. At very small fetch (or very small wave age), the AN generation will also be very small, as the waves are mere ripples. On the other hand, when wave age is large, the waves are swell. These waves typically do not break in the open ocean, implying that they do not generate significant AN. As AN is small at the two extremes, monotonicity is impossible.

The present study is concerned with the prediction of so-called wind noise, though other sources contribute (Sec. I). These sources are not of primary interest here, but there are indirect connections. For example, our data processing ( Appendix B) included removal of data contaminated by sources such as precipitation and shipping. Shipping noise is also connected to sea state insofar as most vessels try to avoid heavy seas. In the case of sea ice, we have noted how ice cover affects the AN generated by wave breaking. Ice also generates AN: thermal-induced cracking creates AN across a broad range of frequencies (e.g., 40–600 Hz or even broader; see Milne, 1972), while cracking, shearing, and other mechanical effects produce AN at lower frequencies (under 200 Hz; Greening and Zakarauskas, 1992; Sheng , 2023). The mechanical effects can be produced by pressure ridges (from horizontal pressure in the ice pack) and motions induced by nonstationary atmospheric pressure. Ridging, in turn, is produced by effects such as ocean currents, wind drag, and wave forcing. In the latter case, as waves enter the marginal ice zone, they are attenuated, and the lost wave momentum is transferred either to the ocean circulation or to the ice. This “push” by the waves on the ice is discussed in Stopa (2018). For further information more generally on AN from sea ice, the reader is referred to Roth (2012), Kinda (2013), Geyer (2016), Ozanich (2017), Tian (2019), Wen (2020), Pajala (2021), and Tollefsen and Buen (2022).

We conclude the following, based on results presented:

• We used correlation analysis (constraining the correlation of AN with wind speed, to identify correlation with wave parameters that is independent of mutual correlation with wind speed) and one method of AN prediction (MLR models). Both methods point to the same conclusion: that prediction of AN can be improved by including wave parameters, but that the degree of improvement is marginal.

• A corollary to above: our experiments robustly indicate that information about AN exists in the ocean wave state that is not contained in the wind speed alone.

• Likely, the process that leads to AN (wave breaking and bubbles) reaches a mature state more quickly than, for example, wave height. This implies that there is less penalty for using “wind only” prediction of AN than for using “wind only” prediction of wave height.

• The correlation analysis indicates that the potential for improving AN prediction using wave parameters is greatest for lower winds (e.g., 2.5–7.0 m/s).

• MLR models are calibrated to specific locations. A consistent improvement via introduction of wave parameters is seen both when applying to the same location (local application) and applying to all locations (cross-application). However, the improvement is marginal. For the local application, typically RMSE is reduced by less than 0.3 dB and/or less than 12% of the original RMSE.

• In the local application of MLR, improvements via addition of wave parameters are typically small at the lowest and highest frequencies, and largest at the middle frequencies (e.g., 0.6–6 kHz for the Antarctic location).

• In the local application of MLR, RMSE of the baseline model (using $U10$ only) is especially high for the Antarctic case. We believe that this is because the baseline model does not account for the impact of sea ice on wind noise. Much like fetch, sea ice is a feature of the ocean that can potentially have strong impact on waves and wave breaking, and the improvement of AN prediction via introduction of wave information is therefore relatively high at this location.

• In the cross-application of MLR models, the models using only two or three input parameters are beneficial over a larger frequency interval than one that relies on six input parameters, which suggests a risk of over-tuning.

See the supplementary material for additional details about hydrophone moorings and model design; band-normalized plots similar to Fig. 4, but for other wave parameters and one other acoustic frequency; correlation results comparable to Fig. 6, for five additional wave parameters; plots demonstrating the non-monotonic dependence (roll-over) of AN with wind speed; and a comparison of MLR results using observed input parameters instead of modeled input parameters.

We are grateful that all observational data used in this study are provided freely to the research community, and acknowledge the significant effort and resources required to make this available. We thank NOAA PMEL (R. Dziak and others), NOAA OST (J. Gedamke), and Ocean Networks Canada for their hydrophone datasets; APL/UW (Jim Thomson and others) for the OSP wave buoy dataset; and the Ocean Climate Stations (OCS) Project Office of NOAA/PMEL (Meghan Cronin and others) for the OSP meteorological dataset. The authors gratefully acknowledge Dr. John A. Hildebrand, UCSD, who was generous with both his data and his time, providing copies of the noise recordings from HARP data loggers. The HARP datasets were each supported by a different group: GOMRI (Gulf Of Mexico Research Initiative), Pacific Fleet, South East Fisheries Center, and collaboration with Chile. David Wang (NRL 7330) contributed to prior discussions of energy moments. Aubrey Geary (Naval Research Laboratory summer intern in 2023) performed machine learning experiments. Gregg Jacobs (NRL 7320) provided useful suggestions for visualization of results. Careful reading by anonymous reviewers resulted in significant improvements to the manuscript. This work was funded by the Office of Naval Research through the NRL Core Program, Program Element 0602435N, Work Unit 6C32, Task Area UW-435-025. The NRL Core 6.2 project is “Ambient noise prediction using wave models.” This paper is contribution NRL/7320/JA–2024/7 and has been approved for public release.

The authors have no conflicts to disclose.

The data that support the findings of this study are available in Mendeley Data at doi:10.17632/4v72pgswc6.1. In the case of the HARP datasets, NRL is able to provide the NRL-processed AN data if the party making the request has explicit permission from the provider of the original AN dataset, Dr. John Hildebrand, UCSD.

Here, we provide definitions of the wave parameters used in this manuscript, and where practical, also describe their physical significance. In Table VII, the parameters that are derived from moments of the non-directional ocean wave spectrum are given.

Other parameters are as follows:

Peak period is defined as $Tp=1/fp$⁠. The peak frequency $fp$ is the frequency $f$ at which spectral density $Ef$ is maximum.

Integrated dissipation is $Dtot=∫Sdsfdf$⁠. This is proportional to the dissipation of wave energy by deep-water breaking (whitecapping). Here, $Sdsf$ is the wave model source term associated with this dissipation. Units of $Sdsf$ and $Dtot$ are m² and m²Hz, respectively.

Here, $Ef,θ$ is the directional form of this spectral density, given in units of m²s/rad or m²s/degree, $Ef,θ=E(f) Df,θ$⁠, and $∫Df,θdθ≡1$⁠.

Frequency narrowness/spread parameter with shorthand FSPR is FSPR =  $∫Eωeiωτdω/m0$⁠, where $ω=2π/f$ and $τ=Tm,0,2$⁠. This is taken from Battjes and van Vledder (1984). Contrary to the name, a narrower spectrum implies a higher FSPR (e.g., Rogers and van Vledder, 2013), but we retain the notation FSPR for consistency with notation of SWAN team (2019). A narrower spectrum is generally associated with a more amplitude-modulated or “groupy” wave trains, e.g., Gemmrich and Thomson (2017).

Steepness with shorthand STEEP following the SWAN team (2019), is computed as $Hm0/λm$⁠, where $λm$ is the mean wavelength $λm=2π(m0/∫kEfdf)$⁠, and $k$ is the wavenumber at frequency $f$⁠, computed from the linear dispersion equation.

Precipitation, with shorthand, precip, here, is taken from ERA5 reanalysis. This is the mean rate of precipitation. Though the units have no bearing on results presented herein, we apply the MLR using “precip” measured in millimeters per hour.

Calibrated spectral estimates were formed from the raw spectra based on Transmission Voltage Response (TVR) curves of the receivers. The calibrated spectral intensities were linearly averaged over 1 min intervals and band levels were derived as the linear mean of the time averaged spectra for the analysis frequency bands (eight for OSP and 18 for ONC and HARP). For the PMEL data, this averaging was done using 1 s FFTs (fast Fourier transforms) with 50% overlap (119 snapshots per minute, 1 Hz frequency resolution) while the ONC and HARP data were processed using 0.1 s FFTs with no overlap (600 snapshots per minute, 10 Hz frequency resolution). Finally, an algorithm was developed for automated removal of wind noise contaminants based on statistical techniques for outlier removal. Outlier selection is based on signal normalization techniques and median-based statistical measures; the wind speed is not used directly. These measures include a sliding split window normalizer (Bentrem , 2018) combined with the median absolute deviation (MAD) statistic (Rousseeuw and Croux 1993) and mean-to-median ratio (MMR) statistic (Schwock and Abadi, 2021a,b), utilized for threshold-based removal of time-frequency localized interference. For all site data, time domain transient detection was applied to each bin of normalized band averaged results. The MAD statistic was independently computed, in decibels, across all relevant frequency bins of the averaged spectra, and the MMR statistic was formed as the ratio (decibel difference) of the averaged power spectra and correspondingly computed median spectra. Thresholds applied for each detection criteria were empirically chosen to maximize interference removal while minimizing collateral damage due to the erroneous exclusion of valid wind noise data. Figure 8 shows example spectra from the Clayoquot Slope site along with excluded points (dark blue and magenta bins/regions). Indicated regions of precipitation and shipping noise, along with other sources of interference, were largely removed across all datasets by automated pruning. Figure 9 displays corresponding band wind noise estimates and coincident wind speed (magenta curves).