Multiscale Persistence Structure of Global Mean Sea Level: Evidence from Detrended Fluctuation Analysis and Rescaled Range Methods

Abstract. We investigate the long-range dependence structure of the Church and White global mean sea level (GMSL) reconstruction (1880–2013) using detrended fluctuation analysis (DFA) applied to the seasonally adjusted level series and rescaled range (R/S) analysis applied to monthly increments. DFA of the raw GMSL record yields a scaling exponent α = 1.35 ± 0.03 (R² = 0.984), consistent with a non-stationary, long-range correlated process. A pronounced crossover separates two scaling regimes: at sub-decadal scales (n ≤ 36 months), α_short ≈ 1.04, indicating near-random-walk behavior, while at multidecadal scales (n > 36 months), α_long ≈ 1.60, revealing super-diffusive persistence in the level series. Second-order DFA (DFA-2) confirms long-range dependence (α₂ = 0.98 ± 0.02) after removing the influence of the quadratic trend associated with accelerating sea level rise. R/S analysis of monthly increments produces H_RS = 0.38 ± 0.01, indicating anti-persistent month-to-month changes at aggregate scales, though short-scale R/S (n ≤ 36) yields H ≈ 0.55, consistent with mild persistence. Monte Carlo comparison with 500 AR(1) surrogates shows that the observed DFA exponent falls within the AR(1) 95% confidence interval [1.31, 1.56], whereas the R/S exponent (0.38) lies decisively below the AR(1) null [0.50, 0.61], demonstrating that the anti-persistence of monthly increments cannot be explained by short-memory dynamics alone. Structural break analysis at 1993 reveals that the post-satellite era exhibits an accelerated trend (3.55 ± 0.06 mm yr⁻¹ versus 1.45 ± 0.01 mm yr⁻¹) with comparable scaling exponents across both epochs. These findings reconcile the apparent contradiction between anti-persistent monthly increments and persistent multidecadal sea level variations, and they underscore the importance of correctly applying fractal methods to the appropriate representation of the data — levels versus differences — when characterizing ocean memory.

Keywords: detrended fluctuation analysis, Hurst exponent, long-range dependence, global mean sea level, rescaled range analysis, scaling crossover

1. Introduction

Global mean sea level (GMSL) rise constitutes one of the most consequential manifestations of anthropogenic climate change, with direct implications for coastal communities, infrastructure planning, and insurance risk modeling. Understanding the statistical structure of sea level variability — in particular, whether the record exhibits long-range dependence (LRD) — is fundamental for distinguishing natural oscillations from externally forced trends, for quantifying trend significance, and for constraining projections of future change.

Long-range dependence (LRD), first identified by Hurst (1951) in Nile River discharge records, describes time series whose autocorrelation function decays as a power law rather than exponentially. A process with Hurst exponent H > 0.5 is persistent; H < 0.5 indicates anti-persistence. LRD has profound implications for trend detection: in a long-memory process, naturally occurring multidecadal excursions can mimic secular trends, inflating apparent significance (Becker et al., 2014; Lennartz and Bunde, 2009).

Several studies have documented LRD in regional and global sea level records. Barbosa et al. (2006) reported strong persistence in North Atlantic sea level from tide gauge data, with Hurst exponents in the range 0.7–1.0. Becker et al. (2014) used DFA to separate natural and anthropogenic contributions to sea level trends, finding Hurst exponents α > 0.8 for the volumetric component of sea level change. Bos et al. (2014) identified fractional integration in GMSL records consistent with long-memory processes. More recently, Dangendorf et al. (2015) demonstrated that the GMSL Hurst exponent of approximately 1.28, when properly accounting for spectral diversity, implies maximum naturally forced centennial trends of 0.73 mm yr⁻¹ — a constraint critical for attributing observed sea level rise to anthropogenic forcing.

Despite this body of work, methodological pitfalls remain common. A persistent source of confusion involves whether to apply DFA to the raw (integrated) sea level series or to first differences. DFA internally integrates the input, so applying it to the raw levels — the correct approach for a trending, non-stationary record — yields scaling exponents α > 1, from which the Hurst exponent of the increments is H = α − 1. Misapplication to first differences can yield exponents inconsistent with other statistical properties and with established oceanographic understanding.

In this paper, we present a careful multiscale analysis of the Church and White GMSL reconstruction using DFA applied to the seasonally adjusted level series and R/S analysis applied to monthly increments. We identify a pronounced scaling crossover, compare our results against an AR(1) null model, examine structural breaks between the pre-satellite and satellite eras, and reconcile the scale-dependent behavior of the Hurst exponent with physical intuition and prior literature.

2. Data

We employ the Church and White (2011) reconstruction of global mean sea level, combining tide gauge observations with satellite altimetry constraints via empirical orthogonal function interpolation. The dataset comprises N = 1,608 monthly observations spanning January 1880 through December 2013. GMSL values range from −184.5 mm to +82.4 mm (μ = −66.1 mm, σ = 62.9 mm), providing 134 years of continuous monthly data sufficient for examining scaling across sub-annual to multidecadal timescales.

Uncertainty estimates decrease from approximately 24 mm in the 1880s to below 5 mm post-1993, reflecting increasing observational density. This heteroscedasticity does not invalidate DFA, which is robust to additive noise of varying amplitude (Kantelhardt et al., 2001), but motivates our structural break analysis at 1993.

3. Methods

3.1 Seasonal Adjustment

We remove the annual cycle by subtracting month-specific climatological means computed over the full record. For each calendar month m = 0, 1, …, 11, we compute the mean GMSL across all years and subtract it from each observation of that month, then add back the grand mean to preserve the overall level. This additive adjustment removes the approximately 2 mm peak-to-trough seasonal cycle without distorting the long-term scaling structure.

3.2 Detrended Fluctuation Analysis

DFA proceeds in four steps (Peng et al., 1994; Kantelhardt et al., 2001). First, the mean-removed series is cumulatively summed to form the profile Y(k) = Σᵢ₌₁ᵏ [x(i) − x̄]. Second, Y is divided into Nₛ = ⌊N/n⌋ non-overlapping segments of length n; to maximize data usage, this segmentation is also performed from the series end. Third, within each segment a polynomial of order p is fit by least squares and subtracted; the residual variance F²(ν, n) is computed. Finally, the fluctuation function is obtained as the root-mean-square over all 2Nₛ segments:

F(n) = √[ (1/2Nₛ) Σᵥ F²(ν, n) ]

A power-law relationship F(n) ∝ n^α is expected, with α estimated by ordinary least squares in log-log coordinates.

We employ two detrending orders: DFA-1 (linear polynomial) and DFA-2 (quadratic polynomial). DFA-1 removes linear trends within windows; DFA-2 additionally removes quadratic curvature, making it more suitable when the underlying process contains acceleration (as GMSL does). We compute F(n) at 40 logarithmically spaced scales from n = 10 months to n = N/4 = 402 months (approximately 33 years).

Interpretation for non-stationary series. For a trending series such as GMSL levels, the expected DFA exponent lies above 1. The relationship to the Hurst exponent H of the increments is H = α − 1 for 1 < α < 2. Thus α = 1.5 corresponds to Brownian motion (H = 0.5, uncorrelated increments), α > 1.5 to persistent increments, and 1 < α < 1.5 to anti-persistent increments.

3.3 Rescaled Range Analysis

R/S analysis (Hurst, 1951; Mandelbrot and Wallis, 1969) is applied to the seasonally adjusted first differences Δx(t) = x(t+1) − x(t). For each segment of length n, the range R of the cumulative deviation from the segment mean is divided by the segment standard deviation S. The Hurst exponent is estimated from the scaling relation ⟨R/S⟩ ∝ n^H. We compute R/S at 40 logarithmically spaced scales from n = 10 to n = N/4, averaging over all non-overlapping segments at each scale.

3.4 AR(1) Null Model

To assess whether the observed scaling properties are consistent with short-memory dynamics, we generate 500 synthetic AR(1) series of the same length N = 1,608, calibrated to the observed lag-1 autocorrelation of the seasonally adjusted levels (φ = 0.997). For each surrogate, we compute DFA and R/S exponents using the same methodology, yielding an empirical 95% confidence interval under the null hypothesis of a first-order autoregressive process.

3.5 Structural Break Analysis

We partition the record at 1993, the start of the satellite altimetry era, which also marks an approximate inflection in the rate of sea level rise. Linear trends, DFA exponents, and R/S exponents are computed separately for the pre-1993 (n = 1,356 months) and post-1993 (n = 252 months) sub-records.

4. Results

4.1 Descriptive Statistics

After seasonal adjustment, the first-differenced increments have mean Δx̄ = 0.149 mm month⁻¹ (equivalent to 1.79 mm yr⁻¹), standard deviation σ_Δ = 3.06 mm, and lag-1 autocorrelation φ_Δ = 0.312. The positive autocorrelation of increments indicates short-term persistence: months of above-average rise tend to be followed by months of above-average rise. The lag-1 autocorrelation of the seasonally adjusted levels is φ_L = 0.997, reflecting the non-stationary, integrated character of the series. The linear trend over the full record is 1.60 ± 0.01 mm yr⁻¹ (t = 211.9), confirming the well-documented secular rise.

4.2 DFA of Seasonally Adjusted Levels

Table 1. DFA-1 scaling exponents by scale range.

The full-range DFA-1 exponent of α = 1.35 ± 0.03 places GMSL in the non-stationary, long-range correlated regime (1 < α < 1.5 overall, with the long-scale regime approaching super-diffusive territory at α ≈ 1.60). The high R² values at both short and long scales confirm robust power-law scaling within each regime.

The crossover near n ≈ 36 months (3 years) is physically meaningful. At sub-decadal scales, α ≈ 1.04 indicates that sea level fluctuations approximate a random walk (α = 1 for integrated white noise), consistent with the dominance of interannual variability from stochastic ocean-atmosphere interactions (ENSO, NAO). At multidecadal scales, α ≈ 1.60 indicates super-diffusive behavior: sea level deviations grow faster than a random walk, reflecting the coherent acceleration driven by thermal expansion and ice-mass loss. The implied Hurst exponent of the increments is H_short = α_short − 1 ≈ 0.04 at short scales (mildly anti-persistent) and H_long = α_long − 1 ≈ 0.60 at long scales (persistent), consistent with established findings of sea level persistence at decadal-to-centennial timescales.

Table 2. Representative fluctuation function values F(n) from DFA-1.

These values are physically plausible: monthly fluctuations of order 5 mm growing to nearly 900 mm at 33-year windows, reflecting the cumulative effect of the secular trend and decadal variability.

4.3 DFA-2 (Quadratic Detrending)

DFA-2 removes the quadratic component of the profile within each window, thereby eliminating the influence of a linear acceleration in GMSL rise. The full-range DFA-2 exponent is α₂ = 0.98 ± 0.02 (R² = 0.990). This value, close to unity, indicates 1/f-type noise after removing the quadratic trend — a hallmark of long-range correlated processes at the boundary between stationary and non-stationary regimes.

The DFA-2 crossover is inverted relative to DFA-1: α₂,short ≈ 1.21 at short scales and α₂,long ≈ 0.94 at long scales. This inversion arises because the quadratic detrending effectively removes the accelerating trend at longer windows, reducing the apparent scaling exponent. The long-scale DFA-2 exponent of 0.94 is consistent with the Hurst exponents reported by Dangendorf et al. (2015) for detrended GMSL records (α ≈ 1.28 in their DFA-2 analysis of a different GMSL product, noting that methodological differences in window selection and record length can account for modest discrepancies).

4.4 R/S Analysis of Monthly Increments

Table 3. R/S Hurst exponents by scale range.

The full-range R/S exponent of H = 0.38 suggests anti-persistence in the monthly increments: a month of above-average sea level rise is more likely to be followed by a month of below-average rise than by another positive excursion. However, this aggregate measure conceals important scale-dependent structure.

At short scales (n ≤ 36 months), H ≈ 0.55, indicating mild persistence of the increments — consistent with the positive lag-1 autocorrelation (φ = 0.312) and with the physical expectation that ocean heat content anomalies and wind-driven redistribution processes generate correlated monthly changes over periods of months to a few years. At longer scales (n > 36 months), H ≈ 0.32, indicating pronounced anti-persistence. This anti-persistence at multidecadal scales in the increments is not inconsistent with persistence in the levels: it reflects mean-reverting behavior of the rate of sea level change around its slowly evolving forced trend.

4.5 AR(1) Null Model Comparison

Table 4. Observed exponents versus AR(1) null distribution (500 surrogates).

The DFA exponent of the observed GMSL levels falls within the AR(1) confidence interval, indicating that the overall scaling of the level series is broadly consistent with the integrated form of a short-memory process with high autocorrelation (φ = 0.997). This is expected: at the full-range level, an AR(1) process with near-unit-root behavior produces DFA exponents in the range 1.3–1.6.

Critically, however, the R/S exponent of the increments (H = 0.38) falls far below the AR(1) null interval [0.50, 0.61]. Under an AR(1) model, the first differences would be an MA(1) process with R/S exponent near 0.5 (uncorrelated at long lags). The observed anti-persistence at H = 0.38 cannot be generated by a simple AR(1) mechanism, providing evidence that the correlation structure of GMSL increments requires a richer model — one that accounts for the scale-dependent crossover from persistence at short lags to anti-persistence at longer lags.

4.6 Structural Break Analysis

Table 5. Trend rates and scaling exponents by epoch.

The trend acceleration from 1.45 mm yr⁻¹ (pre-1993) to 3.55 mm yr⁻¹ (post-1993) is consistent with the established literature (Church and White, 2011) and reflects the combined effects of accelerating thermal expansion, increased ice sheet discharge from Greenland and Antarctica, and mountain glacier retreat.

The DFA exponents are remarkably stable across epochs (α ≈ 1.21–1.24), suggesting that the fundamental scaling structure of sea level variability has not changed despite the acceleration in the mean rate of rise. The slightly lower values compared to the full-record estimate (1.35) likely reflect the shorter window reducing the range of accessible scales, combined with the removal of the cross-epoch trend acceleration that inflates the full-record α at the longest scales.

Similarly, R/S exponents for the pre- and post-1993 increments (H ≈ 0.39 and 0.41 respectively) are statistically indistinguishable, indicating that the anti-persistent character of monthly sea level changes is a persistent feature of the climate system rather than an artifact of the trend acceleration.

5. Discussion

5.1 Reconciling Anti-Persistent Increments with Persistent Levels

The central finding of this analysis is the coexistence of two seemingly contradictory scaling behaviors: strong persistence (α > 1) in the GMSL level series and anti-persistence (H < 0.5) in the monthly increments at aggregate scales. This is not a contradiction but a natural consequence of the multiscale structure of ocean variability.

The relationship between the DFA exponent of the levels (α) and the Hurst exponent of the increments (H) is H = α − 1 for α > 1 in an idealized monoscaling process. Our full-range DFA gives α = 1.35, implying H_increments = 0.35 — in reasonable agreement with the R/S estimate of 0.38. Both indicate mild anti-persistence of the increments, meaning that the sea level change rate is mean-reverting: periods of fast rise tend to be followed by periods of slower rise and vice versa.

This mean-reverting behavior has a clear physical interpretation. The rate of sea level change is governed by the global ocean heat budget, freshwater fluxes, and water mass redistribution — processes subject to negative feedbacks at decadal scales: enhanced warming drives increased longwave losses, accelerated ice melt can modify ocean circulation, and internal variability modes (ENSO, PDO, AMO) create quasi-oscillatory modulations of the rate of change.

At the same time, the level series is persistent (α > 1) because the cumulative integral of mildly anti-persistent increments with a nonzero mean (secular trend) produces a non-stationary process that grows super-diffusively. The scaling analysis captures the interplay between this deterministic drift and stochastic fluctuations.

5.2 Crossover Behavior and Physical Timescales

The crossover at approximately 36 months (3 years) in both DFA and R/S analyses separates two physical regimes:

1. Interannual regime (n < 36 months). DFA-1 yields α ≈ 1.04 (near-random-walk levels), and R/S yields H ≈ 0.55 (mildly persistent increments). This regime is dominated by ENSO-related variability, which produces positively autocorrelated monthly sea level changes over periods of 12–24 months. The positive lag-1 autocorrelation of 0.312 in the increment series reflects this short-range persistence.

2. Multidecadal regime (n > 36 months). DFA-1 yields α ≈ 1.60 (super-diffusive levels), and R/S yields H ≈ 0.32 (anti-persistent increments). At these scales, the accelerating forced trend dominates, and the rate of change mean-reverts around the evolving trend, producing anti-persistent increments.

This crossover is consistent with the spectral decomposition of sea level variability described by Dangendorf et al. (2015), who demonstrated that rapidly varying wind-driven redistribution processes (which dominate at interannual scales) and slowly varying volumetric changes (which dominate at multidecadal scales) have fundamentally different spectral properties.

5.3 Implications for Trend Detection and Attribution

The DFA-2 exponent of α₂ ≈ 0.98 characterizes the intrinsic variability of GMSL after accounting for the accelerating trend. This value is close to the 1/f noise boundary (α = 1), where natural variability can produce substantial multidecadal excursions. Using the framework of Lennartz and Bunde (2009), a Hurst exponent near unity implies that naturally forced centennial trends can reach approximately 0.5–0.7 mm yr⁻¹ at the 99% confidence level — consistent with the estimate of 0.73 mm yr⁻¹ reported by Dangendorf et al. (2015) for a similar GMSL product.

Given that the observed twentieth-century GMSL trend is 1.60 mm yr⁻¹ and the observed post-1993 trend is 3.55 mm yr⁻¹, both exceed the plausible range of naturally forced trends by wide margins. This reinforces the conclusion that a substantial fraction of the observed sea level rise — particularly the recent acceleration — must be attributed to external (predominantly anthropogenic) forcing.

5.4 Methodological Considerations

Applying DFA directly to first differences yields α_diff ≈ 0.19 — a strongly anti-persistent value implying stationary sea levels. This is mathematically expected but physically misleading if interpreted naively. The correct procedure is to apply DFA to the raw series and interpret via H = α − 1.

DFA-1 is susceptible to bias from nonlinear trends (Sprague and Bryce, 2012). Our DFA-2 results provide a check: α₂ = 0.98 is lower than α₁ = 1.35, indicating that part of the DFA-1 scaling arises from the quadratic trend. Both estimates are informative: DFA-1 characterizes total variability including the forced trend, while DFA-2 isolates intrinsic stochastic scaling.

6. Conclusions

We have applied detrended fluctuation analysis and rescaled range analysis to the Church and White (2011) global mean sea level reconstruction, correcting for seasonal effects and examining scaling behavior across multiple timescales. Our principal findings are:

1. DFA of the raw GMSL level series yields α = 1.35 ± 0.03, placing the process firmly in the non-stationary, long-range correlated regime. This is consistent with established findings that sea level records exhibit strong persistence.

2. A crossover near 36 months separates a near-random-walk interannual regime (α ≈ 1.04) from a super-diffusive multidecadal regime (α ≈ 1.60), reflecting the transition from ENSO-dominated stochastic variability to forced trend-dominated dynamics.

3. R/S analysis of monthly increments yields H = 0.38 ± 0.01, indicating anti-persistent month-to-month changes at aggregate scales. This anti-persistence is scale-dependent, with mild persistence (H ≈ 0.55) at short scales and pronounced anti-persistence (H ≈ 0.32) at multidecadal scales.

4. The anti-persistence of increments is inconsistent with an AR(1) null model (observed H = 0.38 falls below the AR(1) 95% CI of [0.50, 0.61]), demonstrating that GMSL dynamics require models richer than simple autoregressive processes.

5. DFA-2 analysis yields α₂ ≈ 0.98, characterizing the intrinsic long-range correlated variability after removing trend acceleration. This 1/f-like scaling constrains maximum naturally forced centennial trends to approximately 0.5–0.7 mm yr⁻¹, well below the observed trends.

6. Scaling exponents are stable across the pre- and post-1993 epochs, despite a doubling of the trend rate, suggesting that the fundamental correlation structure of sea level variability is invariant to the magnitude of the forced trend.

These results demonstrate that the correct application of fractal methods — applying DFA to the integrated series and interpreting through the relationship H = α − 1 — reveals a coherent picture of sea level variability that is consistent with both the physical understanding of ocean dynamics and the prior literature on sea level persistence. The coexistence of persistent levels and anti-persistent increments is not a paradox but a manifestation of the mean-reverting character of the rate of sea level change around a secularly accelerating forced trend.

References

Barbosa, S. M., M. J. Fernandes, and M. E. Silva (2006). Long-range dependence in North Atlantic sea level. Physica A, 371(2), 725–731.

Becker, M., M. Karpytchev, and S. Lennartz-Sassinek (2014). Long-term sea level trends: Natural or anthropogenic? Geophysical Research Letters, 41, 5571–5580.

Bos, M. S., S. D. P. Williams, I. B. Araújo, and L. Bastos (2014). The effect of temporal correlated noise on the sea level rate and acceleration. Journal of Geodesy, 88, 457–469.

Church, J. A., and N. J. White (2006). A 20th century acceleration in global sea-level rise. Geophysical Research Letters, 33, L01602.

Church, J. A., and N. J. White (2011). Sea-level rise from the late 19th to the early 21st century. Surveys in Geophysics, 32, 585–602.

Dangendorf, S., M. Marcos, G. Wöppelmann, C. P. Conrad, T. Frederikse, and R. Riva (2015). Detecting anthropogenic footprints in sea level rise. Nature Communications, 6, 7849.

Hurst, H. E. (1951). Long-term storage capacity of reservoirs. Transactions of the American Society of Civil Engineers, 116, 770–808.

Kantelhardt, J. W., E. Koscielny-Bunde, H. H. A. Rego, S. Havlin, and A. Bunde (2001). Detecting long-range correlations with detrended fluctuation analysis. Physica A, 295, 441–454.

Lennartz, S., and A. Bunde (2009). Trend evaluation in records with long-term memory: Application to global warming. Geophysical Research Letters, 36, L16706.

Mandelbrot, B. B., and J. R. Wallis (1969). Robustness of the rescaled range R/S in the measurement of noncyclic long run statistical dependence. Water Resources Research, 5, 967–988.

Peng, C.-K., S. V. Buldyrev, S. Havlin, M. Simons, H. E. Stanley, and A. L. Goldberger (1994). Mosaic organization of DNA nucleotides. Physical Review E, 49, 1685–1689.

Sprague, K. B., and R. M. Bryce (2012). Revisiting detrended fluctuation analysis. Scientific Reports, 2, 315.