A statistical significance test for sea-level variability

A statistical test is presented to address the null hypothesis that sea-level fluctuations in the open ocean are solely due to additive noise in the wind stress. The effect of high-frequency wind-stress variations can be represented as a correlated additive and multiplicative noise (CAM) stochastic model of sea-level variations. The main novel aspect of this article is the estimation of parameters in the discrete CAM model from time series of sea surface height observations. This leads to a statistical test, similar to the red noise [or AR(1)] test for sea surface temperature, which can be used to attribute specific sea-level variability to other effects than wind-stress noise. We demonstrate the performance of this test using altimeter data at several locations in the open ocean.


Introduction
Sea level varies on many temporal and spatial scales (Harrison, 2002). Tide gauges and altimetry are currently providing most of the instrumental data; some of the tide gauge records are over 100 years long, while the altimetry data record starts in 1992. This study is motivated by the problem of attributing sea surface height (SSH) variations to specific mesoscale ocean features, such as eddies (Firing and Merrifield, 2004) or rings (Beal et al., 2011). This is problematic as much of the sea-level variability is due to other processes, for example, high-frequency wind-stress fluctuations and Rossby waves (Hughes and Williams, 2010).
For sea surface temperature (SST) variability, a similar problem occurs of distinguishing SST changes due to specific large-scale phenomena (e.g. El Niño) from those caused by high-frequency fluctuations in atmospheric temperature and the associated heat fluxes. This problem was addressed over 40 years ago by Hasselmann (1976) who introduced a stochastic model of the ocean mixed layer. In this model, the SST anomalies can be modelled by an Ornstein-Uhlenbeck process with a Gaussian probability density function (PDF). The discrete variant of this process is the AR(1) or red noise process, which is serving as a null hypothesis for SST variability. Indeed, when applied to SST variability in the Eastern Tropical Pacific, the El Niño variations are such that this null hypothesis can be rejected (Neelin et al., 1998).
In a series of studies, extensive statistical analyses of daily observed SST (Sura et al., 2006;Sardeshmukh and Sura, 2009) and SSH (Sura and Gille, 2010) variations were performed. There are many areas of the globe where the skewness (S) and excess kurtosis (K) values of SST anomaly time series are far from Gaussian values (S = K = 0). Examples are large areas in the Southern Ocean (e.g. S ~ −0.4) and the Eastern Tropical Pacific (e.g. S ~ −0.4), with also substantial seasonal dependence. Much larger deviations from Gaussian distributions are found in SSH anomaly time series (Sura and Gille, 2010), where S values range from −2 to +2. Analysis of PDFs for SSH indicates that these have a piecewise power law behaviour, also suggested in earlier studies (Thompson and Demirov, 2006). In Hughes and Williams (2010), it was found (using 12 years of altimetry data) that an AR(5) spectral fit provides an adequate representation of the shape of the spectrum over much of the ocean; this information was used for trend (rate of sealevel change) detection. In Bos et al. (2014), the effect of the specific statistical models on the error in the estimated sea-level trend was calculated.
To explain the non-Gaussian behaviour of SST variability, the stochastic Hasselmann model (from now on, red noise model) was extended to include (in addition to the additive noise due to the heat flux) a multiplicative noise term, due to the dependence of the heat flux noise on the SST anomaly (Sura et al., 2006, Sardeshmukh andSura, 2009). The resulting model of correlated additive and multiplicative (CAM) noise can be analysed through its equilibrium PDF by solving the corresponding Fokker-Planck equation. It was shown (Sardeshmukh and Sura, 2009) that this model gives a power law PDF for which / K S 3 2 2 ≥ , in agreement to what is found from SST observations. A CAM noise stochastic model was also used to explain the non-Gaussian behaviour of SSH anomalies (Sura and Gille, 2010) and to show why From all these studies, it is clear that the CAM noise stochastic model may serve as a null hypothesis, in particular for SSH variability where the deviations from Gaussian behaviour are much larger than that for SST. While the use of such a null hypothesis is mentioned in Sura and Gille (2010) this has, to our knowledge, not yet resulted in a statistical test that can be used to reject it. Such a test would replace the red noise test that is often carried out to draw conclusions about the significance of certain frequencies of SSH variability.
The main purpose of this study is to develop a statistical test for wind-stress-driven sea-level variability based on a stochastic model for which the parameters can be estimated from the altimetry (or tide gauge) time series (like the lag-1 autocorrelation in the AR(1) process). In section 2, we present a reduction of the stochastic shallow-water model to the CAM noise model as an alternative to the reduction (Sura and Gille, 2010). Then, in section 3, we further analyse the properties of the CAM model for noise-driven SSH variability. In section 4, we develop a new statistical test for which parameters can be estimated from observed time series. A summary and discussion (section 5) concludes the article.

An alternative motivation for the CAM model
A model that has been traditionally used to study the propagation of surface waves in the ocean is a barotropic shallow-water model in a spherical sector basin with a flat bottom. It is well-known that such a model contains the whole spectrum of waves, e.g. Poincaré and Rossby waves as well as Kelvin waves along the boundaries. It should, therefore, serve as an elementary model to study wind-stress noise-induced variability in sea-level variations.
In Sura and Gille (2010), the CAM model for sea-level variability was motivated by considering small-scale nonlinear (advective) interactions in the shallow-water model. However, the precise interactions leading to the multiplicative noise component are unclear. Below, an alternative (or additional) motivation for the CAM noise model is considered just through the wind-stress variations.

The stochastic shallow-water model
Consider an ocean basin with a horizontal domain bounded by a closed contour G. The density of the ocean is constant, and the flow is driven by a wind stress ( , ) 0 τ φ θ τ τ τ * = φ θ , where 0 τ is its amplitude and ( , ) τ τ φ θ provides the spatial pattern. We denote (U, V) as the horizontal velocities and h as the depth of the basin. The governing shallow-water equations are non-dimensionalized using scales r 0 , H, U, r 0 /U and 0 τ for length, layer depth, velocity, time and wind stress, respectively, where r 0 is the radius of the earth and become (for a flat-bottom ocean) where the dimensionless parameters are given by

Rossby number
Inverse Froude number Ekman n number At the horizontal boundaries, no-slip boundary conditions are prescribed, i.e., The wind-stress forcing has a deterministic steady part ( , ) d d τ τ φ θ and a stochastic part where the latter, following Sura et al. (2001), can be prescribed in terms of atmospheric horizontal surface velocities U a and V a anomalies. Hence, where σ is the amplitude of the noise. In addition, where z( , ) φ θ is a given spatial dimensionless pattern and the quantities u η and v η represent uncorrelated (white) noise and satisfy The same holds for t s ( , , ), ( , , ) The weight function z( , ) φ θ can, for example, be a Gaussian shape with circular symmetry to the atmospheric variability (Sura et al., 2001). The white-noise structure of the atmospheric forcing is an idealization: it reflects the time-scale separation between the slow 'climate system', represented by the ocean, and the fast 'weather system', represented by the atmosphere (Hasselmann, 1976).

Linearized system
We can write the solutions in the model as When we discretize the system on a two-dimensional grid in ( , ) φ θ , we can write it in the following matrix form is the state vector, with n representing the number of grid points, A, M and B are 3n × 3n matrices, D is a 3n × 3n × 3n third-order tensor and η is a 3n-dimensional vector of uncorrelated white noise with unit intensity. The origin of the multiplicative noise term, involving the tensor D, has a direct physical interpretation: the input of momentum by the wind stress depends on the thickness of the layer.
If we represent the state vector X of the equation (2.8) in the form (u,v,h) (here we drop the primes to simplify the notation), and recalling that the matrix M is a diagonal matrix composed, respectively, of ε for the u and v parts and 1 for the h part, we can rewrite the SDE (2.8) in the form where each ij α is a squared matrix (part of the matrix A) h β and h δ are the only non-zero components, respectively, of the matrix B and the tensor D. The quantities u η and η v are the two independent sub-vectors of the vector η. The term hh h u δ has to be interpreted as h ( ) ( ) h ijk j u k δ η , with summation over repeated indices. On the large scale,  1 ε , while F (1) ε =  , and hence, we can neglect the left-hand side in both the first two equations and rewrite the vectors u and v in terms of h and the noise (assuming invertibility of the block matrix formed by , , uu uv vu α α α and vv α ). Substituting these equations in the time evolution for h, we have with an appropriate definition of the parameters. Thanks to the independence of h u and h v , it is guaranteed that h h is also a white-noise term.

Analysis
As in the model of Hasselmann (1976), we now neglect the interaction between the sea-level height at different locations so capture only the local effects of the wind-stress noise (which are assumed to be the strongest). This is the principle of diagonal dominance of  α in Sardeshmukh and Sura (2009), and hence, we focus on the scalar version of (2.10), given by where now h,  A,  B and  D are scalars. The physical interpretation of the constants can be interpreted as the effect of the local near-geostrophic response on the sea level. The multiplicative noise term arises because the input of momentum by the wind stress depends on the total layer thickness. This is different from the SST model in Hasselmann (1976), where the atmospheric heat flux into the mixed layer does not depend on the total mixed layer depth. Note that a more extended mixed-layer model as in Sardeshmukh and Sura (2009), where also the dependence of the perturbative heat flux on SST is taken into account, does include a multiplicative noise term.
The equation (3.1) captures the effects of the wind-driven noise on sea-level variability in its simplest way. However, without any indication of the appropriate choice of the interpretation of the noise term, we are not able to solve it. In fact, when we try to solve the equation, we have to decide whether to use the Itô or the Stratonovich calculus (Pavliotis, 2016). To make this choice, we have to look at the origin of the noise in the system. As we pointed out in section 2.1, the white-noise structure of the wind forcing is the limit case of a fast decorrelating atmospheric variability. It can be shown (Van Kampen, 1981) that these kinds of systems can be approximated by equations with white noise using the Stratonovich interpretation. Therefore, a more correct way to write the equation (3.1) is the following one where the symbol  is conventionally chosen to represent a Stratonovich SDE and the subscript t is referring to the time dependence of the stochastic process. However, from a mathematical point of view, it is usually easier to deal with Itô SDEs, especially because they constitute Markovian processes, in contrast with Stratonovich SDEs. It is easy to show (Pavliotis, 2016) that the equation (3.2) is equivalent to the Itô SDE where the solution p p h t ( , ) = is the PDF of the stochastic variable h t at time t. Furthermore, the boundary conditions for the distribution and its derivative in space are We are interested in the stationary distribution, that is which can be solved to give and hence the integral of p h ( ) s over the whole domain diverges. By definition, a solution of the Fokker-Planck equation has to be smooth and integrable over the domain (Pavliotis, 2016). Hence, for D > 0, the solution is given by where N is the normalization constant. One can show that p h ( ) s satisfies the stationary Fokker-Planck equation and the requirements for its solution. In case D < 0, the solution is found as In Figure 1a, a series of trajectories is seen to converge to an equilibrium distribution, and the histogram in Figure 1b can be well represented by the analytically determined p s .

Moments of the distribution
It is useful to calculate the moments of the distribution. In principle, we could determine the moments of the stationary distribution, starting from the distribution itself The integrals cannot be solved analytically, but ODEs for the moments can be obtained using an Itô formula. Given a one-dimensional SDE and a smooth function t n → , the following (Itô) formula holds (Pavliotis, 2016) Defining the nth moment as we can write an evolution equation for M n by combining the equations (3.3) and (3.15), resulting in  which corresponds (with a different notation) to those obtained in Sardeshmukh and Sura (2009). Hence, we can compute some interesting statistical quantities relative to CAM noise (mean µ, variance C, skewness S and kurtosis K) As expected, the statistical properties of the solution h t are quite different from the ones obtained from an Ornstein-Uhlenbeck process. For instance, the distribution is asymmetric (skewness different from zero), and the kurtosis (different from three) indicates a heavy-tail behaviour. It is easy to show that, for D = 0, we can recover the corresponding moments for the Ornstein-Uhlenbeck process. Furthermore, in order for the variance to be positive, we have to impose the condition 2A + D 2 < 0 for the SDE (3.3).

Temporal discretization
The easiest way to perform a temporal discretization of the SDE (3.3) is using the Euler-Maruyama scheme (Gardiner, 2009), resulting in with Δt the time step and η t time uncorrelated white noise, with zero-mean unitary variance. The stochastic process (3.21) is equivalent to the equation (3.3), assuming that Δt is small. The same result is valid in the comparison between an AR(1) and an Ornstein-Uhlenbeck process. A simple argument to prove this statement can be given calculating the second and the third moment of the discrete process (3.21) and expressing them in terms of the parameters A, B, D and Δt. In the limit → t 0 ∆ , we can neglect the higher order terms, and the equations turn out to be equivalent to the moments of an SDE with CAM noise (3.19).
Suppose we have a time series of SSH data and want to determine whether the stochastic process (3.21) can describe the statistics of the time series, how do we estimate the parameters in this equation? While the literature concerning the estimation of the parameters of an Ornstein-Uhlenbeck process is quite broad, to our knowledge there is no systematic approach existing for the process (3.21). The first thing to notice is that we are not able to perform the estimation basing it only on the stationary distribution of the corresponding SDE or on its related statistical quantities. In fact, one can always rescale the parameters of the continuous CAM process and also rescale the time. In other words, the distribution is invariant under the following transformation with k > 0.

Parameter estimation
As for the red noise case, to make the estimation unique, it is necessary to include some information about the temporal behaviour of the time series in the procedure. For instance, we could use the autocorrelation structure of the time series. To this effect, we first have to determine the autocorrelation function for an SDE with CAM noise (3.3).
Using the general solution of the SDE (Pavliotis, 2016) ∫ Knowing that the term in front of the exponential function is the variance of the SDE, we can rewrite the normalized autocorrelation function for the corresponding stochastic process as where n represents the time lag. Through the last equation, substituting the autocorrelation function with its estimated value and choosing the value of ∆t (providing that we use the same value in the simulation of the process), it is possible to determine the parameter A. To estimate the other two parameters, we are going to use two different methods: a moments estimation and a fit of the data with the PDF. One of the most well-established methods to estimate the parameters of an Ornstein-Uhlenbeck process consists in writing the autocorrelation function and the second moment (that is the variance, for a zero-mean process) in terms of the parameters, then estimating these two quantities from data and substituting the values in the equations. In our case, since we have an additional parameter (compared with the Ornstein-Uhlenbeck process), we need another equation, namely the third moment. Therefore, after having found A with the autocorrelation function, we are going to make use of the moment equations where the quantities  M 2 and  M 3 are respectively the estimated second and third moment. In Figure 2a, the result of the estimation for a simulated time series is shown (see the figure caption for further details), and it is shown to be successful.
The same data can be used to estimate the parameters A and B via the autocorrelation function and the second moment of the Ornstein-Uhlenbeck process, and the result is shown in Figure 2b. It is clear that the estimated parameters differ substantially from the true ones. The result is a distribution which does not resemble at all the true one. The reason behind this behaviour has to be found in the statistical properties of the stochastic processes under consideration: while the Gaussian distribution is completely determined by its first two moments, this is not the case for more complex processes, such as the SDE with CAM noise (3.3). Therefore, even if the autocorrelation function and the variance of two processes are the same, it does not mean they have the same distribution, assuming that at least one of them is not Gaussian.
Therefore, to make the test of the null hypothesis meaningful, one should consider, together with the power spectrum analysis, an evaluation of the accuracy of the parameters estimation. To this effect, we introduce a χ 2 measure ( ) where z i is the midpoint of the ith interval obtained by binning the data (L being the number of bins), F(z i ) is the corresponding value of the normalized histogram, and p x ( ) s is the stationary PDF under consideration (depending on which model is used for the estimation). The χ measure represents the distance between the distribution of the data and the one predicted by the model. It has to be noticed that this measure is very useful to evaluate which distribution is the closest one to the one obtained from the data. Nevertheless, it cannot be used to compare the accuracy of a certain estimation for different time series: its value, indeed, depends on how much the data are spread around the mean, assuming a constant k. For the simulated time series in Figure 2, we obtain 10 sim RED 4 χ ∼ − , 10 sim CAM 7 χ ∼ − , where the subscript RED refers to red noise, i.e. Ornstein-Uhlenbeck process with a Gaussian distribution. As expected, the comparison between the two measures allows us to state that the CAM estimation works well, whereas the Ornstein-Uhlenbeck process is not an accurate model for our simulated data.
Another method of estimation consists in fitting the histogram of the time series under consideration with a certain distribution (which in the CAM case is represented by equation (3.11)). That means solving the following nonlinear least-squares problem (for unknowns A, B and D)  where z i is the midpoint of the ith interval obtained by binning the data (L being again the total number of bins), F(z i ) is the corresponding value of the normalized histogram, and p A B D z ( , , , ) s is the stationary PDF for the chosen model. To perform the least-squares problem, we used the MATLAB function lsqnonlin, which computes L functions and finds the minimum of the norm (3.28). The choice of the initial condition seems to be crucial for the convergence of the algorithm, especially for a high number of functions: in our code, we performed the optimization starting from the computation of only two functions (corresponding to a fit with a histogram with two bins), choosing a certain initial condition. Later, we incremented the number of bins and chose, as initial condition, the parameters found in the previous step. Then, we iterated the procedure until the maximum number of bins (100) was reached.
In Figure 3, the result of the estimation for the simulated time series under consideration is shown (see the figure caption for further details). Again we can compute the χ measures for the two estimates, which turn out to be very similar to the previous ones, i.e. 10 7 . For long simulated time series, at least for the CAM noise case, the two methods (moments and fit) give similar results; this is not the case for short ones, with possible consequences in the test of the null hypothesis. We consider the distribution fitting method the most effective one to estimate the parameters: by the method, we take into account also the other moments, while with the moment estimation, we are only able to consider the second moment for a red noise process and the second and third moments for a CAM noise process.

Application: sea-level variability
The distribution fitting method of the previous section is applied to SSH anomaly data in this section. The data set we use consists in daily gridded multi-mission SSH anomalies, measured over ~20 years, collected by the Archiving, Validation and Interpretation of Satellite Oceanographic (AVISO) database (http://www.aviso.altimetry.fr/en/data. html; accessed on November 2016).
In Figure 4, we first give a representation of the accuracy of the CAM noise estimation with respect to the red noise one. For every location in the ocean, we performed the two estimations and compared the values of the respective values of the χ 2 measure (see figure caption for further explanation). It is clear that the CAM estimation works (cyan colour) better in many locations of the oceans. It is also interesting to compare this result with the analysis of the skewness of the SSH distributions obtained from data, as determined in Sura and Gille (2010). From the comparison between Figures 4 and 3 in Sura and Gille (2010), we see that, when the distribution of data is skewed, the CAM estimation works well, especially in the cases where the skewness is positive.
From the results of section 3, we can now formulate the null hypothesis for the sea-level variability. Wind-stress (white) noise causes sea-level variability according to the CAM noise process, and if sea-level variability is attributed to different causes, the null hypothesis should be rejected. To test this hypothesis, the following procedure is used, as applied to sea surface anomaly data (mean zero, detrended, seasonal cycle removed): (i) estimate the parameters A, B and D in the CAM noise process, (ii) run many realizations of the stochastic process with the same time step as that in the data and (iii) determine power spectra for each realization. After choosing a significance level p ∼ , plot the p 1 − ∼ quantile of the sampled spectra, for each frequency, together with the power spectrum of the data. If there are peaks extending above the chosen spectrum, it means that the null hypothesis is rejected at a p 100(1 ) − ∼ % significance level. We applied this procedure for several time series, corresponding to different locations of the ocean (see stars in Fig. 4), which represent the most interesting cases.
We first consider a location (label 1) where, according to Figure 4, the sea-level statistics are approximately Gaussian. In Figure 5 , and hence indeed the red noise estimation performs slightly better than the CAM process one. It has to be noticed that, when the distribution of the data is approximately Gaussian, the CAM process that solves equation (3.28) will have a very small (positive) value of the parameter D. This is the reason why the difference between the two estimations, and the correspondent χ measures, is small. The situation is well represented in Figure 5, where the curves relative, respectively, to red and CAM noise almost overlap, both in regard to the distribution and the power spectrum. The peaks extending above the CAM power spectrum (~100-day period) are generated by the mesoscale eddy field (Holland, 1978) in this region.
Next, we look at a region where the distribution of the sea-level anomalies is definitely skewed. In Figure 6, we present the results for a location in the North Pacific Ocean (location 2). The values of c are 0.91 sim RED χ = and 0.14 sim CAM χ = , clearly showing that a CAM noise spectrum is much better than a red noise spectrum: Figure 6b shows that some peaks, which seem to be significant under the red noise test, are not significant under the CAM noise test. The only remaining significant peak is again generated by the mesoscale eddy field in that region. As a third example, we consider the region west of South Africa (location 3). The Agulhas Current affects the South-West Indian Ocean, along the east coast of South Africa and, where it turns back on itself, it periodically releases an eddy into the South Atlantic Gyre. By investigating the variability of sea level in the area around the east coast of South Africa, we expect to observe a peak in the power spectrum, with a period of ~50 days (Beal et al., 2011). Furthermore, the CAM noise model should not be able to represent the dynamics of the sea level, at least in the region of the spectrum where these phenomena occur (Fig. 7). The values of the χ measure are 0.17 sim RED χ = and 0.07 Again, the ratio between the two measures is low, indicating that the two estimations give similar results: the distribution of the data, indeed, does not differ much from a Gaussian distribution. The 50-day peak is significant under both the red noise and CAM noise tests. A final case is a location where both the estimations fail. Observations over the last century indicate that Kuroshio current in the North Pacific Ocean shows bimodal behaviour off the southern coast of Japan (Qiu and Miao, 2000). While the mechanisms leading to such bimodality have not been completely understood, it is clear that they must be part of the internal dynamics of the ocean, and thus they cannot be explained by wind-stress noise. In Figure 8, the result of the estimation of the parameters is shown for one of these time series. The estimation has been performed for both an Ornstein-Uhlenbeck process and a CAM noise SDE. Even without investigating the power spectrum of the time series, we can conclude that our null hypothesis has to be rejected: this example clarifies the importance of the accuracy of the estimations, as a preliminary investigation before proceeding with the null hypothesis test.

Summary and discussion
In this article, we formulated a statistical test associated with the null hypothesis that sea-level variability is only forced by random (white noise) fluctuations in the wind stress. In many studies on sea-level variability, one wants to distinguish the effects of specific flow phenomena, such as eddies, on sea-level variability from the noisy wind-stress driven variability. In Sura and Gille (2010), it was already argued that the basic SSH model (represented by the stochastic process h t representing sea-level anomalies) is given by the Itô SDE with parameters A, B and D and constraint 2A + D 2 < 0 (which implies that A < 0). It was shown here that the choice of a stochastic linear model with CAM noise for the SSH variability also arises naturally from the shallow-water equations under a random wind forcing, when only the local geostrophic flow effect due to the wind-stress forcing is considered. The multiplicative noise term arises because the momentum input by the wind is inversely proportional to the layer depth, which makes the situation more complicated than in the classical Hasselmann model (Hasselmann, 1976), which captures only additive noise processes. Whereas the estimation of parameters in the red noise process is relatively straightforward, the non-Gaussian behaviour of the SSH statistics requires a more sophisticated method of estimation (i.e. distribution fitting) with respect to the usual autocorrelation and moments evaluation. The reason is the fact that the knowledge of a few moments does not necessarily provide reliable information about the statistical distribution of the data, especially in the case the length of the time series is limited. Furthermore, we defined a χ 2 measure, which is an indicator of the accuracy of our estimation of the PDF. For similar reasons, and relevant to the null hypothesis test, we argued that the power spectrum analysis is not enough: a preliminary study of the distribution is necessary to guarantee that the PDF of the chosen stochastic process (with estimated parameters) can resemble the distribution of data. In other words, the evaluation of the accuracy of the estimation tells us whether the estimation with the chosen model is meaningful or not; if yes, then we can proceed with the power spectrum analysis.
The comparison between the χ measures for the different estimations (red and CAM) indicates that, in many locations of the ocean, the stochastic model with red noise is a good approximation. By contrast, for clearly skewed distributions, the CAM model is much more reliable than the red noise one. The comparison between the power spectra tells us that the CAM power spectrum extends always above the red noise one, and the distance between the two depends on the non-Gaussianity of the data.
Once the null hypothesis was chosen, we performed a statistical test to reject it based on SSH data. From the power spectrum analysis of several selected time series, we can conclude that the CAM noise process under investigation can explain most of the variability of the sea level. When there are some peaks extending above the chosen 99% interval of confidence, they are attributable to physical processes different from wind-stress noise. In the Agulhas case, we clearly see a peak corresponding to a period of ~50 days, extending above the spectrum of the CAM process, as expected. Despite the fact that differences with red noise may be small in practice, the CAM noise test is in principle a better test and is also more stringent to attribute sea-level variability to specific phenomena. Although in the Agulhas case considered, this may not matter much, it can certainly matter in more extreme sea-level fluctuations as CAM has a power law distribution and the red noise an exponential one.
It has to be noticed that the null hypothesis test performed in this work has the desired significance level only in the case it is carried out for each single frequency. If one is interested in multiple testing, for instance when the same process is responsible for multiple peaks in the power spectrum, then a higher frequency pointwise confidence level has to be used (Mudelsee, 2010, section 5.2.5.1). The consequence of not employing the correct multiple testing confidence interval is to detect false positives in the significant peaks.
Obviously, wind-stress noise is not the only source of noise induced by small-scale processes as is the case for SST, where noise is not only due to that in the atmospheric heat flux (Sura et al., 2006) and hence may also exhibit non-Gaussian behaviour. For the SSH case, also the heat flux is likely to be important through steric influences, as are small-scale processes in ocean mixing, such as sub-mesoscale processes. This defines a possible future line of study, where more and more complicated stochastic models can serve to define null hypotheses for small-scale induced variability in relevant observables in the climate system.

Declaration
Funding: European Union's Horizon 2020 research and innovation programme for the ITN CRITICS (643073; D.C. and H.D.); the Mathematics of Planet Earth programme of the Dutch Science Foundation (NWO). Ethical approval: none. Conflict of interest: none.