For the convenience of the readers, our basic results are shown in Figure 1. We identified five individual oscillations (upper panel), including a sea-level amplitude of 70 mm (top–bottom [t-b]) of the 18.6-year oscillation caused by the lunar nodal oscillation (LNO), whereas Schmith, Thejll, and Nielsen, with their method, found other spectra and found that the amplitude effect of the LNO would be only on the order of 10 mm (t-b) (i.e. one-seventh of what we found). These differences are neither strange nor inexplicable but are caused by the two fundamentally different methods (see section below).
Figure 1.
Upper panel: The five harmonic oscillations identified by Hansen et al., 2015. Solid black curve shows the directly observable quasi-oscillations in normally applied 19-y mean curves detrended for the 18.6-year oscillation, i.e. the sum of the 74-, 56-, 111-, and 28-year oscillations. Lower panel: Comparison of sea-level curve as found by addition of the five found harmonic oscillations and general sea-level rise. Nineteen-year floating mean of annual means.
Together with a general sea-level rise of 1.18 mm/y, the sum of these five sea-level oscillations constitutes a reconstructed or theoretical sea-level curve of the eastern North Sea to the central Baltic Sea (Figure 1, lower panel), which correlates very well with the observed sea-level changes of the 160-year period (1849–2009), from which 26 long tide gauge time series are available from the eastern North Sea to the central Baltic Sea.
Such identification of oscillators and general trends over 160 years would be of great importance for distinguishing long-term, natural developments from possible, more recent anthropogenic sea-level changes. However, we found that a possible candidate for such anthropogenic development, i.e. the large sea-level rise after 1970, is completely contained by the found small residuals, long-term oscillators, and general trend. Thus, we found that there is (yet) no observable sea-level effect of anthropogenic global warming in the world's best recorded region.
KNOWN SEA-LEVEL EFFECTS OF THE LUNAR NODAL OSCILLATION
Schmith, Thejll, and Nielsen generally found agreement—within uncertainties—with our identified parameters (shown in their table 1), but with the important exception that the sea-level amplitude of the 18.6-year oscillation is maximally one-seventh (i.e. around 10 mm) of our finding (70 mm; t-b). This is important, because the amplitude of the 18.6-year oscillation, in our opinion, is larger than any other oscillation of the system.
However, such small amplitudes of the 18.6-year LNO as suggested by Schmith, Thejll, and Nielsen are strongly contrasting both theoretical and observed amplitudes of the 18.6-year oscillation. It is generally agreed that the theoretical, mean amplitude is around 40–50 mm (t-b) if the globe was completely covered with one ocean (cf. Baart et al., 2012) and that the amplitudes are largest at high latitudes and smallest near the equator. At North Sea–Baltic Sea latitudes, the theoretical amplitude should be on the order of 30–40 mm (t-b; cf. Woodworth, 2012). For Stockholm at the outskirt of our test area, Wróblewski (2001) calculated a theoretical equilibrium amplitude of 14.4 mm (t-b).
However, when turning from theory to measurements Woodworth, Shaw, and Blackmann (1991) found real amplitudes between 88 and 402 mm (t-b) in 13 records around the British Isles corresponding to a mean amplitude of 211 mm (t-b), whereas Yndestad, Turrell, and Ozhigin (2008) found by wavelet analysis an amplitude of 80–90 mm at Aberdeen (84 mm in Woodworth, Shaw, and Blackmann, 1991). In the long Stockholm curve, Yndestad, Turrell, and Ozhigin (2008) found dominant imprints of the 18-, 55-, and 75-year oscillations, as we did. In a satellite-based study of the Pacific, Cherniawsky et al. (2010) found real amplitudes of 30–70 mm (t-b), i.e. 1.5 to 3.5 times higher than Pacific theoretical amplitudes.
Thus, where real amplitudes have been studied, they are considerably larger than anticipated theoretical amplitudes. If the 18.6-year oscillation's amplitude would be only around 10 mm in the North Sea and Baltic region as suggested by Schmith, Thejll, and Nielsen, it would be highly surprising and in contrast to theoretical models, as well as in contrast to the empirical study by Woodworth, Shaw, and Blackmann (1991) of the British Isles, the studies of Yndestad (2006) and Yndestad, Turrell, and Ozhigin (2008) of the North Atlantic region, and the study by Baart et al. (2012) of the Dutch North Sea coast.
Our finding of an amplitude of 70 mm (t-b) is within the center range of other empirical findings of the region; Schmith, Thejll, and Nielsen's proposal is not. Figure 2 shows the annual mean residuals we found for identification of the dominant 18.6-year oscillation.
Figure 2.
Residuals left for identification of the 18.6-y oscillation according to our least residual sine regression iteration method for the lowest possible number of oscillators (five) that in coherence can produce the lowest possible root mean square residual of the theoretically reconstructed and the observed sea-level curves.
KNOWN GENERAL SEA-LEVEL RISE
Identification of a linear or polynomial general sea-level rise over the 160-year record is not a general anticipation, but another result, of the method we applied for least residual sine regression iteration. We have shown that the smallest residuals are produced at a general sea-level rise of 1.16 mm/y if no synchronization of the oscillations is anticipated, and 1.18 mm/y if the oscillations are completely synchronized to rational factors of the 18.6-year LNO. The method's finding of these two parametric values is completely independent of anything else than being an integral part of our least residual regression method. The identified general sea-level trend in the range of 1.1 to 1.2 mm/y is within the center range of what has been found in most other studies of the region (referred to and discussed in Hansen, Aagaard, and Binderup, 2012; Hansen, Aagaard, and Kuijpers, 2015). Consequently, we consider our method's capability to reproduce empirically found trends to support strongly the reliability of our method.
OTHER SPECIFIC COMMENTS
In the beginning of the abstract, Schmith, Thejll, and Nielsen say that we have claimed that our theory is supported by tidal theory. This is not true. We have not written about tidal theory, but about classical wave theory as described in textbooks on wavelets and superposition (e.g., Hubbard, 1996) and on entrainment of waves (frequency and amplitude locking) as reviewed by, e.g., Pikovsky and Rosenblum (2007). Thus, as should be expected from classical wave theory, we discovered that our findings show synchronization to rational factors of the oscillations' period lengths and locking to rational rates of the oscillations' amplitudes (see also Hansen, 2015).
Schmith, Thejll, and Nielsen claim that our method is unnecessarily complicated and that the many iterations are “inviting accumulation of round-off errors,” which may explain our different results. This is not true, because each individual iteration applies the original sea-level data set and do not apply any calculations thereof or recalculated figures. Every iteration performs a modulation of the original data set around the oscillation parameters found by the previous iteration, such that still smaller residuals will occur. This leads to a still more precisely determined set of oscillation parameters, which in the next iteration is modulated on the original data set. This procedure is continued until smaller residuals will not occur by modulation of the original data set around the oscillation parameters found by the previous iteration. Thus, accumulation of round-off errors is not possible.
Schmith, Thejll, and Nielsen claim that our results cannot be reproduced. This is not strange, but a consequence of the fact that they did not apply our method as we did and explained in our detailed method section as well as in several letters. Schmith, Thejll, and Nielsen have not redone our time-consuming iteration procedure but, instead, applied a surrogate method, which they erroneously believe simulates our method. As explained in the next section, our method differs fundamentally from Schmith, Thejll, and Nielsen's method and the statistical limitations adhered to it. Our method and basic idea attempt to identify a theoretical curve by identification of the smallest possible number of harmonic oscillations (and a general polynomial or linear trend) that in coherence will produce the smallest possible residuals (i.e. difference between the observed and the sum-curve of all identified elements).
This is a crucial point, and the two methods are not comparable when more than one oscillation characterizes an oscillatory system.
FUNDAMENTAL METHODOLOGICAL DIFFERENCES
In the present discussion, the readers should primarily be aware that they witness two fundamentally different approaches to reconstruction of sea levels (and any other kind of superpositioned, harmonic oscillations). The method applied by Schmith, Thejll, and Nielsen basically claims—within broad limits of statistical uncertainty—that the parameters of singular sea-level oscillations can initially be identified from continuous time series by harmonic amplitude spectrum analysis. In contrast, our method (as described in great detail in Hansen, Aagaard, and Kuijpers, 2015) claims that the parameters of a truly harmonic oscillation cannot be identified without determination of all superimposed sea-level oscillations and general trends. Thus, our method differs from the often applied spectral analysis method of, e.g., Schmith, Thejll, and Nielsen, which will only be roughly able to identify a fraction of the oscillations searched for and mostly those of clearly different periods, as well as relatively large quasi-oscillations formed by superposition of two or more harmonic oscillations.
Because such more directly observable quasi-oscillations are composed by superposition of two or more individual oscillations of different period lengths, the underlying individual oscillations cannot be distinguished by spectral analysis unless the available time series is many times longer than the period lengths of the longest underlying oscillations.
In the present context this means, e.g., that the 85-year peak found in the spectrum shown in Schmith, Thejll, and Nielsen's figure 1 is a quasi-oscillation produced by superposition of “our” 55.8- and 74.4-year harmonic oscillations, which quasi-oscillation will peak every ca. 80 to 90 years (cf. Hansen, Aagaard, and Kuijpers, 2015, figure 4). The exact period and amplitude of the quasi-oscillation will depend on the actual, exact phases of the two superimposed, harmonic oscillations. If, for instance, these two oscillations were the longest of the system, separation by traditional spectral analysis of the 55.8- and 74.4-year oscillations would require a continuous time series of 220 to 300 years (dependent on the two oscillators' phase distribution in the represented time series) to obtain complete separation of the two oscillators, whereas our method is independent of the phases of oscillations shorter than the represented time series and would therefore only require a continuous time series as long as the longest oscillation (75 y in this example).
By our method of least residual sine regression iteration we wished to find a way to decompose more directly observable quasi-oscillations into their underlying components (Figure 1, upper panel). Thus, we found it is possible to identify five oscillators that constitute an ensemble of significant harmonic oscillations plus one general trend, which in coherence will produce very small residuals (i.e. root mean square = 2.5 mm) and great similarity (correlation = 0.997) with the observed sea-level changes (Figure 1, lower panel).
As stated by Schmith, Thejll, and Nielsen, neither they nor we “know of any parametric approach to correctly evaluate confidence levels with trending and cyclic data.” This implies that although the reconstructed or theoretical sea-level curve correlates extremely well with the observed curve, the general lack in statistical science of methods and possibilities for calculating realistic confidence levels of cyclic data superimposed on a general trend can neither strengthen nor weaken any theory on such matters.
The Monte Carlo simulation applied by Schmith, Thejll, and Nielsen relies on surrogate data and neglects our basic principle that the individual oscillators' residuals must be individually and iteratively regressed in the oscillations' order of significance (amplitude). This means that the result of every iteration must be surveyed continuously and that the order of successive regressions must be changed when the advance of iterations shows change in the oscillations' order of significance. Such built-in surveillance is neither present in Schmith, Thejll, and Nielsen's Monte Carlo simulations nor in their amplitude spectrum analyses and may explain their misleading results.
We wish to add these observations illustrating the capacity of our method: (1) Smaller residuals cannot be produced by adding a sixth oscillation in the period interval of 18.6 to 160 years. (2) Larger residuals are produced by individual removal of one of the five oscillations followed by modulation of the remaining four oscillators or by further modulation of the ultimately identified parameters of the five oscillations. (3) Second- and third-order polynomial modulations with accelerations above 0.0004 mm/y2 of the general trend of 1.18 mm/y would also produce larger residuals. (4) Because the oscillations are iteratively modulated in their successively identified order of significance (amplitude) and individually in the order of period, amplitude, and phase, the reconstructed curve in Figure 1 (lower panel) cannot be produced by a different ensemble of oscillators (number of oscillators and different oscillation parameters) without producing larger residuals. (5) Our method is not able to (i) propose oscillations smaller than the smallest parametrically identified oscillator (i.e. 8 mm, t-b), (ii) distinguish general trends from possible oscillators that are longer than the available time series (i.e. 160 y), (iii) propose parametric changes over time of a general trend beyond the repetition period of entrained oscillators (i.e. 223 y).
CONCLUSION
We conclude that the above facts strongly support our hypothesis. The general lack in statistical science of methods for calculating realistic confidence levels of cyclic data superimposed on a general trend can neither strengthen nor weaken any such hypothesis.
However, our method produces results that are in concordance with empirical findings in all other studies of the region's 18.6-year oscillation. Our method identifies a regional general sea-level rise by 1.16–1.18 mm/y corresponding to the center range of 1.1–1.2 mm/y in other studies of the region. Also, our colleagues' finding of a large 85-year peak strongly supports our finding of quasi-oscillatory peaks every 80–90 years formed by superposition of “our” 55.8- and 74.4-year oscillations. We have shown that the identified oscillations and general trend reproduce the observed sea levels extremely well (correlation 0.997). We also have shown that backward prolongation of the five oscillations will reproduce the large preinstrumental peaks (culminating ca. AD 1790 and ca. AD 1850) originally found by Hansen, Aagaard, and Binderup (2012) in the region's salt marsh shorelines. Moreover, it has been shown that a 19-year cycle is present in these salt marsh shorelines in a 500-year record (AD 1150–1652) of coastal salt production (Hansen, 2010). When, furthermore, the three largest oscillations have identified plausible physical causes (sea-level effects of the LNO, North Atlantic Oscillation, and Atlantic Multidecadal Oscillation), we see no theoretical or empirical reasons for compromising our hypothesis.
