Tipping point analysis of ocean acoustic noise

We apply tipping point analysis to a large record of ocean acoustic data to identify the main components of the acoustic dynamical system and study possible bifurcations and transitions of the system. The analysis is based on a statistical physics framework with stochastic modelling, where we represent the observed data as a composition of deterministic and stochastic components estimated from the data using time-series techniques. We analyse long-term and seasonal trends, system states and acoustic fluctuations to reconstruct a one-dimensional stochastic equation to approximate the acoustic dynamical system. We apply potential analysis to acoustic fluctuations and detect several changes in the system states in the past 14 years. These are most likely caused by climatic phenomena. We analyse trends in sound pressure level within different frequency bands and hypothesize a possible anthropogenic impact on the acoustic environment. The tipping point analysis framework provides insight into the structure of the acoustic data and helps identify its dynamic phenomena, correctly reproducing the probability distribution and scaling properties (power-law correlations) of the time series.


Introduction
The Preparatory Commission of the Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO) has established a global network of underwater hydrophones as a part of its hydroacoustic observations (others being seismic, infrasound, and radionuclide), with the goal of continuous monitoring for possible nuclear explosions (CTBTO 2013). The CTBTO database provides several unique and large oceanic acoustic records, covering more than ten years of continuous recording with a high temporal resolution of 250 Hz. In this article, we study the records of the hydrophone H01W1 at the Cape Leeuwin station. The hydrophone is located at a depth of about 1 km off the south-west shore of Australia. We apply tipping point analysis and identify the main components of this acoustic dynamical system, which we then model, with reconstruction of the probability distribution and 25 scaling properties (power-law correlations) of the observed data. Both the probability distribution and scaling properties are important for ensuring that the model correctly represents the observed data, because probability distribution characterises the range and frequency of time series values, while scaling properties characterise their temporal arrangement (Kantelhardt et al 2002), . 30 Tipping points in climatic subsystems have become a widely publicised topic of high societal interest related to climate change; see, for example, (Lenton et al 2008). Applications of tipping point analysis have been found in geophysics (Livina and Lenton 2007), (Lenton et al 2009), (Livina et al 2010), (Livina et al 2011), , (Lenton et al 2012a), (Lenton et al 2012b), , (Cimatoribus et al 2013), , statistical physics (Vaz 35 Martins et al 2010), , ecology , and structure health monitoring (Livina et al 2014), (Perry et al 2016).
A stochastic model combining deterministic and stochastic components is a powerful yet simple tool for modelling time series of real-world dynamical systems. Given a one-dimensional trajectory of a dynamical system (the recorded time series), the system dynamics can be modelled by the 40 stochastic equation with state variable z and time t: whereż is the time derivative of the system variable z(t), and D and S are deterministic and stochastic components, respectively. Component D(z,t) may be stationary or dynamically changing (for instance, containing long-term and/or periodic trends). Many geophysical variables, which follow 45 seasonal variability, can be approximated by a stochastic model (1): D(z,t) = −U (z,t). If the structure of the potential (the number of potential wells) changes, the tipping point is a bifurcation. If the potential structure remains the same, while the trajectory of the system samples various states, such a tipping point is transitional (Livina et al 2011). The stochastic component, in the simplest case, may be Gaussian white or red noise, with possible multifractality 60 and other nonlinear properties. The tipping point methodology is currently based on the techniques of degenerate fingerprinting and potential analysis, which are described below.

Methodology
The tipping point analysis consists of the following three stages: 1) anticipating (pre-tipping, or analysis of early warning signals), 2) detecting (tipping), and 3) forecasting (post-tipping).

65
Anticipating tipping points (pre-tipping) is based on the effect of slowing down of the dynamics of the system prior to critical behaviour. When a system state becomes unstable and starts a transition to another state, the response to small perturbations becomes slower. This critical slowing down can be detected as increasing autocorrelations (ACF) in the time series (Held and Kleinen 2004).
Alternatively, the short-range scaling exponent of Detrended Fluctuation Analysis (DFA) (Peng et 70 al 1994) may be monitored (up to 100 units, which in the case, for example, of daily data correspond to 3.5 months, see (Livina and Lenton 2007)). The lag-1 autocorrelation is calculated in sliding windows of fixed length (conventionally, half of the series length) or variable length (for uncertainty estimation) along the time series, which produces a curve of an early warning indicator. This indicator describes the structural dynamics of the time series. If the curve of the indicator remains flat and 75 stable, the time series does not experience a critical change (whether bifurcational or transitional). If the indicator rises to a critical value of 1 (the monotonic trend can be estimated, for instance, using Kendall rank correlation), it provides a warning of critical behaviour.
Lag-1 autocorrelation is estimated by fitting an autoregressive process of order 1 (AR1): where η is a Gaussian white noise process of unit variance, σ is the noise level, and c = e −κ∆t is the ACF-indicator with κ the decay rate of perturbations. Then, c → 1 as κ → 0 when a tipping point is approached. In addition, the DFA method utilises built-in detrending of a chosen polynomial order, which allows one to distinguish transitions and bifurcations in the early-warning signals. This can be done by comparing several early-warning indicators, with and without detrending data in sliding For the purposes of potential analysis, the dynamics of the system is approximated by a potential 95 stochastic model with a polynomial U (which, in general, may depend on both state variable z and time): whereż is the time derivative of the system variable z(t), η is Gaussian white noise of unit variance and η is the noise level. In the case of a double-well potential, U can be described by a polynomial 100 of 4th order (assuming its quasi-stationarity, with dependence on the state variable z only): The Fokker-Planck equation for the probability density function p(z,t), has a stationary solution given by The potential can be reconstructed from time series data of the system using the following relation to the probability density function: which means that the empirical probability density p d (kernel distribution) has a number of modes This stochastic approximation of the system structural dynamics has remarkable accuracy for data subsets of length as short as 400 to 500 data points, demonstrating above 90% rate of successful 120 detection, as was shown in an experiment with double-well-potential artificial data (Livina et al 2011). For data subsets of length greater than 1000 data points it correctly detects the structure of the potential with a success rate of over 98%. Because of the data gaps, we interpolate the SPL data to achieve equidistant 1-minute temporal 150 resolution. We then remove the seasonal periodicity by subtracting the averaged seasonal cycle over the 14 years of observation to obtain the fluctuations where S i is the interpolated SPL data, andS i is the mean 1-minute interpolated SPL data. The

Results and discussion
We analyse the global trends of these five datasets, assuming the simplest linear model in a leastsquare regression. To estimate the uncertainty in the trends, we apply the "jackknife" technique (see (Efron 1982) and (Wu 1986)). We use the "delete-d" variation of the method, with random where d is the number of the excluded datapoints in each sample ("delete-d"), r = n − d, T are statistics of the trend estimator, see further details in (Shao and Tu 1995).

165
The resulting trends show a small annual decline in SPL for all five datasets, as shown in Figure 2. The above trend analysis was applied to deseasonalised fluctuations (SPL broadband). It is interesting that the average annual cycle of the initial broadband data, too, has declining trend, which is illustrated in Figure 3. We next apply the pre-tipping analysis (early warning signals) to analyse lag-1 autocorrelations and variance of the broadband SPL record, with estimation of uncertainty. We vary the length of the sliding windows for calculating these indicators between 1/4 and 3/4 of the record length to obtain the averaged curves and standard uncertainties and display the indicator values at the end of each 185 window, as shown in Figure 5. Similarly to CTBTO data, the effect of increasing autocorrelation and decreasing variance was earlier observed in bifurcating artificial data changing from white noise to random walk, in . The acoustics dynamics may be undergoing a similar tipping.
Note that this analysis of early warning signals is performed with large enough windows (starting from length of 3 years up to 9 years), which identify large-scale variability, with possible dynamics 190 on the scale of decades ahead.
Further, we apply potential analysis to identify smaller-scale variability, varying the length of the sliding window from three days to one year. The resulting potential plot is shown in Figure 5.
El Niño/Southern Oscillation (ENSO) can be monitored using several indices. We show in Fig-8 Nonlin  To understand better what dynamical changes occur in the acoustic fluctuations, it is useful to plot 200 the histograms of the corresponding subsets of data. Figure 6 demonstrates the difference between the 2-well-potential (first part of the year 2016) and the 3-well-potential (second part of the year 2016) subsets, which correspond to green and cyan areas in Figure 5.
The variability of the potential can be understood as appearance and disappearance of the SPL fluctuations, which are present in the 3-well-potential subsets and disappear in sub-periods of 2-  will predominantly see seismic waves that impinge close-to-vertically from below (have a small slowness, high apparent velocity across the array) which are subsequently scattered by the wavy sea surface instead of propagating coherently onwards. Hence steep slopes couple better.
The detailed analysis of directional acoustic propagation is out of scope of the current paper and 230 may be analysed later elsewhere.
Finally, we analyse the scaling properties of the deseasonalised fluctuations of the broadband SPL to identify the type of noise present in this dynamical system. When the noise is white, the Detrended Fluctuation Analysis (DFA) scaling exponent has value 0.5, whereas red noise has values of the exponent higher than 0.5 (Peng et al 1994), with nonstationary red noise having exponent 235 higher than 1 (random walk has exponent 1.5). The scaling exponent is estimating by fitting the fluctuation curve F (s) ∼ s α in a log-log plot, as shown in Figure 7. When we apply the scaling analysis to the deseasonalised broadband SPL, in both short and long temporal range it has a high exponent (about 0.9), which means that the acoustic fluctuations are stationary red noise, and this is how they should be modelled to represent accurately the stochastic term in Eq. (1).

Model
Based on the above analysis, we can formulate a stochastic model for the acoustic oceanic noise. We adopt an additive model with the following terms: where U (t) is the system potential, T (t) is a long-term linear trend, P (t) is a seasonal trend, and 245 Φ are red-noise fluctuations. In Equation 11, we use time as the main variable of the time series, assuming that only the shape of the potential U is defined by the state variable z as described above.
The parameters of the model (the global trend slope, the amplitude of the seasonality, the coefficients of the potential, the scaling exponent of the fluctuations and the dynamics of its autocorrelation and variance) can be derived from the data and used for simulating artificial data for comparison. Such 250 stochastically-modelled artificial data can be used for a long-term forecast of acoustic data and for testing various hypotheses of the hydroacoustic dynamical system. For simplicity, we illustrate this with a double-well potential term U (t), and further parameters derived from the broadband SPL. We show five subsamples of the SPL broadband data in Figure 8, where observed and modelled time series are compared (left column), as well as their fluctuation curves (right column).

Conclusion
We have applied tipping point analysis and identified deterministic and stochastic components of the ocean acoustic data. We have discovered a possible signature of El Niño in the deep-ocean acoustic data, which is an interesting observation confirmed by both potential analysis and direct estimation of the probability density function of the broadband SPL. Given that the hydrophones are located at 260 depth, and the number of factors influencing the hydroacoustic system in conjunction with the global climate system is large, the investigation of the transitional mechanisms between the surface multiannual phenomena and deep-water acoustic processes may be a subject of a separate paper. The current dynamics of the acoustic fluctuations, which demonstrate slow but steady changes in early warning indicators, gives indication of an upcoming tipping point in this hydroacoustic system, with possible 265 appearance/disappearance of system states, which in this context denote higher/lower SPL fluctuations. Because Cape Leeuwin is a busy shipping junction in world trade, and as trading processes intensify (at the same time requiring more modern ships, with more efficient and less noisy engines), we hypothesise that frequency ranges of the oceanic acoustic noise will be affected unequally, due to multiple factors related to anthropogenic activities. Some frequency bands may decrease in level be- Agency (EMSA 2015).

275
Other potential causes of trends in ocean noise levels include changes in the frequency of other anthropogenic sources such as geophysical surveying, changes in the number and distributions of biological sources such as large cetaceans, changes in natural sources of sound such ice breaking and ice formation, and changes in the ocean environment which may affect the propagation of sound (for example, sea temperature). Whatever the potential causes, our analysis presents a possible 280 approach for monitoring and modelling such processes in future.

Acknowledgement
The