It Often Howls More than It Chugs: Wind versus Ship Noise Under Water in Australia’s Maritime Regions

: Marine soundscapes consist of cumulative contributions by diverse sources of sound grouped into: physical (e.g., wind), biological (e.g., ﬁsh), and anthropogenic (e.g., shipping)—each with unique spatial, temporal, and frequency characteristics. In terms of anthropophony, shipping has been found to be the greatest (ubiquitous and continuous) contributor of low-frequency underwater noise in several northern hemisphere soundscapes. Our aim was to develop a model for ship noise in Australian waters, which could be used by industry and government to manage marine zones, their usage, stressors, and potential impacts. We also modelled wind noise under water to provide context to the contribution of ship noise. The models were validated with underwater recordings from 25 sites. As expected, there was good congruence when shipping or wind were the dominant sources. However, there was less agreement when other anthropogenic or biological sources were present (i.e., primarily marine seismic surveying and whales). Off Australia, pristine marine soundscapes (based on the dominance of natural, biological and physical sound) remain, in particular, near offshore reefs and islands. Strong wind noise dominates along the southern Australian coast. Underwater shipping noise dominates only in certain areas, along the eastern seaboard and on the northwest shelf, close to shipping lanes.


Introduction
The oceans abound with natural physical sounds (from wind, rain, polar ice, and seismic activity), biological sounds (from crustaceans, fishes, and marine mammals), and anthropogenic sounds (from transport, construction, offshore exploration, and mining). Soundscapes naturally change over time because of temporal cycles in weather (e.g., cyclones and annual monsoon [1,2]) and animal behaviour (e.g., diurnal foraging patterns, lunar spawning, seasonal mating, and annual migration [3][4][5][6]). However, in many habitats, soundscapes further change with patterns of human presence (e.g., temporary construction or summer recreation [7]) and some have changed steadily over time with increasing intensity of anthropophony (e.g., due to shipping [8]).
In 1996, the European Commission identified air-borne noise as one of the main terrestrial environmental issues in Europe, having been neglected compared to chemical pollution [9]. Subsequently, the Commission enacted sound mapping as an important step to assess and manage sound exposure levels in urban areas [10]. A little later, the issue of underwater ocean noise received similar attention, being declared a pollutant, and with underwater sound monitoring and mapping being suggested [11]. Nowadays, underwater noise footprints of individual anthropogenic operations are commonly mapped for environmental impact assessments (e.g., [12][13][14]). Longer-term, large-scale marine into five classes based on their length only (and not by type or function; e.g., tanker versus passenger ferry): <25 m, ≥25-<50 m, ≥50-<100 m, ≥100-<200 m, and ≥200 m. The regularity of the AIS location reporting depends on location and time, and the data we used was subsampled to provide locations, at most, every 5 min. From this data, ship tracks were interpolated by dead reckoning to intervals of 60 s if two successive AIS positions met criteria based on the time between the polls and the straightness of the vessel's path (as per Appendix B in [39]). For each ship track, the time spent in each grid cell was computed, and time was summed over all ships belonging to the same class, yielding a grid of cumulative time spent in each cell over the 6-month period, by ship class.

Ship Data
Data on ship type, size, position, and speed were obtained from Automatic Identification System (AIS) logs managed by the Australian Maritime Safety Authority (AMSA). AIS data were extracted for the period 1 April 2015-30 September 2015. Ships were grouped into five classes based on their length only (and not by type or function; e.g., tanker versus passenger ferry): <25 m, ≥25-<50 m, ≥50-<100 m, ≥100-<200 m, and ≥200 m. The regularity of the AIS location reporting depends on location and time, and the data we used was subsampled to provide locations, at most, every 5 min. From this data, ship tracks were interpolated by dead reckoning to intervals of 60 s if two successive AIS positions met criteria based on the time between the polls and the straightness of the vessel's path (as per Appendix B in [39]). For each ship track, the time spent in each grid cell was computed, and time was summed over all ships belonging to the same class, yielding a grid of cumulative time spent in each cell over the 6-month period, by ship class.

Wind Data
Hourly data on surface wind speed (10 m altitude) were obtained over a similar 6-month period (1 April 2012-30 September 2012) from the Bureau of Meteorology and CSIRO [41], based on the NCEP Climate Forecast System [42]. The data varied in spatial resolution (4, 10, and 24 arcminute grids), which we projected and re-sampled to a 5 km × 5 km UTM grid. Over all grid cells, wind speed varied between 0.5 and 30 m/s (i.e., 1-58 kn). Wind speeds were binned to match the sea states represented in the 'Wenz curves' [43] and noise spectra were assigned to each wind speed bin (Figure 1b). The 'Wenz curves' were converted to linear power spectral density, then integrated over frequency, before applying 10log 10 to yield broadband mean-square sound pressure levels. Expressed as root-mean-square sound pressure levels, the associations became: 79 dB re 1 µPa for ≥1-<4 kn, 87 dB re 1 µPa for ≥4-<7 kn, 92 dB re 1 µPa for ≥7-<11 kn, 99 dB re 1 µPa for ≥11-<22 kn, 105 dB re 1 µPa for ≥22-<48 kn, and 113 dB re 1 µPa for ≥48 kn.

Acoustic Zones
The Australian EEZ had previously been broken up into 28 'acoustic zones' (Figure 2), whereby each zone was characterised by a unique set of hydroacoustic parameters of the water (i.e., sound speed profiles), geoacoustic parameters of the seafloor (i.e., thickness of the sediment layer, density, compressional sound speed and attenuation, and shear sound speed and attenuation), and bathymetric parameters (i.e., water depth and slope) [37]. The idea was to set up one sound propagation environment for each zone and then model all of the ships in that zone. ≥1-<4 kn, 87 dB re 1 µPa for ≥4-<7 kn, 92 dB re 1 µPa for ≥7-<11 kn, 99 dB re 1 µPa for ≥11-<22 kn, 105 dB re 1 µPa for ≥22-<48 kn, and 113 dB re 1 µPa for ≥48 kn.

Acoustic Zones
The Australian EEZ had previously been broken up into 28 'acoustic zones' (Figure  2), whereby each zone was characterised by a unique set of hydroacoustic parameters of the water (i.e., sound speed profiles), geoacoustic parameters of the seafloor (i.e., thickness of the sediment layer, density, compressional sound speed and attenuation, and shear sound speed and attenuation), and bathymetric parameters (i.e., water depth and slope) [37]. The idea was to set up one sound propagation environment for each zone and then model all of the ships in that zone.

Source-Receiver Transects
Within each zone, all of the grid cells that had a cumulative ship time > 0 (no matter the ship class) were identified. From each of these 'source cells', 36 radials were cast at 10degree intervals. The bathymetry was extracted along each of these radials in 5 km steps out to a maximum range of 100 km, yielding thirty-six 100 km source-receiver transects around each source. A bathymetry reading at 2.6 km range from the centre of each source cell was added, representing the mean distance between two random points inside a square (i.e., 0.5214 times the edge length [44]). If a transect hit land, all subsequent bathymetry samples were set to 'not a number' along this 100 km transect. If a source sat near a zone boundary, then the 100 km transects were extracted with bathymetry from the neighbouring zone or from a 100 km buffer around the outside of the EEZ.
All of the transects from all of the source cells (of all ship classes) within a zone were passed to an unsupervised Kohonen neural network (i.e., a self-organizing map, SOM) with 900 neurons [45] (also see [13], where this SOM was previously used to cluster bathymetry transects). The neural net sorted the transects into 900 groups, based on their bathymetric shape. Further grouping was achieved by k-means clustering allowing for 64 clusters [46]. Cluster centroids were computed as the arithmetic mean of all transects

Source-Receiver Transects
Within each zone, all of the grid cells that had a cumulative ship time > 0 (no matter the ship class) were identified. From each of these 'source cells', 36 radials were cast at 10-degree intervals. The bathymetry was extracted along each of these radials in 5 km steps out to a maximum range of 100 km, yielding thirty-six 100 km source-receiver transects around each source. A bathymetry reading at 2.6 km range from the centre of each source cell was added, representing the mean distance between two random points inside a square (i.e., 0.5214 times the edge length [44]). If a transect hit land, all subsequent bathymetry samples were set to 'not a number' along this 100 km transect. If a source sat near a zone boundary, then the 100 km transects were extracted with bathymetry from the neighbouring zone or from a 100 km buffer around the outside of the EEZ.
All of the transects from all of the source cells (of all ship classes) within a zone were passed to an unsupervised Kohonen neural network (i.e., a self-organizing map, SOM) with 900 neurons [45] (also see [13], where this SOM was previously used to cluster bathymetry transects). The neural net sorted the transects into 900 groups, based on their bathymetric shape. Further grouping was achieved by k-means clustering allowing for 64 clusters [46]. Cluster centroids were computed as the arithmetic mean of all transects within one cluster (see examples for Zone 16 in Figure 3). Sound propagation was modelled along each of the 64 centroids within one zone.
within one cluster (see examples for Zone 16 in Figure 3). Sound propagation was modelled along each of the 64 centroids within one zone. Figure 3. Graphs of all 98,532 source-receiver transects of Zone 16 plotted by cluster (grey), with centroid bathymetries shown in black. X-axes are range (km) and Y-axes are depth below the sea surface (m). Note the changing Y-ranges.

Sound Propagation Model
Sound propagation over each centroid bathymetry was modelled with RAMGeo in AcTUP V2.8 [47] (https://cmst.curtin.edu.au/products/underwater/ accessed on: 27 March 2021) based on zone-specific acoustic environments consisting of three layers: the water column, an unconsolidated surface sediment layer, and a consolidated calcarenite sediment layer as a half space. Water column parameters included the zone's mean sound speed profile [37] and water density profile. Representative temperature and salinity data were extracted from the World Ocean Atlas [48,49] to calculate water densities based on the UNESCO formula for sea water density [50]. Unconsolidated surface sediment layers throughout the EEZ comprised predominantly fine material (silt-sand) with sufficiently low shear wave speeds (<250 m/s) to allow modelling as a fluid. Hence, unconsolidated surface sediment parameters only included the zone-specific layer thickness, compressional sound speed, compressional wave attenuation, and density [37]. Surface sediment layer thickness was estimated as 0.5 m for zones within the sediment-starved carbonate Figure 3. Graphs of all 98,532 source-receiver transects of Zone 16 plotted by cluster (grey), with centroid bathymetries shown in black. X-axes are range (km) and Y-axes are depth below the sea surface (m). Note the changing Y-ranges.

Sound Propagation Model
Sound propagation over each centroid bathymetry was modelled with RAMGeo in AcTUP V2.8 [47] (https://cmst.curtin.edu.au/products/underwater/ accessed on: 27 March 2021) based on zone-specific acoustic environments consisting of three layers: the water column, an unconsolidated surface sediment layer, and a consolidated calcarenite sediment layer as a half space. Water column parameters included the zone's mean sound speed profile [37] and water density profile. Representative temperature and salinity data were extracted from the World Ocean Atlas [48,49] to calculate water densities based on the UNESCO formula for sea water density [50]. Unconsolidated surface sediment layers throughout the EEZ comprised predominantly fine material (silt-sand) with sufficiently low shear wave speeds (<250 m/s) to allow modelling as a fluid. Hence, unconsolidated surface sediment parameters only included the zone-specific layer thickness, compressional sound speed, compressional wave attenuation, and density [37]. Surface sediment layer thickness was estimated as 0.5 m for zones within the sediment-starved carbonate platform [51]. Surface sediment thickness in the remaining zones appears variable [52][53][54][55][56][57], and so was modelled as 2 m.
In contrast to the surface sediment layer, calcarenite acts as an elastic material with a shear sound speed of 1400 m/s resulting in important propagation effects [51]. While RAMGeo is a parabolic equation for fluid seabeds, reasonable results can be obtained by using an equivalent fluid approximation with reflection coefficients representative of the elastic model [58]. The procedure to find an equivalent fluid approximation for each zone included (a) creating an environment with calcarenite as an elastic material (see [51] for geoacoustic properties), (b) creating an environment with calcarenite as a fluid layer starting with a compressional sound speed of 1250 m/s and a compressional wave attenuation of 4.5 dB/λ c , (c) calculating the reflection coefficients for each modelled frequency for both environments with the reflection coefficient model BOUNCE [59], and (d) adjusting the fluid layer parameters until a representative equivalent had been reached ( Figure 4). platform [51]. Surface sediment thickness in the remaining zones appears variable [52][53][54][55][56][57], and so was modelled as 2 m. In contrast to the surface sediment layer, calcarenite acts as an elastic material with a shear sound speed of 1400 m/s resulting in important propagation effects [51]. While RAMGeo is a parabolic equation for fluid seabeds, reasonable results can be obtained by using an equivalent fluid approximation with reflection coefficients representative of the elastic model [58]. The procedure to find an equivalent fluid approximation for each zone included (a) creating an environment with calcarenite as an elastic material (see [51] for geoacoustic properties), (b) creating an environment with calcarenite as a fluid layer starting with a compressional sound speed of 1250 m/s and a compressional wave attenuation of 4.5 dB/λc, (c) calculating the reflection coefficients for each modelled frequency for both environments with the reflection coefficient model BOUNCE [59], and (d) adjusting the fluid layer parameters until a representative equivalent had been reached ( Figure 4). RAMGeo modelled sound propagation loss for the centre frequencies of eight fulloctave bands between 10 and 2000 Hz over a range of 100 km and up to 7100 m depth, which is more than the maximum water depth (6388 m) of the EEZ. The depth and range resolutions were 10 m. The source depth for all ships was chosen as 5 m below the sea surface. An example RAMGeo output at 640 Hz for the 64 centroid bathymetries of Zone 16 is shown in Figure 5. The bathymetry itself is just visible as a black line, below which propagation loss was greatest. The colours vary from 60 dB propagation loss (dark red) to 110 dB (dark blue). Several patterns are obvious: Convergence zones appeared over all deep bathymetries, leading to low propagation loss (i.e., high received levels) near the sea surface about 60 km from the source (i.e., clusters 6, 11, 26, 32, 41, 52, 53, 55, 59, 60, and 64). Over upwards-sloping bathymetries, propagation loss was greatest (e.g., clusters 1, 4, 19, 35, 37, 54, and 58). Over downwards-sloping bathymetries, sound may reflect into the deep-ocean sound channel, which has an axis at about 1 km depth off Australia. Once inside the channel, sound may propagate over vast ranges at very little additional loss because of no further seafloor (and to a lesser extent, sea surface) interactions (i.e., clusters 12, 13, 25, 36, and 49). Finally, RAMGeo does not include frequency-dependent absorption and so this additional loss was applied outside of and after RAMGeo [60].  RAMGeo modelled sound propagation loss for the centre frequencies of eight fulloctave bands between 10 and 2000 Hz over a range of 100 km and up to 7100 m depth, which is more than the maximum water depth (6388 m) of the EEZ. The depth and range resolutions were 10 m. The source depth for all ships was chosen as 5 m below the sea surface. An example RAMGeo output at 640 Hz for the 64 centroid bathymetries of Zone 16 is shown in Figure 5. The bathymetry itself is just visible as a black line, below which propagation loss was greatest. The colours vary from 60 dB propagation loss (dark red) to 110 dB (dark blue). Several patterns are obvious: Convergence zones appeared over all deep bathymetries, leading to low propagation loss (i.e., high received levels) near the sea surface about 60 km from the source (i.e., clusters 6, 11, 26, 32, 41, 52, 53, 55, 59, 60, and 64). Over upwards-sloping bathymetries, propagation loss was greatest (e.g., clusters 1, 4, 19, 35, 37, 54, and 58). Over downwards-sloping bathymetries, sound may reflect into the deep-ocean sound channel, which has an axis at about 1 km depth off Australia. Once inside the channel, sound may propagate over vast ranges at very little additional loss because of no further seafloor (and to a lesser extent, sea surface) interactions (i.e., clusters 12, 13, 25, 36, and 49). Finally, RAMGeo does not include frequency-dependent absorption and so this additional loss was applied outside of and after RAMGeo [60].

Accumulation of Received Levels
Within each zone, one ship class was treated at a time. The source cells corresponding to one ship class were found, thirty-six 100 km radials were cast at 10-degree intervals around each source cell, and bathymetry was extracted along each radial and sampled in 5 km steps; the mean distance between a ship and a receiver of 2.6 km within the source cell was inserted at the beginning. In other words, source cells were assigned a received level at 2.6 km. Then, stepping through the source cells for this ship class in this zone, for each of the 36 radial transects, the best matching centroid bathymetry was found. In fact, as the SOM had been trained with all source-receiver transects from all ship classes in this zone, finding the 'best matching centroid' reduced to simply looking up into which cluster this transect originally went. Then, for each frequency modelled, propagation loss (PL)

Accumulation of Received Levels
Within each zone, one ship class was treated at a time. The source cells corresponding to one ship class were found, thirty-six 100 km radials were cast at 10-degree intervals around each source cell, and bathymetry was extracted along each radial and sampled in 5 km steps; the mean distance between a ship and a receiver of 2.6 km within the source cell was inserted at the beginning. In other words, source cells were assigned a received level at 2.6 km. Then, stepping through the source cells for this ship class in this zone, for each of the 36 radial transects, the best matching centroid bathymetry was found. In fact, as the SOM had been trained with all source-receiver transects from all ship classes in this zone, finding the 'best matching centroid' reduced to simply looking up into which cluster this transect originally went. Then, for each frequency modelled, propagation loss (PL) along this centroid was recovered and subtracted from the corresponding octave band source level (SL) plus the cumulative time in dB re 1 s (with T in units of second) of this ship class in this source cell, yielding received levels (RL): RL = SL + 10 log 10 (T) − PL In this equation, SL is a number, the octave band source level expressed as a meansquare pressure level [dB re 1 µPa 2 ] at the modelled frequency. The duration term is also a number [dB re 1 s]. PL [dB], however, is a matrix with values in depth and range. Therefore, RL is a matrix containing received sound exposure levels (SEL) [dB re 1 µPa 2 s] as a function of range and depth. There was one such RL matrix for each frequency.
To change from the polar coordinates (in which sound propagation was modelled) to the Cartesian grid of the EEZ, RL was interpolated to the 5 km grid of the EEZ, at each depth. Including frequency as an additional dimension, this yielded a 4-dimensional matrix of longitude, latitude, depth, and frequency, covering the entire EEZ. The matrix was populated cumulatively by summing sound exposure (i.e., in linear, not logarithmic terms) over all 36 transects about each source cell, and then over all source cells, before taking 10 log 10 again to yield cumulative sound exposure levels (C-SEL).

Ship Noise Map
Broadband sound exposure levels were computed by summing sound exposure over frequency, thereby reducing the matrix to three dimensions, then converting to dB. A further reduction to two dimensions was achieved by finding the maximum sound exposure level over the top 200 m, representing the 'worst case' of exposure for animals that dive over this depth [61]. One such map is presented for each ship class, as well as cumulatively over all five classes. These maps of cumulative sound exposure level were accumulated over six months encompassing the austral winter. They can be read as average mean-square sound pressure level (SPL) maps by subtracting the 6-month duration in dB re 1 s: SPL = C-SEL − 10 log 10 (183 d × 24 h/d × 60 min/h × 60 s/min/s) = C-SEL − 72 dB re 1 s

Wind Map
The wind map was produced by converting the hourly root-mean square sound pressure levels to linear mean-square sound pressures, then integrating over time, and converting back to decibel. 10 log 10 (3600) was added to account for the number of seconds per hour, yielding cumulative sound exposure levels from wind at each cell over the 6-month winter period.

Comparison between Ship and Wind Noise
For comparison, the cumulative sound exposure levels from wind were subtracted from those of ships (summed over all classes) in every grid cell, and plotted, to show in which geographic regions one dominated over the other. We also added the modelled sound exposures from ships and wind, then converted to decibel, to plot the combined ambient noise exposure levels over the 6-month period.

Validation
Archival underwater acoustic recordings from the northwest, west, south, and southeast of Australia were used in an attempt to validate the modelled noise maps. These data were collected by autonomous recorders [62]  Long-term spectral averages (LTSA) were computed in 5 min windows and integrated over frequency (10-2000 Hz) and time (1 April-30 September) to yield cumulative sound exposure. LTSAs were visualised in the software CHORUS [63] to provide an overview of the soundscape and its contributors over multiple weeks to months at a time. Spectrograms with a resolution of 1 s (50% overlap) were computed to zoom into any 5 min sample when the sound sources were not immediately identifiable in the LTSAs. Power spectral density percentile plots (in which the nth percentile gives the power spectral density level exceeded n% of the time, at each frequency) helped identify the dominant contributors to the winter soundscapes [64].

Results
Australia-wide maps of ship noise C-SEL (by class) over the period 1 April 2015-30 September 2015 are shown in Figure 6. Cumulative sound exposure levels over all classes are also plotted. Wind noise C-SEL, the level difference between ship noise C-SEL and wind noise C-SEL, and the combined C-SEL of ship and wind noise are shown in Figure 7. Hyperlinks to the data can be found in the Data Availability section.

Validation
Modelled C-SELs of ship and wind noise are compared to measured C-SELs from the validation datasets in Table 1. There is good agreement (to within 3 dB) between model and measurement noise levels at sites 9 (NW Shelf, WA, Australia), 16 (Bremer Canyon, WA, Australia), and 25 (Tuncurry, NSW, Australia). The former two were dominated by strong wind, the latter by ships. Figure 8A shows almost continuous strong wind at the Bremer Canyon site, a faint Antarctic blue whale (Balaenoptera musculus intermedia) chorus throughout winter, peaking in May, and distant passes of ships. The cumulative energy from wind dominates and is the reason for the good agreement between model and measurement. Figure 8B shows briefer periods of strong wind off Tuncurry and a distant Antarctic blue whale chorus. The dominant feature of this soundscape were numerous passes of ships at close range, and this is the reason for the good agreement between model and measurement at this site.
Disagreement between model and measurement noise levels at other sites was due to unaccounted, additional, non-targeted noise contributions to the soundscape: marine animals and industrial operations. Figure 9 provides an overview of the biological contributors to the soundscape. The stereotypical sounds of Omura's whales (Balaenoptera omurai); Antarctic blue whales; pygmy blue whales (Balaenoptera musculus brevicauda); the unidentified source of the spot call, fin whales (Balaenoptera physalus); dwarf minke whales (Balaenoptera acutorostrata), and humpback whales (Megaptera novaeangliae) have been well described in the literature; as have Australian fish choruses [65][66][67]. These animals dominated the winter soundscapes near islands and reefs (sites 1, 4, 6, 9, 12, 14, 15, 17 and 18). Examples of soundscapes almost free from ships but noisy with animals are shown in Figures 10 and 11. Examples of soundscapes affected by anthropogenic noise are shown in Figure 12. At the time of recording, seismic surveying was the most common anthropogenic source that we did not model.
We were able to determine C-SEL variability over time at nearby sites. Winter recordings at sites 17 (2016)

Discussion
The aim of our project was to develop a model for underwater ship noise in the Australian EEZ that could be used by industry and government to manage marine zones, their usage, stressors, and potential impacts. To put ship noise into context, we also modelled natural noise from wind under water. The models are based on numerous assumptions and involve a lot of averaging in space and time, leading to uncertainty. We therefore validated the models as a whole by comparing modelled sound exposure levels to measured underwater sound levels from 25 acoustic data sets collected over a 12-year span. Agreement was good when the underwater soundscape mostly contained the two sources modelled: ships and wind. Agreement was poorer when sound sources were missed (i.e., not modelled): seismic surveying, whales, and fishes.
Ship presence and movement were based on AIS data from the winter of 2015. Ships logged their positions at irregular time intervals, requiring that we interpolate between successive logs. We applied criteria for speed and direction continuity before straight-line interpolation, and where these were not met, we accepted holes in tracks, leading to an underestimation of ship time in the corresponding cells. We further had very few vessels in the smallest class (<25 m), as these mostly private recreational vessels do not log AIS positions. We therefore clearly underestimated their underwater noise contribution, in particular to coastal soundscapes. In addition, we did not take into account ships just outside of the EEZ and so underestimated noise levels near the EEZ boundary. Given that most AIS data were available for the larger and noisier vessels, we chose a monopole source depth corresponding to larger vessels (5 m) and applied this to all vessels in the model, for simplicity. The introduced uncertainty in modelled received levels is perhaps greater in winter (which we modelled) than summer, given all of our sound speed profiles exhibited a shallow surface duct of variable depth. Accounting for different source depths for the different vessel classes would require modelling sound propagation over the 64 cluster centroids in each zone multiple times, which we did not do, but could be done to improve accuracy. This might be desirable for more localised applications and modelling over smaller areas than the entire EEZ (e.g., regional seismic surveys or coastal developments). Placing the monopole at deeper depth than the propeller depth of small vessels during sound propagation modelling will likely enhance long-range received levels of the smaller, hence quieter, vessels, which are possibly underrepresented in the AIS data, meaning the errors do not add but work in reverse. Finally, the source levels produced by the RANDI model fall within the broadband quartiles reported recently [68]; however, the spectral shapes might differ. MacGillivray and de Jong [69] very recently showed that the RANDI model overpredicted source power spectral density below~250 Hz for bulk carriers, vehicle carriers, tankers, container ships, and cruise ships, yet underpredicted source power spectral density above~250 Hz. This might lead to differential errors in different regions (deep versus shallow water), depending on the efficiency with which sound below and above 250 Hz propagates in each environment. Other studies reported RANDI to overestimate [70,71] or underestimate, particularly above 200 Hz [72]. Underprediction of source levels by the RANDI model might be more common for the smallest vessels, in particular those with powerful motors, such as whale-watching boats and tugs [69,[73][74][75][76].
In terms of the underwater sound propagation model used, the most common source of uncertainty is a lack of data on the seafloor composition and thus, acoustic properties. We used typical values from [51], but geoacoustic properties may vary from place to place. Hydroacoustic data (i.e., temperature, salinity, and sound speed profiles) were missing in some coastal zones and thus required spatial extrapolation. The equivalent fluid model applied is only approximate up to grazing angles of 50 • and thus, more accurate for long-range propagation modelling. Modelling sound propagation only along bathymetry cluster centroids, instead of every source-receiver transect, introduced additional uncertainty. However, with a median water depth of 1809 m for all source cells in the entire EEZ, deviations of individual bathymetries from centroid bathymetries are likely to affect modelled received levels more in shallow and coastal rather than offshore waters. While deviations in bathymetry from cluster centroids may change the pattern of constructive and destructive interference and thus yield rather variable received levels at specific locations in space and depth, there will not be a consistent bias in modelled received levels. Modelling along centroids will lead to both over-and underprediction, depending on range, depth, and frequency. These effects will be important on a fine spatial scale, but disappear on a coarse grid. Finally, the received level depends greatly on receiver depth. We chose to plot maximum received levels over the top 200 m, corresponding to the water layer in which most baleen whales travel. Any receiver depth (or depth range) may, of course, be picked from the model results, corresponding to specific animal depths.
The wind model we used was based on the classic review done by Wenz [43]. Other models, such as the Cato model [77] extend to lower frequencies and thus, yield higher levels (up to 2 dB) in high sea states. The Cato model would reduce the model-versusmeasurement difference (i.e., improve the wind noise prediction) at the wind-dominated sites (9,16).
The map of underwater ship noise was based on AIS data from the year 2015, the map of underwater wind noise was based on wind data from 2012, and the in situ measurements were from various years (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). For a fine-scale model (i.e., small grid size), the exact positions and types of vessels would matter and therefore, validation with measurements from different years might be less successful. However, on a coarse grid, fine-scale variability averages out. For the ship noise map to differ by 3 dB, twice the number of ships (i.e., twice the power) would be needed. We showed close agreement in measured levels over consecutive years at the same sites, except when strong temporary sources occurred in some sets (e.g., industrial exploration) or when more variable, biological sources dominated in some years.
The geographic grid size chosen for the model might affect the received levels in some cells and change the ship-to-wind noise ratio. We modelled on a 5 km × 5 km grid, and so the source cells were assigned a received level at 2.6 km range. A 2.5 km × 2.5 km grid would have a mean receiver range of 1.3 km. If ships are evenly distributed within a 5 km × 5 km cell, then halving the grid size will increase received levels within the source cells by 20 log 10 (2) = 6 dB. The time spent in the source cell, however, will decrease by a factor 4, or, 10 log 10 (4) = 6 dB, making up for the increase in received level (i.e., decrease in propagation range and thus, propagation loss). If ships are unevenly distributed within the larger grid cell, then changing to a finer grid will yield a net increase in modelled received noise levels within source cells. In comparison, the modelled noise levels from wind, being a sheet (rather than monopole) source, will not vary with grid size as wind speed changes on a much larger spatial scale offshore. Therefore, in areas where shipping lanes are well-defined and narrow (<5 km wide), ship noise levels may exceed wind noise levels by more than modelled in this article. Based on our model and its 25-point validation, the Australian EEZ has a higher proportion of natural underwater noise from wind over ship noise than the North Sea and likely other northern hemisphere oceans [18,32]. Part of the Australian marine soundscape appears pristine, if pristine is defined as an absence of anthropogenic noise and a richness of biological noise (see also [78]). We have shown that accurate models of the Australian marine soundscape must include biological sources (i.e., primarily whales and fishes). Natural biological and physical noise ought to be considered in management frameworks to provide context (e.g., for noise management in the Southern Ocean [79]).
Our recommendations for future work include the establishment of a databank of Australian ship source spectra as started by [80], which will allow replacing the RANDI model with monopole source spectra from actual measurements. We have shown that other anthropogenic noise sources cannot be excluded in areas and years where these dominate and their contribution to the marine noise budget should be assessed. Comparing long-term cumulative sound exposure might not be the quantity most useful to managers. Instead, sound energy could be integrated over much shorter time frames and maps of % time above certain management thresholds be plotted [81], which is likely more relevant to biological receptors than an annual or seasonal integral or average. The different sound sources have different acoustic features (e.g., ship and wind noise are continuous, while seismic surveying and pile driving are pulsed) and bioacoustic impact is likely driven by different acoustic quantities (e.g., sound exposure versus peak pressure [82]). Therefore, different quantities will have to be mapped for different types of impact. Moreover, these sources exhibit fundamentally different sound radiation fields, where an underwater explosion is a monopole, a ship is a dipole, pile driving a line source, and wind a sheet source, requiring different modelling approaches.

Data Availability Statement:
A spreadsheet with the acoustic parameters that characterise each zone (i.e., mean winter sound speed profile; mean water depth and slope; sediment thickness; and compressional sound speed, shear sound speed, compressional absorption coefficient, shear absorption coefficient, and mean seafloor density) is available for download, as is a shape file of the spatially separated 28 acoustic zones (see https://tinyurl.com/3webp3pn). A spreadsheet with the sound speed profiles, water density profiles, and geoacoustic properties of the seafloor is available as Supplementary Material. Maps and data of cumulative sound exposure levels (C-SEL) from shipping over all ship classes in the Australian EEZ corresponding to Figure 6 can be found at Seamap Australia (https://seamapaustralia.org/; https://tinyurl.com/ahbs6nwr) and the AODN (https://catalogue.aodn.org.au/geonetwork/srv/eng/metadata.show?uuid=480847b4-b6 92-4112-89ff-0dcef75e3b84). Map and data of modelled wind noise from Figure 7 can also be found at Seamap Australia (https://tinyurl.com/a477j3b9) and the AODN (https://catalogue.aodn.org.au/ geonetwork/srv/eng/metadata.show?uuid=0d3c7edc-463a-4fa0-8039-4d5a779035c3). All sites last accessed on: 30 March 2021.

1.
Beginning with a GIS layer of the Australian marine bathymetry.

2.
Add layers to the grid with ship positions, grouped by ship size (i.e., ship length), yielding one layer per ship class.

3.
Split the EEZ grid into 28 previously determined acoustic zones.

4.
For each zone: a. Find all grid cells that contain ships of any class, cast 36,100 km radials in 10-degree intervals, and extract bathymetry along the radials. b.
Cluster all extracted bathymetries (over all radials around all cells with ships) with a neural network and subsequent k-means into 64 clusters. c.
Compute sound propagation along each cluster centroid, for the centre frequencies of adjacent octave bands. d.
For each ship size class: i. Find the cells that contain ships of this class (source cells), cast 36,100 km radials in 10-degree intervals, and extract bathymetry along the radials. ii.
For each source cell: • For each radial: • Look up into which cluster this radial went; • For each frequency: Retrieve propagation loss as a function of range and depth. Add octave band source level for this ship class. Add cumulative time that a ship of this class spent in this source cell to yield sound exposure level as a function of range and depth. Regrid from polar to Cartesian coordinates.
iii. Accumulate sound exposure over all radials and source cells to yield a 4-d matrix of cumulative sound exposure level as a function of longitude, latitude, depth, and frequency for each ship class.
e. Accumulate sound exposure over all ship classes.
Sum over frequency to yield a 3-d matrix of cumulative sound exposure level as a function of longitude, latitude, and depth. 7.
Pick the maximum cumulative sound exposure level over depth to yield a 2-d map of cumulative sound exposure level versus longitude and latitude.