Dangendorf et al. (2021)
- Title:
Data-driven reconstruction reveals large-scale ocean circulation control on coastal sea level
- Key Points:
Regional coastal sea levels are largely dominated by sterodynamic sea-level variations
Steric sea level signals originate from open ocean regions often thousands of kilometers away
Better understanding of sterodynamic sea-level variations is key for robust near-term coastal sea-level estimates
- Keywords:
steric dynamic sea level, coastal sea level, tide gauge, Empirical Orthogonal Functions
- Corresponding author:
Dangendorf
- Citation:
Dangendorf, S., Frederikse, T., Chafik, L., Klinck, J. M., Ezer, T., & Hamlington, B. D. (2021). Data-driven reconstruction reveals large-scale ocean circulation control on coastal sea level. Nature Climate Change, 11(6), 514–520. doi:10.1038/s41558-021-01046-1
- URL:
Abstract
Understanding historical and projected coastal sea-level change is limited because the impact of large-scale ocean dynamics is not well constrained. Here, we use a global set of tide-gauge records over nine regions to analyse the relationship between coastal sea-level variability and open-ocean steric height, related to density fluctuations. Interannual-to-decadal sea-level variability follows open-ocean steric height variations along many coastlines. We extract their common modes of variability and reconstruct coastal sterodynamic sea level, which is due to ocean density and circulation changes, based on steric height observations. Our reconstruction, tested in Earth system models, explains up to 91% of coastal sea-level variability. Combined with barystatic components related to ocean mass change and vertical land motion, the reconstruction also permits closure of the coastal sea-level budget since 1960. We find ocean circulation has dominated coastal sea-level budgets over the past six decades, reinforcing its importance in near-term predictions and coastal planning.
Introduction
The ability to predict and estimate future coastal sea-level changes depends on a thorough understanding of the observational record. Recent studies provide evidence that the global sea-level budget — the comparison of observed sea-level rise with the sum of individual contributions — can be reasonably closed with independent observations from remote sensing techniques and Argo floats since 2005 (refs. 1,2,3) but also with a combination of observations and reanalysis estimates since 1900 (ref. 4). This points to notable progress since the 2013 fifth assessment report of the Intergovernmental Panel on Climate Change [5]. The budget has also been closed regionally for larger areas off the coast since 2005 (refs. 6,7) and at basin scales since the 1950s [4,8]. However, at the coast, where information is most needed, budget assessments have higher uncertainty due to the number of additional processes that cancel out in the global mean. These include: changes in the sterodynamic component [9], due to ocean circulation change and steric (density-driven) expansion; gravitational, rotational and deformational effects accompanying global barystatic (ocean mass-driven) sea-level changes; and vertical land motion induced by glacial isostatic adjustment, tectonics or local subsidence.
For barystatic sea level and vertical land motion, a number of observations and model estimates now provide constraints at decadal [10] to century scales11. Specifically, the barystatic estimates due to ice-mass loss, terrestrial water storage and glacier melting and their representation at individual locations have considerably improved [4,8].
However, the role of steric sea-level changes and their signature along shallow coastlines is still poorly understood [12,13]. The steric component, defined as the depth integral of density changes in the ocean, is the dominant driver of sea-level variability in the (deep) open ocean (Fig. 1a) with mostly minor contributions coming from ocean-bottom pressure14 (Fig. 1b). In contrast, at the coast where the ocean’s depth decreases to zero, the resulting steric contribution vanishes as well (Fig. 1a). Thus, changes in ocean density will necessarily lead to larger expansion and contraction of the water column in the deep ocean than at the shallow coast15, driving a sea-surface height gradient between the two. This gradient needs to be balanced and hence forces changes in ocean-bottom pressure, with the latter becoming the dominant contributor to sea-level variability at the coast (Fig. 1b). However, due to a wide variety of ocean dynamic processes, coastal sea-level variability cannot be directly derived from steric sea-level changes in the nearby open ocean12,16,17 (Fig. 1c). Due to the lack of this direct relation between the open-ocean steric effects and coastal sea level, both the geographical origin and the underlying processes of the resulting coastal sterodynamic sea-level variations (hereafter SDSL) are poorly understood, inhibiting their proper quantification6,18. Particularly challenging in many coastal regions is the presence of strong boundary currents and coastally trapped waves that mediate open-ocean steric effects to the coast and introduce, due to their coast-parallel flow and propagation, substantial alongshore coherence in sea level often extending over thousands of kilometres17. Indeed, sea-level signals at the continental coasts often originate from perturbations at remote locations13,16, hindering the simple across-shelf integration of steric height from the nearby deep ocean as previously described by ref. 12. While global ocean reanalyses should provide a physically consistent solution for simulating coastal SDSL, they currently face challenges such as drifts and biases due to air–sea fluxes and ocean mixing errors, coarse model resolution or the assimilation of observations that become sparse back in time19. As a consequence, they poorly reproduce observations at the coast20,21 (particularly spatial trends; Extended Data Fig. 1). These knowledge and modelling gaps about past SDSL changes are critical to resolve given that sea-level information is most needed at the coast. Together with the improved information about non-steric processes, these factors call for an alternative, and ideally data-driven, approach for estimating SDSL in the coastal zone.
Fig. 1: Relative roles of ocean-bottom pressure and steric components for coastal SDSL. a,b, Shown is the explained SDSL variance through steric height (a) and ocean-bottom pressure (b) in a historical ocean simulation with the CNRM-CM6-1-HR ESM. c, Schematic of the relationship between deep-ocean steric expansion and ocean-bottom pressure changes on the shallow shelf. The exchange is largely controlled by boundary currents and coastally trapped waves.
Here, we reassess coastal sea-level changes in light of increasing data and knowledge of the non-steric contributions and by introducing a two-step statistical framework for estimating coastal SDSL changes. In the first step, we investigate the statistical link between coastal SDSL variability and open-ocean steric height to identify the source regions that communicate with the coast and discuss our current understanding of the physical processes that are involved. In a second step, we use this statistical link to reconstruct coastal SDSL solely on the basis of open-ocean steric height by applying a regression approach based on empirical orthogonal functions (EOFs). Finally, we use the reconstructed SDSL estimates to analyse the coastal sea-level budget for nine selected regions.
Identification of regions of common variability
We start by identifying the geographic source regions of coastal SDSL variability. We select 89 tide-gauge records spanning the period from 1960 to 2012 (the chosen period is constrained by the availability the different datasets; Methods), grouped into nine regions of coherent variability (North Sea, Northwest Atlantic north/south of Cape Hatteras, Northeast Pacific, Hawaii, Japan Sea, West Australia, New Zealand and Western Pacific) and shown in Fig. 2a and Extended Data Fig. 2. These records are first corrected for vertical land motion and barystatic sea-level components to isolate density-related SDSL signals (Methods). Then we follow the principal idea presented in refs. 22,23 for the North Sea and the Northwest Atlantic north of Cape Hatteras and correlate the tide-gauge residuals (hereafter SDSLTG) with steric height time series from improved (via mapping methods and correction schemes for instrumental biases) gridded temperature and salinity fields over the upper 2,000 m (ref. 24) in the adjacent open ocean. A comparison with other steric products can be found in Methods. As source regions, we define all locations showing a significant correlation (P ≤ 0.05) with the corresponding regional mean of SDSLTG.
Fig. 2: Origin and reconstruction of coastal SDSL changes. a, Spatial correlation patterns between the detrended SDSL residuals at tide gauges (SDSLTG, grey dots) and detrended steric height from the open ocean for nine coastal regions. Only significant correlations (P ≤ 0.05) are shown. Note that in some regions, correlation patterns of different regions may overlap. b, Illustration of the reconstruction of SDSLEOF in observations (top) and in validation experiments with ESMs (here illustrated by the CNRM-CM6-1-HR ESM, bottom) for the Northwest Atlantic south of Cape Hatteras. Shadings represent the 95% confidence intervals (CI) from all Monte-Carlo samples. c, Correlation and trend differences between reconstructed SDSLEOF and the a priori known SDSLTG in the validation experiments with 12 ESMs. Also shown is the correlation between SDSLEOF and SDSLTG in actual observations (orange dots). TG, tide gauge; NW, northwest; NE, northeast; NCH, north of Cape Hatteras; SCH, south of Cape Hatteras.
For all regions, significant correlations are widely spread over the oceans and often extend thousands of kilometres away from the tide-gauge sites (Fig. 2a globally and Extended Data Fig. 3 for expanded views into each region). Along eastern boundaries the correlation patterns range from the (sub)tropics to high latitudes with narrowing bands of high correlations along the continental slope moving poleward. In the Atlantic (Extended Data Fig. 3a), steric sea-level variations near the Strait of Gibraltar have recently been linked to both local wind forcing and Atlantic Meridional Overturning Circulation changes16,25, while further north along the Portuguese coastline, sea-level variability is tightly connected to longshore wind forcing13,26,27,28. The along-shelf coherence (Extended Data Fig. 3a) is consistent with the hypothesis that coastally trapped waves communicate SDSL variations from the eastern boundary into the North Sea16,17,28. A similar, but much more pronounced, correlation pattern can be found in the Northeast Pacific, where the largest correlations, r > 0.9, occur around the Equator and off the California coast (Extended Data Fig. 3d). This is consistent with ref. 29, who demonstrated that coastal sea levels along the Californian coastline vary in concert with fluctuations in equatorial trade winds and longshore winds generated around the Aleutian Low. Westerlies in the western and central Pacific generate equatorial Kelvin waves, which first propagate eastward along the Equator and then, after being reflected at the eastern boundary, travel northward as coastal Kelvin waves16,30.
For the western regions in the Atlantic and Pacific, correlations indicate a dynamic connection to the western boundary currents (Kuroshio and Gulf Stream), although with an interesting distinction in the Atlantic north and south of Cape Hatteras31 (Extended Data Fig. 3b,c). North of Cape Hatteras the correlation pattern follows the continental slope into the Labrador Sea and the Subpolar Gyre. This agrees with ref. 23, who provided evidence that Labrador Sea density anomalies (driven by both atmospheric processes in the upper layers and variations in the Deep Western Boundary Current at intermediate and deeper levels) propagate southward as coastally trapped waves32. South of Cape Hatteras, however, highest correlations are directly centred on the Gulf Stream pathway. This is consistent with the suggestion that coastal sea-level variability in this region is linked to a fast barotropic response of the coastal ocean to (overturning-related) large-scale heat divergence in the open ocean33,34 and/or subtle variations in the strength and position of the Gulf Stream35.
In the Indo-Pacific regions (West Australia and Western Pacific), coastal SDSL variations are highly coherent even between widely separated regions leading to overlapping correlation patterns extending from the Bay of Bengal into the Central Pacific (Extended Data Fig. 3g,i). These signals are linked to the El Niño/Southern Oscillation and primarily driven by equatorial trade winds36. These trade winds induce westward propagating equatorial Rossby waves, which first cross the Indonesian throughflow region and then travel poleward as coastally trapped waves along the west coast of Australia37,38. The coasts of New Zealand show coherence with the larger South Pacific Gyre region as well as the warm subtropical currents in the Tasman Sea extending onto the southwestern Australian continental shelves. Sasaki39 used an eddy-resolving ocean general circulation model to show that long baroclinic Rossby waves, forced by wind stress curl over the subtropical gyre, are an important driver of New Zealand’s decadal sea-level variability. For gauges at Hawaii, a typical example of an open-ocean site surrounded by steeply sloping topography, highest correlations are pronounced locally northeast of the islands, indicating a more local control of coastal sea-level variations.
In summary, for all nine regions the identified correlation patterns are consistent with known processes induced by ocean dynamics (like wave guides and propagation characteristics) that have also been identified in basic ocean simulation experiments16. This underpins the physical nature of the statistical relationship and suggests that coastal SDSL changes are primarily remotely forced by perturbations in the open ocean and transferred to the coast (Fig. 1c) through the action of (coastally trapped) Kelvin and Rossby waves.
Reconstruction of SDSL at the coast
The large-scale coherence between coastal and open-ocean steric sea level indicates that a purely data-driven reconstruction of SDSL signals along the coast might be possible. Here, we apply an EOF regression approach40 to compute the covariance relationship between coastal SDSLTG and steric height from each source region over 1960 to 2012. This relationship is then used to analyse SDSLTG solely on the basis of steric height in the open ocean (hereafter, SDSLEOF; Methods). The EOF effectively filters the common signal and removes noise from both datasets. The regression character further allows for a scaling of the steric height, which becomes necessary if signals are either amplified (for example, due to resonance processes generated by winds) or damped (for example, due to bottom friction) when travelling towards the coast. We note, however, that the EOF does not consider any time lags, which might play a role if baroclinic Rossby waves with long travel times are involved in the transfer.
To test the robustness of the reconstruction against overfitting and to assess uncertainties, we use ocean simulations from 12 Earth system models (ESM) over the period 1850–2012 with historical forcing (Supplementary Table 2). In the ESMs coastal SDSL, open-ocean steric height and ocean-bottom pressure variations are a priori known and therefore form an ideal testbed for the EOF approach (using source regions similarly determined in each ESM as for the observational record; Methods). We extract in total 12,000 realizations of 53-yr periods to calibrate the reconstruction (with randomly varied model parameters; Methods) and use an independent, random 53-yr period (from the remaining model years) for validation. The resulting time series are displayed, as an example, for the Northwest Atlantic south of Cape Hatteras in Fig. 2b and for all regions in Extended Data Fig. 4. In all model runs, the EOF reconstruction mimics coastal SDSL variations reasonably well, which is reflected in significant correlation coefficients over all models and subsamples providing only subtle differences between calibration and validation periods (Fig. 2c). EOF reconstructions in tropical and subtropical coastal regions with known links to the El Niño/Southern Oscillation show slightly larger correlations (median r ≥ 0.86, P ≤ 0.05) than those at higher latitudes (for example, the Northwest Atlantic north of Cape Hatteras) or in large continental shelf seas (for example, the North Sea) (median r = 0.64–0.77, P ≤ 0.05). This is probably related to increased and atmospherically forced barotropic variability in these records, which might mask the remotely forced and density induced SDSL variations28,41,42. Similarly well reconstructed are linear trends in the ESMs (Fig. 2d). These usually agree with the model internal true coastal SDSL signal within a range ±0.7 mm yr–1 in all subsamples and regions and show little evidence for systematic trend biases (median trend differences are of the order of ±0.06 mm yr–1 over both calibration and validation periods) (Fig. 2d). Overall, this demonstrates the skill of the EOF approach in reconstructing coastal SDSL from open-ocean steric height and motivates its application to observations.
The SDSLEOF reconstructions based on open-ocean steric height observations from the source regions are again displayed, as an example, for the Northwest Atlantic south of Cape Hatteras in Fig. 2b and for all regions in Extended Data Fig. 4. We generated an ensemble of 5,000 SDSLEOF reconstructions at each site, which is based on varying barystatic sea-level corrections at tide-gauge records (Methods) and randomly varied parameter choices in the EOF approach (Methods). In agreement with the ESM experiments, the ensemble SDSLEOF reconstructions show significant median correlations to SDSLTG that range from r = 0.62 (P ≤ 0.05) in the Northwest Atlantic south of Cape Hatteras to r = 0.95 (P ≤ 0.05) in the Western Pacific (Fig. 2c). In the three regions that are directly influenced by western boundary currents (Northwest Atlantic south of Cape Hatteras, Japan Sea and New Zealand), SDSLEOF displays slightly lower (but still good) agreement with coastal residual sea level than in the other regions. This is probably related to the more complex dynamics involving mesoscale eddy activity and slowly propagating Rossby waves39,43,44 that are hard to capture in the EOF approach without any consideration of time lags. Such time lags seem to play a minor role in regions dominated by coastally trapped waves, which typically have travel times below a month from the source regions to the coasts26,30. Overall, the EOF reconstructions display an encouraging performance in transferring open-ocean steric height signals toward the coast that exceeds those from standard reanalysis products (Extended Data Fig. 1) and even holds for longer timescales as further discussed in the budget assessment below.
Implications for the coastal sea-level budget
Taking advantage of the improved representation of SDSL variability at the coast, we finally reassess the total coastal sea-level budget for each region by adding the barystatic component back to the reconstructed SDSL component and comparing it to the vertical land motion corrected coastal sea level in each region. This budget is displayed in Fig. 3 for linear (Fig. 3a,b and Table 1) and nonlinear trends (Fig. 3c) during 1960 to 2012. In all regions the trend budget can be closed within the respective uncertainties (Fig. 3b and Table 1) that are, particularly in Northern Europe and North America, dominated by vertical land motion (Fig. 3b). Median differences to observations are everywhere within a range ±0.2 mm yr–1. The only exception is Hawaii, where the median trend is overestimated by 0.4 mm yr–1. We note, however, that there are also substantial uncertainties in the underlying steric products (Extended Data Fig. 5 and Methods) that may account for these differences. At all locations, except Hawaii and the Japan Sea, SDSL components explain the major fraction of the total trend budget (Fig. 3a).
Fig. 3: Coastal sea-level budget over 1950 to 2012. a, Shown are (as pie diagrams for each region) the fractions of the coastal trend budget that are explained by SDSLEOF and barystatic components. b, The linear trend budget of coastal sea level for nine regions after correcting each tide gauge for glacial isostatic adjustment and residual vertical land motion. Blue bars represent regional averages, while the stacked blue/orange bars represent the budget estimated with SDSLEOF and barystatic sea level. c, The rates derived from a nonlinear trend for each component and region. All shadings and error bars represent 95% CI derived from the 5,000-member ensemble.
Table 1 Linear trend budget. Linear trends for coastal sea level (corrected for vertical land motion), the total budget and each contributor over the period 1960 to 2012. All trends are given as a median of the 5,000-member ensemble with the 95% CIs provided in brackets. GRD, gravity, rotation and deformation. NW Atlantic NCH/SCH, Northwest Atlantic north/south of Cape Hatteras; NE Pacific, Northeast Pacific.
Region |
Observations (mm yr–1) |
Budget (mm yr–1) |
SDSLMEFOF (mm yr–1) |
Barystatic GRD (mm yr–1) |
|---|---|---|---|---|
North Sea |
2.01 (1.30;2.76) |
2.09 (1.58;2.52) |
1.73 (1.28;2.09) |
0.36 (0.17;0.56) |
NW Atlantic NCH |
1.43 (0.78;2.40) |
1.42 (1.02;1.87) |
0.76 (0.50;1.14) |
0.65 (0.39;0.87) |
NW Atlantic SCH |
1.40 (0.52;2.32) |
1.48 (1.12;1.80) |
0.76 (0.56;0.98) |
0.72 (0.42;0.97) |
NE Pacific |
1.23 (0.20;2.26) |
1.40 (0.97;1.80) |
1.02 (0.89;1.20) |
0.38 (–0.02;0.71) |
Hawaii |
1.55 (1.24;1.89) |
1.18 (0.71;1.58) |
0.47 (0.27;0.67) |
0.71 (0.31;1.05) |
Japan Sea |
1.29 (0.89;1.75) |
1.13 (0.64;1.60) |
0.48 (0.12;0.81) |
0.65 (0.35;0.96) |
West Australia |
1.48 (1.11;1.85) |
1.60 (1.10;2.04) |
0.89 (0.61;1.19) |
0.71 (0.34;1.03) |
New Zealand |
1.74 (1.16;2.35) |
1.66 (1.12;2.11) |
0.99 (0.64;1.25) |
0.66 (0.32;0.99) |
The budget terms also explain the observed interannual variability, which is expressed by an explained variance ranging from 38% at sites in the Northwest Atlantic south of Cape Hatteras to 91% along the Californian coast in the Northeast Pacific (Fig. 2c and Extended Data Fig. 4). The temporal variability is largely dominated by the SDSL component at all sites (Fig. 3c), although barystatic sea level through ice-mass loss and natural and anthropogenic terrestrial water storage variations (Methods) adds an accelerating signal to the rates. Largest accelerations in barystatic sea level are found around the Hawaiian Islands, where rates have been increasing from close to 0 to >2 mm yr–1 since the late 1990s. However, in the total budget, this acceleration has been reversed by strong negative anomalies in the SDSL components since the early 2000s (Fig. 3c). In contrast, Northern European coastlines only exhibit a small barystatic sea-level signal with an average rate of 0.4 (0.2; 0.6) mm yr–1 over 1960 to 2012 but this is (despite the large uncertainties in different steric products) counterbalanced by the largest average SDSL rates of all sites (Fig. 3b). In the Northeast Pacific, SDSL shows the well-known sea-level suppression since the 1970s45, while at gauges in the Western Pacific and West Australia, SDSL is characterized by a sustained acceleration from close-to-zero rates before the 1980s to >14 mm yr–1 in the 2010s (in an opposite direction to the negative SDSL rates at Hawaii); a value that is six times larger than the simultaneous barystatic sea-level rise. This acceleration has previously been attributed to anthropogenic forcing46 and represents the most pronounced SDSL sea-level change signal of all sites. Further accelerations in the SDSL terms can be seen in the Japan Sea and New Zealand since the early 1980s and along the coasts of Northeast America north of Cape Hatteras since the 1990s but their magnitudes are far smaller than those seen in the tropical Indo-Pacific regions and an attribution to natural and anthropogenic forcing has not yet been achieved.
The large site-specific rates also suggest substantial regional SDSL deviations from the simultaneous global mean (Fig. 4b) that are considerably larger than the spatial variations that have been contributed by barystatic sea level. Such deviations are important for coastal planning purposes, particularly at engineering-relevant timescales of a few decades. To generalize the potential of coastal SDSL variations to deviate from the global mean at varying timescales, we calculate moving trends for window sizes between 10 and 53 yr for both coastal SDSL and global mean steric sea level and assess their respective ratio (Fig. 4). At decadal timescales (~10 yr), local variations can be several hundred times larger or smaller than the global mean in any region considered here. This number decreases with increasing timescale but regional deviations can still be as large as ten, six and three times the simultaneous global mean for periods of 20, 30 and 50 yr, respectively. Given the dominance of SDSL variations at these timescales (Fig. 3c), this indicates that the key to more robust near-term coastal sea-level estimates for the coming decades clearly lies in a better understanding of SDSL variations and their geographical origin.
Fig. 4: Scaling of local and global SDSL. a, The maximum ratio between coastal SDSLEOF and global mean steric sea level observed for different periods and each of the nine regions between 1960 and 2012. The ratios have been derived from moving trends calculated for window sizes ranging from 10 to 53 yr. Three periods at 20, 30 and 50 yr have been highlighted, demonstrating that SDSL can differ by up to ten, six and three times from the global mean. The tick black line highlights the factor 1 indicating variations of magnitude equal to the global mean. The corresponding SDSL time series, smoothed with a singular spectrum analysis using an embedding dimension of 10 yr, are shown in b.
Our results represent a notable advance in estimating drivers of regional coastal sea level over the past 53 yr. Regardless of the accelerating barystatic sea-level terms around the globe, coastal sea-level variations are still largely dominated by the SDSL terms. This highlights the urgent need for a better understanding of the underlying physics to provide robust, near-term sea-level predictions. Our analysis represents a step in this direction as it demonstrates that for all analysed regions, coastal SDSL signals originate from the open ocean that are often thousands of kilometres away from individual sites. This reinforces the key role of large-scale ocean circulation in transferring the steric signal to the coast. The transfer can accurately be approximated using the EOF technique introduced here and allows for, together with the other components, improved closure of the coastal sea-level budget. Further modelling studies, such as those recently undertaken for the Northern European Shelf28,47, are required to clarify the involved oceanographic processes and to develop dynamical downscaling procedures that allow for more robust projections along the coast48.
Methods
Tide-gauge data and corrections
We make use of a set of 89 annual tide-gauge records from nine coastal regions taken from the online portal of the Permanent Service for Mean Sea Level in Liverpool49 and listed in Supplementary Table 1. The nine coastal regions have been selected on the basis of previous studies demonstrating that their tide-gauge records show similar dynamic sea-level variability4,21,25,29,36,39,43. In each region, individual records have been chosen on the basis of criteria such as data availability (>75%, except of a few exceptions primarily in the Indo-Pacific) and homogeneity (visual inspections and data flags). The individual records, together with their corresponding virtual stations and a cross-correlation matrix, are shown in Extended Data Fig. 2. Nine clusters of pronounced positive correlations appear, thus underpinning the coherence within each region. Existing data gaps have been filled with the local realizations from the hybrid sea-level reconstruction from ref. 50. The local realizations from the hybrid reconstruction include a local residual process from the Kalman Smoother that accounts for local effects such as vertical land motion (but does not contribute to the sea-level fields in the ocean and the corresponding global mean sea level) and therefore almost perfectly mimics their long-term trends50,51. They are also highly correlated with the real tide gauges (Supplementary Table 1). The gap-filling is a requirement for the EOF approach outlined below and avoids issues with benchmark differences when being merged to regional means52. We note, however, that the median percentage of data gaps is only 6% over all sites (Supplementary Table 1) and thus has a negligible influence on our results. Our major aim here is to investigate the SDSL variability in each of these regions. A challenge in this regard stems from the fact that tide-gauge records are affected by numerous different additional processes (vertical land motion and barystatic sea level) that may mask the actual SDSL signals. To isolate the SDSL variability at each site, we therefore initially remove vertical land motion and barystatic components. Each of these corrections comes with notable uncertainties. These uncertainties are considered here and build the basis for a probabilistic ensemble assessment with 5,000 members.
The first correction corresponds to vertical land motion, which primarily affects coastal sea level at lowest frequencies and induces spatially varying trends between individual locations. To correct for vertical land motion, we fit a linear trend to the differences between each tide-gauge record (before the gap-filling) and the nearest-neighbour time series from the hybrid reconstruction (that does not include the residual component and is corrected for its median glacial isostatic adjustment field) from ref. 50, which combines a process-based Kalman Smoother51 with ordinary EOF reconstructions53. As the Kalman Smoother aligns a series of predescribed processes (barystatic fingerprints, dynamic sea-level changes from climate models and 161 glacial isostatic adjustment models but excluding non-climatic local vertical land motion; ref. 51) to a global set of tide gauges, the difference between the hybrid reconstruction and tide gauges should, next to model errors, predominantly be driven by local vertical land motion51. The basic idea is similar to the Gaussian process approach used in ref. 54 for future projections and estimates based on the difference between tide-gauge records and satellite altimetry10 but it has the advantage of providing much longer residual series covering the entire length of tide-gauge records back to 1900 (that is, the period over which the hybrid reconstruction is available). As a preliminary cross-validation, we compare our estimates to the vertical land motion dataset from ref. 4. Their dataset is based on Global Navigation Satellite System observations and differences between tide gauges and nearby satellite altimetry but corrected for the nonlinear crustal components of present-day barystatic sea-level change and, therefore, is directly comparable to ours. The corrections are only available at 65 out of our 89 tide-gauge sites. Both datasets are significantly correlated (r = 0.78, P ≤ 0.05) with a root mean square difference of 0.8 mm yr–1 (which is substantially smaller than that from satellite altimetry minus tide gauge (1.22 mm yr–1); ref. 55) but vertical land motion estimates from the hybrid reconstruction show a higher correlation to linear trends from tide-gauge observations (r = 0.91, P ≤ 0.05) than those from ref. 4 (r = 0.67, P ≤ 0.05) (Extended Data Fig. 6). For this reason, as well as the fact that the hybrid reconstruction provides vertical land motion estimates at more stations than the ref. 4 dataset, we proceed with the hybrid reconstruction-based vertical land motion estimates. The largest uncertainty in the hybrid reconstruction-based vertical land motion estimates stems from the 161 glacial isostatic adjustment models used in the Kalman Smoother from ref. 53 that consider different solid-Earth parameters (lithosphere thickness and mantle viscosity) and varying global deglaciation histories over the past 20,000 yr. The second uncertainty in the vertical land motion estimates is related to fitting of the linear trend to noisy data. This uncertainty has been modelled considering that the residuals are normally distributed but temporally correlated. To consider both uncertainties in the budget assessment, we perturb the vertical land motion estimates with random noise from the fitting uncertainty as well as the uncertainty from the 161 glacial isostatic adjustment models (considering their spatial correlation). This results in a 5,000-member ensemble of plausible vertical land motion corrections at each site.
The next components that we remove from the tide-gauge records are the barystatic terms due to contemporary mass redistribution. Here, we make use of a 5,000-member ensemble by ref. 4 that combines sea-level contributions from ice-sheet56,57,58,59,60, glacier11,61 and terrestrial water storage62,63,64,65 observations and reanalysis estimates and accounts for observational and model uncertainties. The barystatic effects have little effect on interannual variability of sea level but they are highly nonlinear and produce large spatial variability between individual locations (Fig. 3b). Removing the two 5,000-member ensembles of vertical land motion and barystatic sea level results in an ensemble of SDSLTG residuals including the uncertainties resulting from the initial corrections.
In addition to vertical land motion and barystatic effects, tide-gauge records are also affected by the barotropic response of the ocean to atmospheric forcing. This consists of barotropic wind forcing and the inverted barometer response to sea-level pressure fluctuations over the oceans and represents a dominant fraction of the sea-level spectrum from intra-annual to decadal scales42. While the barotropic wind forcing term is usually considered as an inherent process of SDSL, it may introduce large variability of opposite sign at different locations27,28 and mask the actual SDSL variability seen by a tide gauge41. The latter may particularly be the case, when SDSL signals have, for instance, been initiated by wind forcing in remote regions affected by different atmospheric circulation systems than those responsible for the local barotropic signals. For this reason, we initially remove these wind and pressure effects from each tide-gauge record for the correlation and EOF analyses. Note, however, that we later add this component back to SDSL for the overall validation and the budget assessment. Thus, the SDSLEOF components shown in Figs. 1–3 in fact represent SDSL plus atmospheric pressure effects. As an estimate of these effects, we use the outputs from barotropic simulations with the MIT global circulation model forced with atmospheric reanalysis winds and sea-level pressure from the twentieth century reanalysis66 over the period 1871–2012. Further details on the model configuration and detailed validations can be found in refs. 21,42. The wind and pressure components are most important at higher latitudes, where stratification is weak, and they introduce both spatially varying temporal trends as well as pronounced intra-annual to decadal variability [42].
Steric height data
We estimate steric changes from a gridded temperature24 and salinity reconstruction67 on the basis of hydrographic profiles covering the upper 2,000 m of the ocean over the 1960 to 2012 period using the TEOS-10 GSW software for MATLAB68. This temperature and salinity reconstruction is based on an optimal interpolation approach, which uses the CMIP5 multimodel ensemble to derive correlation scales and to set the initial field. We note that there are many other gridded temperature and salinity datasets available that use different mapping methods than ref. 24. Three of these datasets have additionally been analysed here: ref. 69 and EN4 (EN4.2.1.)70 with two different mechanical bathythermograph (MBT) and expendable bathythermograph (XBT) bias corrections schemes following refs. 71,72. There are notable differences in linear trends calculated from each of the steric height datasets over 1960 to 2012 (Extended Data Fig. 5a–d) that are particularly pronounced in areas of major ocean circulation systems (>2 mm yr–1)73. These large differences between individual products may be partly related to the distinct mesoscale eddy activity that is not resolved by hydrographic observations as well as the different mapping approaches in each reconstruction24. They also feed in the SDSLEOF reconstruction (that essentially covers information from these circulation systems) when applied to each dataset individually (Extended Data Fig. 5e,f). However, the SDSLEOF reconstruction based on the ref. 24 dataset has the smallest trend differences to the SDSLTG in most regions (Extended Data Fig. 5) and also provides the best representation of variability (Extended Data Fig. 5f). It is interesting to note that regions with the largest interproduct trend spread also show the largest spread in correlations to SDSLTG with ref. 24 data providing far better agreement than all other products. For this reason, and since the dataset from ref. 24 has been shown to represent a methodological improvement over former reconstructions, we limit our analysis in the main paper to this particular dataset.
EOF reconstruction
The EOF approach assumes that two or more sets of variables share a certain degree of similar variability, which can be expressed by distinct spatial patterns (EOF modes) for each individual variable (here, steric height from the open ocean and SDSLTG at the coast) and a common principal component (PC) corresponding to each spatial pattern. In the literature, similar approaches (but using canonical correlation analysis74 or cyclostationary EOFs40) have been used to reconstruct sparse/short climatic data (for example, precipitation) on the basis of covariates covering much longer periods than the variable of interest (for example, sea-level pressure or sea-surface temperature). Here, we use this approach to transfer steric height estimates from the source regions (identified with the correlation patterns in Fig. 2a) to the coast represented by the tide-gauge residuals after removing vertical land motion, barystatic sea level and (barotropic) wind and pressure effects (SDSLTG). This is done in three steps:
We calculate EOFs over the 1960–2012 period between detrended and smoothed (using a randomly chosen smoothing window between 3 and 5 yr) steric height and SDSLTG. This produces several different modes with variable-specific spatial patterns and common PCs. In our reconstruction, we only use a subsample of these modes that share a certain degree of variance in each region. We randomly varied the corresponding threshold, such that only those modes are considered that cumulatively explain between 92 and 98% of the variance of the entire coupled field. On average, this resulted in the consideration of five (West Pacific) to nine (North Sea) modes in the different regions.
The spatial pattern of the steric height from the EOF analysis is then projected back onto the entire (1960–2012) non-detrended and non-smoothed steric height fields from the source region to produce a series of each PC.
Finally, the PCs from step 2, now solely based on steric height from the open ocean, are projected onto the spatial patterns from tide-gauge residuals (SDSLTG), which produces an EOF-based SDSL reconstruction (SDSLEOF) over the entire 1960–2012 period at each site.
The SDSLEOF reconstructions are then evaluated as spatial averages over all sites in each region. By using varying barystatic corrections and by considering randomly varied parameter choices, we produce 5,000 ensemble members that are used to assess uncertainties of the EOF reconstruction. The EOF approach has two major advantages compared to simple field averages over the source region as formerly used in refs. 22,23 While a simple area average assumes that the steric signal from the open ocean is one on one transferred toward the coast, there are factors at play that may dampen or amplify the signal while propagating from the open ocean toward the coast. The regression character of the EOF approach explicitly takes this into account and scales the signal to match the variability seen at tide gauges. Furthermore, EOF analysis decomposes a field into common modes. As typical for many other atmospheric or oceanographic time series, most of the variance of the entire field is contained in the first few modes, while local effects move as noise into lower EOFs. Thus, the EOF approach is robust against local outliers. To test the stability of the EOF reconstruction against increasing data uncertainties in the steric products before the 1980s, we also calculated a second 5,000-member ensemble on the basis of PCs calculated since solely 1980 (Extended Data Fig. 7). We do not find any significant differences in the reconstructions based on PCs calculated over the entire period or only since 1980. Correlations are, in all cases, >0.98 with the only exception being the three western boundary currents where the correlation coefficients are between 0.92 and 0.94. The trends are also not significantly different between the two ensembles (Extended Data Fig. 7b).
Earth system models
To test the performance of the EOF technique, we apply it to the historical runs of 11 ESMs from the Coupled Model Intercomparison Project Phase 5 (CMIP5) [75], which are all listed in Supplementary Table 2. As CMIP5 models usually have a horizontal resolution that barely exceeds one degree, we additionally included one eddy-permitting model (CNRM-CM6-1-HR; ref. 76) from the HighResMIP initiative. CNRM-CM6-1-HR has a resolution of 25 km, providing a more realistic representation of the coastal zone than its CMIP5 counterparts and many ocean reanalyses19. We decided to use ESMs rather than ocean reanalyses because the latter suffer from (regionally varying) drifts due to sparse data assimilation, uncertainties in air–sea fluxes and mixing errors19, while ESMs are entirely consistent without data assimilation. All simulations have been linearly de-drifted using the pre-industrial control run (PiControl). ESMs have the advantage of providing predescribed components of the sea-surface height (‘zos’), ocean-bottom pressure (‘pbl’) and global mean thermosteric sea-level changes (‘zostoga’). As the ESMs are volume- rather than mass-conserving, we sum zos and zostoga and subsequently subtract pbl to derive steric height from each model (one might also estimate steric height from temperature and salinity profiles in each model, which is, however, more time consuming). We sample coastal sea-level (zos + zostoga) at the grid points closest to the tide-gauge locations and perform a similar correlation analysis with steric height to identify the source regions. Then the EOF technique is applied in a similar fashion (with 1,000 randomly varied smoothing factors and mode-selection thresholds) as for observations. We compare the reconstruction to both the available 53-yr calibration periods as well as a randomly selected and independent validation period of 53 consecutive years that have not been used for calibration. There are three factors to keep in mind when using ESMs as a testbed for the EOF approach. First, the resolution of ESMs is, even in CNRM-CM6-1-HR, still too coarse to fully resolve coastal processes. Second, the model provides ‘perfect observations’ at each location without gaps and measurement errors. Third, we do not have corrections for barotropic wind forcing available for the ESMs as we apply for real observations. While the first two factors spatially smooth the observations and therefore tend to increase the covariance compared to real observations, the missing barotropic correction may lead to a slightly degraded performance in high latitude regions compared to reality.
Statistics
All correlations and regressions were performed after removing linear trends. Statistical significance and error estimates for the correlation analyses are computed using Fisher’s z-transform77. While the test does not account for temporal correlations in the data, it is much more time-efficient and leads to barely different source regions than with repeated simulations using autoregressive processes. Uncertainties of linear trends are calculated assuming that the residuals are temporally correlated following an autoregressive process of the order one. Nonlinear trends in the different budget components have been calculated by a singular spectrum analysis with an embedding dimension of ten using the MATLAB package from ref. 78.
Extended data
Extended Data Fig. 1: Performance of ocean reanalysis in simulating SDSL variability and trends. Shown are linear trends of SDSL as simulated by the a SODAsi.378, b SODA 2.2.480, c ORA S481, d ORA20C82,83, and e GECCO214 reanalysis over the common period from 1960 to 2012. For all reanalysis systems, the model internal global average has been replaced by ref. 24. In b and d, data is only available until 2008 and 2009, respectively. Grey dots show the 89 tide-gauge locations used in this study. f Linear trends in SDSLTG from the virtual stations of the nine coastal regions (grey bars) are compared to trends calculated from SDSL as simulated by the five ocean reanalysis systems (nearest-neighbour series). d, Correlations between detrended SDSLTG and SDSL from the five ocean reanalysis systems. The grey shadings separate the different regions from each other.
Extended Data Fig. 2: Tide-gauge coherence and virtual stations for each region. a, Cross-Correlation matrix for the 89 tide-gauge records ordered by region. Black boxes mark the locations of the selected tide-gauge records for each region. b, The observed tide-gauge records (corrected for vertical land motion; coloured lines) together with the virtual station for each region (thick black line) that has been built based on gap-filled records (see Methods). The percentage of total data availability in each region is given in brackets.
Extended Data Fig. 3: Origin of coastal SDSL changes. Spatial correlation patterns between the SDSL residuals at tide gauges (SDSLTG, grey dots) and steric height24 from the open ocean used to assemble Fig. 1a but for each of the nine coastal regions separately. Only significant correlations (P ≤ 0.05) are shown. a, North Sea, b NW Atlantic north of Cape Hatteras, c NW Atlantic south of Cape Hatteras, d NE Pacific, e Hawaii, f Japan Sea, g West Australia, h New Zealand, and i West Pacific.
Extended Data Fig. 4: Reconstruction of coastal SDSL changes. Extension of Fig. 2b illustrating the reconstruction of SDSLEOF in observations (top) and in validation experiments with ESMs (here illustrated by the CNRM-CM6-1-HR, bottom) for each of the nine coastal regions. a, North Sea, b Northwest Atlantic north of Cape Hatteras, c Northwest Atlantic south of Cape Hatteras, d Northeast Pacific, e Hawaii, f Japan Sea, g West Australia, h New Zealand, and i West Pacific.
Extended Data Fig. 5: Linear trends in steric height and comparison of different observational products. Shown are the linear trends in steric height calculated over the upper 2000m for different gridded observational products. a, ref. 24 and 67, b ref. 68, c EN469 with ref. 70 corrections, and d EN4 with ref. 71 corrections. The grey dots mark the locations of tide-gauge records used in this study. e, Linear trends for the SDSLEOF reconstructions in each region using the four different data products (grey bars = ref. 24 and67; blue = ref. 68; turquoise = EN469 with ref. 70 corrections; yellow = EN469 with ref. 71 corrections. f, Same as e but showing the correlation between SDSLEOF based on the different products and SDSLTG.
Extended Data Fig. 6: Validation of the vertical land motion (VLM) correction. Comparison between observed trends (after the removal of barystatic gravitation, rotation and deformation terms) and residual VLM plus Glacial Isostatic Adjustment from the difference between tide gauges and the hybrid reconstruction from ref. 50 as well as the observed trends and residual VLM from Global Navigation Satellite System plus Glacial Isostatic Adjustment and the difference between tide-gauge and satellite altimetry as calculated by ref. 4 (Fred).
Extended Data Fig. 7: Validation of the EOF approach. a, Shown are the time series of SDSLEOF based on reconstructions using principal components that have been calculated and regressed on the steric height from the open-ocean over the entire 1960 to 2012 period (and as used in the main paper) as well as those based on principal components that have been calculated (and regressed on the steric height from the open-ocean) over the period from 1980–2012. The grey shading marks the corresponding validation period from 1960–1979. b, The corresponding linear trends of SDSLEOF over the common period from 1960–2012. Shadings and error bars represent the 95% confidence intervals.