Atmospheric Sound Propagation Over Rough Sea: Numerical Evaluation of Equivalent Acoustic Impedance of Varying Sea States

1. Introduction

For many decades, atmospheric sound propagation studies have been conducted predominantly on land [1,2,3,4], with only a few studies conducted over sea and small bodies of water [5]. In recent years, however, research in sound propagation in coastal areas has garnered increased interest, driven by the growth of off-shore wind farms [6,7,8], oil and gas platforms, offshore airport construction [9], and other activities. The study of sound propagation over sea is also critical for detection of vessels at sea in cases of low visibility, as well as for echolocation over long ranges, on the order of kilometers.

Atmospheric sound propagation over water is strictly coupled to meteorological conditions, which are responsible for acoustic refraction [10,11]. Another important influence is the geometry of the water surface, which is responsible for acoustic scattering when acoustic wavelengths are comparable to the characteristic size of the water roughness. Parabolic Equation (PE) methods, such as the CNPE (Crank Nicholson PE) [12] and the GFPE (Green’s Function PE) [13], have been commonly used to model acoustic refraction. However, those methods are applicable to ranges with smooth variations of the ground level, and thus not directly able to account for arbitrary geometries. The Beilis–Tappert PE method [14] can efficiently model atmospheric sound propagation over irregular terrain in cases with local slopes less ${20}^{\circ }$ [15]. However, this work considers higher sea states with slopes exceeding ${20}^{\circ }$. To account for local slopes beyond ${20}^{\circ }$ but not exceeding ${30}^{\circ }$, the CNPE model can be modified to obtain the GTPE (General Terrain PE) [16]. The GTPE modification however comes with significant computational cost. The GTPE method can be used to quantify the effect of surface acoustic scattering during atmospheric sound propagation. In the case of propagation over the ocean, the roughness can be considered by generating pseudorandom surfaces consistent with the sea state. The combination of GTPE with random surface roughness has been applied to underwater acoustics [17], but the application to the atmospheric problem has not been reported. The large computational cost of this method makes it unsuitable for real time high frequency (above 1 kHz), long range (on the order of kilometers) atmospheric sound propagation. To reduce computational cost, a flat surface with an equivalent impedance can be used [6,8]. Validation of this approach applied to water roughness is lacking in literature.

In this paper, a numerical method is proposed to estimate equivalent surface impedances for water roughness representing sea states from 2 to 5. One previous approach[6] found an equivalent surface impedance using the Delany and Bazley [18] model. The impedance was used to inform PE predictions compared to long range acoustic measurements. A second approach by Van Renterghem et al. [8] used Finite Difference Time Domain (FDTD) simulations over short distances (source and receiver separated by 5 m) to replicate a method previously reported by Boulanger and Attenborough[19] to find equivalent impedances of pseudorandom sea surfaces. The equivalent impedance determined with the FDTD method was implemented in PE predictions of atmospheric acoustic propagation above water over a 10 km range. The method proposed in this work differs from Van Renterghem et al. in that this work uses GTPE simulations over 1 kilometer range to estimate the equivalent impedance of different sea states. This work overcome limitations noted in previous methods.

The paper is organized as follows. Section 2 details the approach used to model sea surfaces. Section 3 reviews the Monin-Obukhov similarity theory and discusses the meteorological data used in this work to inform a PE solver. Section 4 describes how acoustic excess attenuation is predicted over rough sea using a GTPE solver. Section 5 reviews methods used to estimate equivalent impedances for rough seas and presents this work’s approach. Section 6 and Section 7 present the results and conclusions of this work.

2. Modeling Sea Surfaces

Existing models that relate wave height to wind speed have been studied for open waters, like seas and oceans[20,21,22]. In this paper such models are used to predict water roughness. The assumptions made throughout this paper are: (1) the sea is fully developed, (2) predominant wind waves are aligned with the mean wind direction, and (3) water depth effects are negligible.

2.1. Background on Sea Spectra

PM Spectrum

In the absence of breaking waves, the sea surface can be described by a sea spectrum, a frequency distribution of energy. The most popular sea spectrum for fully developed seas is the Pierson and Moskowitz (PM) spectrum [20], $S_{PM}\left(k\right)$, defined in Equation 1

$$S_{PM}\left(k\right)=\frac{\alpha }{2k^3}exp[-\beta (\frac{g}{k}{)}^2\frac{1}{U_{19}^4}]$$

(1)

where k is the angular spatial frequency, $\alpha =8.1\times {10}^{-3}$, $\beta =0.74$, g is the acceleration of gravity, in m/$s^2$, and $U_{19}$ is the wind speed in m/s measured at 19 meters height. The spectrum depends only on $U_{19}$, which can often be replaced by wind speed measured at 10 meter height with the approximation of $U_{19}\approx 1.026U_{10}$. This model has been applied to sound scattering predictions due to water roughness in both atmospheric and underwater problems[8,17,23]. Despite the broad use and a recent revisitation by Alves et al., the PM spectrum is unsuitable for describing high frequency gravity waves or capillary waves [24]. Water roughness in the high-frequency range can have a non-negligible influence on atmospheric sound propagation.

ECKV Spectrum

A model that is better suited to describe the high frequency sea waves is given by the Elfouhaily, Chapron, Katsaros, and Vandemark (ECKV) [22] sea spectrum. It is also important to note that for a fully developed sea, the ECKV spectrum, like the PM spectrum, is only a function of wind speed. The ECKV spectrum, $S_{ECKV}\left(k\right)$, shown in Equation 2, accounts for two frequency regimes.

$$S_{ECKV}\left(k\right)=\frac{B_l+B_h}{k^3}$$

(2)

where $B_l$ is the low frequency portion of the spectrum that considers the gravity waves, and $B_h$ is the high frequency portion of the spectrum that considers capillary waves. These two portions of the spectrum are similarly modeled as:

$$\begin{array}{ccc}\hfill B_l& =& \hfill \frac{1}{2}{\alpha }_p\frac{c_p}{c\left(k\right)}{F}_p\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill B_h& =& \hfill \frac{1}{2}{\alpha }_m\frac{c_m}{c\left(k\right)}{F}_m\hfill \end{array}$$

(4)

In Equation 3, $B_l$ depends on the generalized Phillips-Kitaigorodskii equilibrium range parameter for long waves[25], ${\alpha }_p$, which in turn depends on the dimensionless wave-age parameter, ${\Omega }_c$, such that ${\alpha }_p=0.006\sqrt{{\Omega }_c}$. In the case of a fully developed sea, ${\Omega }_c=0.84$. The wave phase speed, $c\left(k\right)$, is given by $\sqrt{(g/k)(g+{(k/k_m)}^2)}$. The wavenumber, $k_m$, which corresponds to the gravity-capillary peak, is a constant 370 rad/m. The phase speed of the dominant long waves $c_p$ can be found as $c\left(k_p\right)$, where $k_p=g{\Omega }_c^2/U_{10}^2$ is the wavenumber of the spectral peak. The long-wave side effect function, $F_p$, is defined as:

$$F_p=L_{PM}{J}_{p}exp\left\{-\frac{{\Omega }_c}{\sqrt{10}}\left[\sqrt{\frac{k}{k_p}}-1\right]\right\}$$

(5)

where $L_{PM}$ is the Pierson-Moskowitz shape function and $J_p$ is the Joint North Sea Wave Project (JONSWAP) [21] shape function. These functions are defined as:

$$\begin{array}{lll}{L}_{PM}\hfill & =& exp{\left[-\frac{5}{4}\left(\frac{k_p}{k}\right)\right]}^2\hfill \\ J_p\hfill & =& {\gamma }^{\Gamma }\hfill \end{array}$$

(6)

where $\gamma$ is $1.7$ for a fully developed sea and $\Gamma$ is given by Equation 7.

$$\Gamma =exp\left\{-\frac{{\left(\sqrt{\frac{k}{k_p}-1}\right)}^2}{0.0128{(1+4{\Omega }_c^{-3})}^2}\right\}$$

(7)

In Equation , $B_h$ depends on the generalized Phillips-Kitaigorodskii equilibrium range parameter for short waves, ${\alpha }_m$ given in Equation 8

$${\alpha }_m={10}^{-2}\left\{\begin{array}{cc}\hfill 1+ln(u^{*}/c_m)\hfill & \hfill \mathrm{for}{u}^{*}<c_m\hfill \\ \hfill 1+3ln(u^{*}/c_m)\hfill & \hfill \mathrm{for}{u}^{*}>c_m\hfill \end{array}\right.$$

(8)

where $c_m$ is the minimum short wave phase speed associated with the gravity-capillary peak, and it is given by $c\left(k_m\right)$ = 0.23 m/s. Friction velocity, $u^{*}$, is given by Equation 9

$$u^{*}=\sqrt{C_d}(U_{10}-U_0)\approx \sqrt{C_d}{U}_{10}$$

(9)

where the wind speed at ground level, $U_0$, is usually negligible. As this work is concerned with sound propagation over water, the drag coefficient, $C_d$, is therefore representative of wave-wind interactions. For wind speed $U_{10}$ greater than 1 m/s, Equation 10 has been widely used to estimate drag coefficient over water[26].

$$C_d=(0.8+0.065U_{10}){10}^{-3}$$

(10)

Finally the short-wave side effect function, $F_m$, stated in Equation 4 is given by Equation 11.

$$F_m=L_{PM}{J}_{p}exp\left\{-\frac{1}{4}{\left[\frac{k}{k_m}-1\right]}^2\right\}$$

(11)

Spectra Comparison

Differences between the PM spectrum and the ECKV spectrum for fully developed seas are presented in Figure 1. The largest discrepancies between the spectra can be observed at sea spatial frequencies larger than approximately 10 rad/m. Acoustic scattering becomes relevant when the acoustic wavelength is comparable to the physical wavelength of the sea surface. For a sea spatial frequency of 10 rad/m, the wavelength of the sea matches the wavelength of sound in air at 550 Hz. This study considers acoustic frequencies in the range 125-2000 Hz. The ECKV spectrum is therefore chosen, as it was developed to model high frequency sea waves, with wavelengths comparable to that of atmospheric sound in the frequency range considered in this work.

2.2. Generating Pseudorandom Sea Surfaces

Pseudorandom sea surfaces can be generated using the inverse spatial Fourier transform (FFT) of power sea spectra as described in previous works [8,27] This approach requires a uniform distribution of spatial frequency, k. The resolution, $\Delta k$, and the minimum spatial frequency, $k_{min}$, of this distribution are both equal to $2\pi /L$, where L is the horizontal length of the generated surfaces. As a consequence, the spectra are always unevenly sampled, because of the logarithmic nature of the wave energy distribution. In other words, using a uniform distribution of k leads to an undersampling of the low frequency side of the spectrum, which contains the majority of the energy.

In this work, an alternative approach to generate pseudorandom sea surfaces is proposed, similar to that of Kay et al. [28]. The method consists of replacing the wave height, $z\left(r\right)$, with a discrete sum of N sinusoids, as shown in Equation 12.

$$\begin{array}{ccc}\hfill z\left(r\right)& =& \hfill \sum _{n=1}^N\Re \left[\sigma A_n{e}^{i(k_{n}r+2\pi \rho )}\right]\hfill \\ \hfill \mathrm{where}{A}_n& =& \sqrt{S\left(k_n\right)\Delta k_n}\hfill \end{array}$$

(12)

where the symbol ℜ stands for real part. The sinusoid amplitudes, $A_n$, are given by the spectral density contained in a finite frequency interval. The argument depends on a spatial coordinate r. In Equation 12, $S\left(k_n\right)$ is the average value of a wave spectrum $S\left(k\right)$ over the frequency band centered at $k_n$ with width $\Delta k_n$. The amplitudes are multiplied by $\sigma$, a random number drawn from a normal distribution. This distribution has unit mean and standard deviation of 0.2. The standard deviation of 0.2 satisfies the condition that the maximum wave height, $H_{max}$, is approximately 1.9 times the significant wave height, $H_s$, over 1000 waves[29]. Finally, the phase factor is obtained from $\rho$, a random number drawn from a uniform distribution with range [0 1].

The method presented in this work does not require the resolution $\Delta k$ to depend on the spatial coordinate r. Therefore, it allows generation of pseudorandom sea surfaces for any range L without affecting the spatial frequency resolution. With this method, both linear and non-linear discretizations of k can be used. However, a logarithmic discretization of k is particularly advantageous as it better captures the logarithmic nature of the energy distribution of the waves. An alternative discretization is the equal area method, where $S\left(k\right)\Delta k$ is held constant. A disadvantage of the equal area method occurs at the tail of the distribution where a single sinusoid represents a wider frequency band, leading to possible under-representation of capillary waves.

For each sea state, ten wave profiles of 1 km range are generated using Equation 12, where $S\left(k\right)$ is the ECKV spectrum for a fully developed sea. The spatial frequency k is logarithmically discretized over 1000 bands in the range [0.01 400] rad/m. The upper limit of this range was chosen to be larger than the spatial frequency of a 20 kHz sound in air. The spatial resolution of the sea profile, $\Delta r$, is acoustic frequency dependant, and it is set to be one tenth of the acoustic wavelength. Examples of four sea profiles corresponding to sea states from 2 to 5 are shown in Figure 2.

The mean sea surface elevation is at $z=0$ m. The significant wave heights, $H_s$, of each profile can be measured as $H_s=4{\sigma }_s$, where ${\sigma }_s$ is the standard deviation of the water surface elevation [29]. The average significant wave heights obtained for $U_{10}$ = 3, 5, 7 and 10 m/s found over the ten realizations are $H_s=$ 0.17, 0.47, 0.88 and 1.8 m, respectively. These sets of sea profiles are referred to as sea states 2, 3, 4 and 5, respectively [8]. The slopes of the generated surfaces do not exceed ${30}^{\circ }$ and as such the surfaces are suited for use in a GTPE solver.

To validate the generated profiles, a spectrum can be estimated and compared to the ECKV spectrum. An example of this comparison is shown in Figure 3. It demonstrates good agreement between the generated profiles’ spectra and the ECKV spectrum.

Figure 2. Sea profiles generated from ECKV spectrum for sea states from 2 to 5. Note that x and y-axis are not on the same scale.

Figure 2. Sea profiles generated from ECKV spectrum for sea states from 2 to 5. Note that x and y-axis are not on the same scale.

Figure 3. Comparison between the ECKV spectrum and the averaged spectrum of the generated pseudorandom sea surface profiles, $S_{gen.surf.}$, for $U_{10}$ = 7 m/s. The spectra of the pseudorandom sea surface profiles are averaged over 10 (left), 100 (center) and 1000 (right) realizations.

Figure 3. Comparison between the ECKV spectrum and the averaged spectrum of the generated pseudorandom sea surface profiles, $S_{gen.surf.}$, for $U_{10}$ = 7 m/s. The spectra of the pseudorandom sea surface profiles are averaged over 10 (left), 100 (center) and 1000 (right) realizations.

3. Modeling Refraction

3.1. Monin-Obukhov Similarity Theory

When estimating sound propagation in the atmosphere, refraction effects must be properly modeled. Monin-Obukhov similarity theory (MOST) [30] has been widely used to accurately model wind and temperature profiles in the atmospheric surface layer [31], and has commonly been applied to estimate sound propagation in the atmosphere [11,32,33,34]. The MOST has long been applied to model the marine atmospheric boundary layer (MABL) [35]. According to this theory, wind, temperature and humidity profiles can be obtained from normalized variables: $u_{*}$, $T_{*}$ and $q_{*}$. The wind speed scale is given by the friction velocity, $u_{*}$. The temperature and specific humidity scales, $T_{*}$ and $q_{*}$ are defined as

$$T_{*}=\frac{-Q_H}{{\rho }_0{c}_P{u}_{*}}\mathrm{and}{q}_{*}=\frac{-Q_E}{{\rho }_0{L}_v{u}_{*}}$$

(13)

where $Q_H$ and $Q_E$ are the sensible and latent heat fluxes, ${\rho }_0$ is the density of air, $c_P$ is the specific heat at constant pressure of air and $L_v$ is the latent heat of vaporization for water. Vertical temperature profile, $T\left(z\right)$, can be obtained from MOST as shown in Equation 14.

$$\begin{array}{cc}\hfill T\left(z\right)=& T\left(z_r\right)-\Gamma (z-z_r)+\hfill \\ \hfill & \hfill \frac{P_t{T}_{*}}{{\kappa }_v}\left[ln\left(\frac{z}{z_r}\right)-{\mathrm{\Psi }}_h\left(\frac{z}{L_o}\right)+{\mathrm{\Psi }}_h\left(\frac{z_r}{L_o}\right)\right]\hfill \end{array}$$

(14)

where z is the elevation and $z_r$ is the reference elevation, $\Gamma$ = 0.0098 K/m is the adiabatic lapse rate, $P_t$ = 0.95 is the turbulent Prandtl number in neutral stratification, ${\kappa }_v$ = 0.4 is the Von Kármán constant, ${\mathrm{\Psi }}_h$ is the integral of the universal profile function [11] defined in Equation 20 and in Equation 21 and $L_o$ is the Obukhov length, given by Equation 15

$$L_o=\frac{-u_{*}^3{T}_s{\rho }_0{c}_P}{g{\kappa }_v{Q}_H}$$

(15)

where $T_s$ is the air temperature at the surface.

Within the atmospheric boundary layer, the proportion of turbulence produced by buoyancy as compared to that produced by wind shear is an indicator of the stability of the atmosphere. This proportion can be expressed as the ratio $z/L_o$. For a neutral atmosphere, the effect of buoyancy on density stratification is negligible. In this case, the sensible heat flux $Q_H$ tends to zero as the heat exchange between ground and air is negligible. Thus, $L_o$ will tend toward infinity. Cases in which buoyancy is the driving force can result in either a stable or unstable atmosphere. In either case, $Q_H$ cannot be zero and the magnitude of $L_o$ is therefore small. It also follows that shear effects are relevant only close to the surface. In unstable conditions, the heat flux $Q_H$ from the ground is positive: air parcels tend to rise due to heating from the ground. In stable conditions, it is the opposite: $Q_H$ is negative and air tends to sink. When the stratification is very stable, MOST does not apply. The constant flux assumption is not met as the turbulent mixing is missing. Archer et al. [36] suggests values of 5 m $<L_o<$ 100 m for very stable conditions in the MABL. For heights such that $|z/L_o|\gtrsim 1$, turbulence is mostly suppressed. The Obukhov length can then also be interpreted as the elevation limit for the applicability of MOST [11].

Vertical wind speed profile, $u\left(z\right)$, can be modeled as

$$u\left(z\right)=\frac{u_{*}}{{\kappa }_v}\left[ln\left(\frac{z}{z_0}\right)-{\mathrm{\Psi }}_m\left(\frac{z}{L_o}\right)+{\mathrm{\Psi }}_m\left(\frac{z_0}{L_o}\right)\right]$$

(16)

where ${\mathrm{\Psi }}_m$ is the integral of the universal profile function, defined in Equation 20 and in Equation 21, and $z_0$ is the roughness length. Equation 16 is valid only for heights $z\gg z_0$. The gradient of wind near the surface is primarily driven by the roughness length. Over water, $z_0$ depends on friction velocity according to Equation 17, the Charnock equation [37].

$$z_0={\alpha }_c\frac{u_{*}^2}{g}$$

(17)

where ${\alpha }_c$ is the Charnock parameter. This relationship accounts for an increase in roughness as wave heights grow due to increasing surface stress. Different estimations of ${\alpha }_c$ can be found in the literature, with smaller values for open sea, and slightly larger values for coastal regions [8]. In this work, ${\alpha }_c$ is set equal to 0.0144, as suggested by Garratt [38].

According to MOST, humidity follows the same profile function as temperature. The specific humidity profile, $q\left(z\right)$, is then defined in Equation 18.

$$\begin{array}{cc}\hfill q\left(z\right)=& q\left(z_r\right)+\hfill \\ \hfill \frac{P_t{q}_{*}}{{\kappa }_v}& \hfill \left[ln\left(\frac{z}{z_r}\right)-{\mathrm{\Psi }}_h\left(\frac{z}{L_o}\right)+{\mathrm{\Psi }}_h\left(\frac{z_r}{L_o}\right)\right]\hfill \end{array}$$

(18)

Specific humidity is related to concentration of water vapor, C, as $C=q/(1-q)$. However, since in the atmosphere $C\ll 1$, with typical values of 0.02 in tropical marine atmospheres, the approximation $C\approx q$ is often used[11].

The effective sound speed profile is therefore modeled as [11]

$$c_{eff}\left(z\right)=\sqrt{{\gamma }_{air}{R}_{air}[1+0.511q\left(z\right)]T\left(z\right)}+u\left(z\right)$$

(19)

where ${\gamma }_{air}$ is the ratio of specific heats for dry air and $R_{air}$ is the air gas constant. From Equation 19, it can be noted that the effect of $q\left(z\right)$ on the sound speed profile is significantly smaller than that of temperature, however in this work it will not be neglected.

The universal profile functions in this manuscript follow the notation given by Van Renterghem et al. [8]. For a stable atmosphere[39] ($z/L_o>0$):

$${\mathrm{\Psi }}_m=2ln\left(\frac{1+\chi }{2}\right)+ln\left(\frac{1+{\chi }^2}{2}\right)-2ta{n}^{-1}\chi +\frac{\pi }{2}$$

$${\mathrm{\Psi }}_h=2ln\left(\frac{1+{\chi }^2}{2}\right)$$

(20)

where

$$\chi ={\left(1-16\frac{z}{L_o}\right)}^{1/4}$$

For an unstable atmosphere[40] ($z/L_o<0$):

$${\mathrm{\Psi }}_m={\mathrm{\Psi }}_h=-5\frac{z}{L_o}$$

(21)

3.2. Meteorological Data

To best depict realistic refraction scenarios, parameters needed to estimate sound speed profiles are obtained from the SEAFLUX[41] database. SEAFLUX is a satellite-based dataset of surface turbulent fluxes over the global oceans with 1 hour time resolution and spatial resolution on the order of the kilometer.

Meteorological data from SEAFLUX has been selected by season and time of day. Data has been grouped according to $U_{10}$ wind speed, with bins centered at 3, 5, 7, and 10 m/s, with a tolerance of ± 0.25 m/s. Those bins correspond to sea states 2, 3, 4, and 5. The SEAFLUX parameters corresponding to selected bins are used to then obtain $T\left(z\right)$, $U\left(z\right)$, and $q\left(z\right)$ as shown in Sec. Section 3.1. Equation 19 yields a $c_{eff}$ profile for each case. A ray tracer has been used over 1 km range to determine the highest acoustic path deflected back at the surface by downward refraction, $h_{highestray}$. Instances outside the limit of validity of MOST such that $h_{highestray}>\left|L_o\right|$ are discarded. For each $U_{10}$ bin, the average seasonal diurnal and nocturnal wind, temperature, and specific humidity profiles are obtained from the occurrences that meet all of the criteria. Examples of the average seasonal profiles are shown in Figure 4, along with the corresponding $c_{eff}$ profile.

A location offshore from Duck, NC has been selected to extract meteorological data from SEAFLUX. This location is chosen as data will be used as benchmark for planned measurements in the area. The bulk of field work occurs during summer, therefore summertime meteorological data is considered in this paper. During the day, conditions tend to be very stable at lower wind speeds, with very small positive Obukhov length ($L_o<25$ m) and strongly negative latent heat $Q_H$. The criterion of $h_{highestray}>\left|L_o\right|$ is seldom met at these low wind speeds. Roughly 5% and 30% of the occurrences with $U_{10}=$ 3 and 5 m/s, respectively, do not meet this criterion. Furthermore, MOST is not applicable in very stable conditions. Daytime profiles are for this reason not considered. Sound propagation using night-time meteorological profiles is modeled in this work.

Figure 4. Average temperature, $T\left(z\right)$, wind speed, $U\left(z\right)$, and specific humidity, $q\left(z\right)$ profiles, along with their respective effective sound speed profile, $c_{eff}\left(z\right)$. Profiles are obtained as the average of nightly summer instances of wind speed corresponding to sea states from 2 to 5.

Figure 4. Average temperature, $T\left(z\right)$, wind speed, $U\left(z\right)$, and specific humidity, $q\left(z\right)$ profiles, along with their respective effective sound speed profile, $c_{eff}\left(z\right)$. Profiles are obtained as the average of nightly summer instances of wind speed corresponding to sea states from 2 to 5.

4. GTPE Predictions

Sound pressure level (SPL) predictions are obtained using a GTPE solver for sound propagating in a stratified atmosphere. The sound propagation model can account for surface roughness with local slope less than ${30}^{\circ }$ and refraction. Two sets of GTPE predictions are obtained in this work for every sea state considered. The first set of predictions is obtained with pseudorandom sea surfaces representing sea states from 2 to 5, in a non-refractive atmosphere. The results of these simulations are used to estimate the equivalent impedances presented in Sec. Section 6. A second set of GTPE predictions uses the effective sound speed profiles described in Sec. Section 3 and the identical corresponding pseudorandom sea surfaces used in the first set of predictions. This second set of predictions obtained with a refractive atmosphere are considered as ground truth, and thus used to validate the effective impedances. Examples of excess attenuation predicted in a refractive atmosphere are shown in Figure 5.

Sound pressure level at the receiver, $L_r$, at a distance r from a point source can be found as in Equation 22.

$$L_r=L_s-20{log}_{10}r-L_{\mathrm{abs}}+\Delta L$$

(22)

where $L_s$ is the SPL in free field at 1 m from the source, the term $20{log}_{10}r$ accounts for the spherical spreading, $L_{\mathrm{abs}}$[42] is the air absorption, and $\Delta L$ represents the contribution of all additional mechanisms that can influence sound propagation. This last term, also called excess attenuation, can be expressed as function of refraction caused by wind and temperature profiles, absorption and reflection of the shore, sea roughness, and turbulence. This paper is focused on estimating the contribution of the sea roughness to $\Delta L$. Throughout this work, wind is assumed to blow in the same direction of sound propagation from the source to the receiver. The effect of turbulence is neglected due to its limited contribution on sound propagation in a downward refracting atmosphere.

Sound propagation occurs over water only, and source and receiver are placed at 2 m height from the average surface elevation that is defined as $z=$ 0 m. Sources considered in this work are Gaussian pulses with central frequencies [125, 250, 500, 1000, 2000] Hz. The total length of the numerical domain is 1000 m and the height is one tenth of its length. This range is chosen as it represents the longest GTPE numerical domain with feasible computational times (less than an hour) on the machine used in this work for the highest acoustic frequency (2000 Hz). The domain has been divided into vertical and horizontal steps of one tenth of the acoustic wavelength, $\lambda$. A perfectly matched layer is added outside the domain of interest to eliminate reflections from the upper numerical boundary[10].

Ten sea profiles are generated for each sea state. Sound propagation is modeled above each profile for the frequencies of interest. The excess attenuation $\Delta L$ is predicted along the propagation range ten times for each sea state. A logarithmic average[10] is used to obtain the mean value of $\Delta L$ at each location. Examples of the resulting predictions in a refractive atmosphere over different sea states are shown in Figure 5.

Figure 5 shows that given a sound speed profile $\Delta L$ decreases when the surface is rough. Sea roughness introduces a noise mitigation effect that tends to increase with the sea state. In addition, the standard deviation of $\Delta L$ over ten predictions increases with increasing sea state.

Figure 5. Excess attenuation $\Delta L$ at 500 Hz for $c_{eff}$ profiles corresponding to sea states 2, 3, and 4. See Figure 4. Thick continuous lines are the mean $\Delta L$ values and the shaded areas depict the minimum and maximum predictions. Dashed lines are predictions for smooth surfaces with $c_{eff}$ of sea states 2, 3, and 4.

Figure 5. Excess attenuation $\Delta L$ at 500 Hz for $c_{eff}$ profiles corresponding to sea states 2, 3, and 4. See Figure 4. Thick continuous lines are the mean $\Delta L$ values and the shaded areas depict the minimum and maximum predictions. Dashed lines are predictions for smooth surfaces with $c_{eff}$ of sea states 2, 3, and 4.

5. Equivalent Impedance Methods

To estimate the average noise mitigation effect due to a specific sea state using a GTPE solver is more computationally expensive than a single CNPE simulation. Specifically, the GTPE approach requires a continuous update of the discretized governing equation to account for the change in slope of the ground surface. This update occurs for every horizontal spatial step. In addition, in order to capture the average effect of a specific sea state, the simulation has to be repeated using several pseudorandom sea surfaces. The use of this method to predict SPL in real time over long distances (on the order of kilometers) is impractical. A computationally less intensive alternative is to use an equivalent impedance approach[19]. This method consists of replacing the rough surface with a flat absorbing surface that generates the same effect on the atmospheric sound propagation. In this case, the artificially absorbed acoustic energy is that backscattered by the surface. In specific scenarios of downward refraction, part of this energy may be redirected forward. However, because the PE equation is a one-way wave equation from source-to-receiver [10] this effect cannot be quantified.

Different approaches to estimate the equivalent impedance of a rough sea surface have been adopted by Boulanger and Attenborough[19], Bolin et al. [6,7] and Van Renterghem et al. [8]. The method proposed by Boulanger and Attenborough is the same implemented in the Van Renterghem et al. work. However, the results from Boulanger and Attenborough results are not used for comparison as they have been obtained using surfaces that are not derived from sea state spectra. The results of this work are compared to the experimental observations of Bolin et al. and to the results of Van Renterghem et al.. A brief summary of these two is presented here.

Bolin et al

The authors modeled the sea surface with a frequency-dependent equivalent impedance based on the Delany and Bazley model [18]. Sound propagation predictions obtained with a PE solver were compared to long range acoustic measurements (9 km of open water and 750 m of land). A constant effective flow resistivity of 3×${10}^4$$\mathrm{kPa}\mathrm{s}{\mathrm{m}}^{-2}$ was chosen to represent rough water independently of the sea state[7].

Van Renterghem et al

The authors followed an approach similar to that of Boulanger and Attenborough[19]. A frequency dependent equivalent impedance for sea states 2, 3, and 4 was numerically estimated. Pseudorandom sea surfaces were generated from PM sea spectra. Equivalent impedance was obtained using time domain simulations in which a broad frequency pulse was generated at 0.5 m elevation from the rough surface. The receiver was located at a distance of 5 m at the same elevation as the source. The difference between direct and reflected sound fields was used to estimate a spherical reflection coefficient and thus an effective impedance.

The same experiment described in Van Renterghem et al. has been replicated by the authors of this work with a finite element method solver. Results comparable to those shown in Van Renterghem et al. were found. However, the spectral peak for sea states 2 and larger corresponds to surface wavelengths greater than 5 m. For example, the dominant wavelength for a sea state 2 is about 8 m and for a sea state 5 is about 90 m. Thus, to better capture the contribution of longer sea wavelengths, the same sea surfaces were tested with a 20 m separation between source and receiver. The angle of incidence was kept constant by changing the elevation of source and receiver. This adjustment is also sufficient to avoid acoustic shadowing from waves with heights larger than the source and receiver elevation. Slightly larger values of normalized effective surface impedance, $Z_s$, were found with this modified geometry. It should also be noted that this methodology generates some unphysical results for effective impedance (negative values of $Re\left(Z_s\right)$) with both source-receiver separation distances. The same observation on this method was reported by Boulanger and Attenborough[19]. Negative values of $Re\left(Z_s\right)$ are obtained for surface geometry such that the excess attenuation $\Delta L$ at the receiver is larger than 6 dB and thus the reflection coefficient is larger than unity. As a consequence, the cases that yield a negative real part of the effective impedance must be discarded. However, excess attenuation values larger than 6 dB are physically possible in cases of incoherent scattering. In fact, when the roughness scale is comparable to the acoustic wavelength, more than one ray could reach the receiver simultaneously. As a consequence, the method leads to an underestimation of the effective impedance. In fact, the Boulanger and Attenborough method implemented in the Van Renterghem et al. work is valid for acoustic wavelengths large compared to the roughness scale [19]. This condition is not met at all acoustic frequencies and sea states.

In this work, numerical predictions are used to estimate effective impedance as a function of sea state and frequency. Average excess attenuation predictions for each sea state, at 2 m elevation, were obtained from a GTPE solver with a uniform effective speed of sound as described in Sec. Section 4. The GTPE predictions for each sea state and acoustic frequency were compared to CNPE predictions over a flat absorbing surface in the same non-refractive atmosphere. The flat absorbing surface of the CNPE predictions had varying effective impedances $Z_s$. The real and imaginary parts of the effective impedance were varied in logarithmic steps. The difference between the CNPE predictions and the average of a set of GTPE predictions along the range was calculated. The impedance associated with the minimum difference in excess attenuation was used as the seed for a minimization search algorithm. This two-stage search was needed as the function being minimized contains multiple local minima. Examples of this minimization process are shown in Sec. Section 6. Expressions of effective surface impedance dependent on acoustic frequency and significant wave height were found.

6. Results

Comparisons between the average excess attenuation in a non-refractive atmosphere from a GTPE solver and predictions obtained from a CNPE solver using the effective impedance found with the presented method are shown in Figure 6. The excess attenuation found using an equivalent impedance are in good agreement with those found using the GTPE at all frequencies and sea states. The effective impedances presented in this work, normalized to the acoustic impedance of air, are shown in the legend of Figure 6.

6.1. Surface Impedance Parameterization

The effective impedance of rough sea is frequency and sea state dependent. It can be parameterized as

$$\begin{array}{ccc}\hfill Z_s& =& \hfill \frac{a_r\left(Hs\right)}{f^{m_r}}-b_r\left(Hs\right)+i[17.94-4.43ln\left(Hs\right)]\hfill \\ \hfill m_r& =& 0.38ln\left(Hs\right)+0.76\hfill \end{array}$$

(23)

where $H_s$ is the significant wave height in meters, f is the acoustic frequency in Hertz, and $a_r$ and $b_r$ are fitted parameters as reported in Table 1. The real part of $Z_s$ in Equation 23 follows the form discussed in previous works [8,19]. The fitted parameters $a_r$ and $b_r$ clearly depend on the significant wave height. However, an effective parameterization of $a_r$ and $b_r$ over the significant wave height was not found, as small variations of the two parameters lead to significant deviations of the surface impedance $Z_s$.

Alternatively, the surface impedance can be parameterized as

$$\begin{array}{cc}\hfill Z_s=& \hfill 830.8-120ln\left(f\right)-83.4ln\left(H_s\right)\hfill \\ & \hfill +i[17.94-4.43ln(Hs\left)\right]\hfill \end{array}$$

(24)

where the imaginary part is the same as in Equation 23. However, Equation 24 is only valid for values of $Re\left(Z_s\right)$ larger than roughly 50. It should be noted that the imaginary part of $Z_s$ is substantially smaller than the real part. Variability of the imaginary part of $Z_s$ was not effectively captured in the parameterization. Due to comparatively small magnitude of $Im\left(Z_s\right)$, it is assumed constant with frequency in Eqs. 23 and 24. This assumption, if applied outside the frequency range studied in this work may lead to cases where the imaginary part of $Z_s$ exceed the real part, contrary to results reported in this section. Figure 7 shows a comparison between the results of the minimization approach described in Sec. Section 5 and the parameterization summarized by Eqs. 23 and 24.

Figure 8 compares equivalent impedances obtained in this work to those reported by Bolin et al. [6,7] and Van Renterghem et al. [8]. Bolin et al. limited analysis to discrete frequencies 80, 200 and 400 Hz. The value of $Z_s$ corresponding to 80 Hz is not reported in Figure 8, as it is outside the range of frequencies considered in this work. Additionally, the significant wave height $H_s$ was not measured during the experiments and thus the relative sea state is not reported in Figure 8. Real values of effective impedance used in Bolin et al. are comparable to those obtained in this work for low sea state. The value of flow resistivity (3$·{10}^7$ kg/s/m$m^2$) chosen in Bolin et al. corresponds to that of a highly reflective surface[7]. Effective impedance values presented in Bolin et al. are based on the Delany and Bazley model. The model provides an imaginary part of the effective impedance of the same order of magnitude as the real part, as well as non-physical results at low frequencies [43,44]. However, results obtained in this work show that the imaginary part of the effective impedance is generally much smaller than the real part. The study by Van Renterghem et al. reports on sea states from 2 to 4. As shown in Figure 8 the effective impedance estimated by Van Renterghem et al. is substantially smaller than that obtained in this work and those employed by Bolin et al.. This discrepancy may be attributed to the underestimation of the effective surface impedance in Van Renterghem et al., as discussed in Section 5.

Figure 7. Normalized effective surface impedance $Z_s$, real part (Top) and imaginary part (Bottom). The minimization results were obtained from the method of Sec. Section 5. A magnification of the imaginary part is also shown to better capture the $Z_s$ difference between sea states.

Figure 7. Normalized effective surface impedance $Z_s$, real part (Top) and imaginary part (Bottom). The minimization results were obtained from the method of Sec. Section 5. A magnification of the imaginary part is also shown to better capture the $Z_s$ difference between sea states.

Figure 8. Comparison between $Z_s$ obtained in this work for a non-refractive atmospheres along with results obtained by Bolin et al. and Van Renterghem et al.. Note that the y-axis is on a logarithmic scale.

Figure 8. Comparison between $Z_s$ obtained in this work for a non-refractive atmospheres along with results obtained by Bolin et al. and Van Renterghem et al.. Note that the y-axis is on a logarithmic scale.

6.2. Limitations of the Effective Impedance Method for Rough Sea Surfaces

Scattering from a rough surface can be of two types: coherent and incoherent. Coherent scattering, in the direction of specular reflection, is prevalent when the surface roughness scale is small compared to the acoustic wavelength. Incoherent scattering, in non-specular directions, dominates when the roughness scale is comparable to the acoustic wavelengths[33]. For the sea states considered in this work, the roughness varies from decimeters to a few meters and the acoustic wavelengths are in the same range. Therefore, the dominant scattering is incoherent. While coherent scattering can be seen merely as a surface property and thus decoupled from any other effect, as it does not modify the main direction of the reflected rays, this decoupling is not true for incoherent scattering. In fact, the incoherent scattering depends on the source-receiver positioning as well as refraction. For this reason, effective surface impedance models presented in previous works are limited to surfaces that generate coherent scattering only [19,33].

The effective impedance estimated in Van Renterghem et al. as well as those presented in this work, Figure 7, are based on the assumption that the equivalent impedance of the sea surface is a surface property only, and therefore decoupled from refraction. In both cases, impedances are derived from numerical simulations in a non-refractive atmosphere and used to predict excess attenuation in a refractive atmosphere. The limitations of this assumption are studied in this section. A GTPE solver is able to predict excess attenuation $\Delta L$ accounting for surface roughness and refraction without decoupling the two effects. Averaged GTPE results are compared to predictions of a CNPE solver that uses the same refractive atmosphere and a effective surface impedance obtained in a non-refractive atmosphere.

Figure 9 shows excess attenuation predictions from a GTPE solver (solid line) for acoustic frequencies of 250 and 1000 Hz and sea states 2 to 5 with the corresponding $c_{eff}$ profiles shown in Figure 4. Additionally, predictions from a CNPE solver that uses surface effective impedances derived from Equation 23 in the same refractive atmosphere are shown in dotted lines. The two predictions were in good agreement, within 2 dB, for some acoustic frequencies and sea states (see Figure 9, SS2 and SS3 at 250 Hz and SS3 at 1000 Hz). However, discrepancies up to 10 dB were found in other cases (see Figure 9, SS5 at 250 Hz and SS2 at 1000 Hz). The results suggest that the use of an equivalent impedance for a sea surface that generates incoherent scattering can lead to prediction errors. The discrepancies can be overcome by applying a correction factor on the effective impedance. However, this correction factor must be function of sea state, acoustic frequency, source-receiver geometry and refraction. The use of this correction factor improves the accuracy of the predictions, at the expense of the generalization of a surface impedance dependent on sea state and frequency only.

A way to find the correction factor is to re-estimate the surface impedance in a refractive atmosphere, applying the method discussed in Sec. Section 5. The values of the surface impedance $Z_{sR}$ for source and receiver and 2 m elevation and refraction condition shown in Figure 4 can be parameterized over significant wave height, $H_s$, and frequency, f, as in Equation 25.

$$Z_{sR}=\frac{7575}{H_s^{0.91}{f}^{0.78}}+i\frac{1862}{H_s^{0.64}{f}^{0.66}}$$

(25)

The CNPE predictions obtained with the the effective impedances $Z_{sR}$ from Equation 25 are shown in dotted lines in Figure 9. The use of $Z_{sR}$ instead of $Z_s$ increases the accuracy of excess attenuation predictions. However, at high frequency and high sea state, the, e.g. SS5 at 1000 Hz, the equivalent impedance method fails to produce accurate predictions, highlighting the limits of this method for surface roughness exceeding the acoustic wavelength.

Figure 9. Excess attenuation $\Delta L$ for $c_{eff}\left(z\right)$ profiles corresponding to sea states 2, 3, 4, and 5. Solid lines are the mean $\Delta L$ predictions obtained over a rough surface. Dashed and dotted lines are predictions for smooth surfaces with the same $c_{eff}\left(z\right)$ and effective impedance of Equation 23 (non-refractive), and Eq. 25 (refractive), respectively. (Top) Acoustic frequency, 250 Hz, (Bottom) acoustic frequency of 1000 Hz.

Figure 9. Excess attenuation $\Delta L$ for $c_{eff}\left(z\right)$ profiles corresponding to sea states 2, 3, 4, and 5. Solid lines are the mean $\Delta L$ predictions obtained over a rough surface. Dashed and dotted lines are predictions for smooth surfaces with the same $c_{eff}\left(z\right)$ and effective impedance of Equation 23 (non-refractive), and Eq. 25 (refractive), respectively. (Top) Acoustic frequency, 250 Hz, (Bottom) acoustic frequency of 1000 Hz.

7. Conclusions

In this work, a method to generate pseudorandom sea surfaces consistent with a sea spectrum characteristic of a specific sea state is described. The pseudorandom sea surfaces can be incorporated in a GTPE solver to estimate the effect of roughness on atmospheric sound propagation.

A method to estimate acoustic surface impedance for sea states 2 to 5 is shown. Equivalent impedances can be numerically estimated by comparing sound propagation predictions obtained with a GTPE solver to those of a CNPE solver using an equivalent impedance. Effective surface impedances found in this work can be parameterized on acoustic frequency and significant wave height of the sea state.

Effective impedance methods are mostly limited to cases where surface roughness scale is small compared to the acoustic wavelength. In those cases the scattering is coherent, allowing the effective impedance to be considered an intrinsic property of the surface. However, this condition is not often met for the acoustic frequency range considered in this work, for which sea waves often generate incoherent scattering. This work investigates the limitations of using an effective impedance to emulate the effect of sea surface roughness in a refractive atmosphere. Results show that considering the effective impedance as a surface property can lead to prediction discrepancies that can exceed 10 dB. The use of an effective surface impedance correction factor that considers refraction and source-receiver position substantially increases the accuracy of the predictions.