Timing of emergence of modern rates of sea-level rise by 1863

Sea-level rise is a significant indicator of broader climate changes, and the time of emergence concept can be used to identify when modern rates of sea-level rise emerged above background variability. Yet a range of estimates of the timing persists both globally and regionally. Here, we use a global database of proxy sea-level records of the Common Era (0–2000 CE) and show that globally, it is very likely that rates of sea-level rise emerged above pre-industrial rates by 1863 CE (P = 0.9; range of 1825 [P = 0.66] to 1873 CE [P = 0.95]), which is similar in timing to evidence for early ocean warming and glacier melt. The time of emergence in the North Atlantic reveals a distinct spatial pattern, appearing earliest in the mid-Atlantic region (1872–1894 CE) and later in Canada and Europe (1930–1964 CE). Regional and local sea-level changes occurring over different time periods drive the spatial pattern in emergence, suggesting regional processes underlie centennial-timescale sea-level variability over the Common Era. Sea-level rise is a significant indicator of climate changes and it is important to identify the time of emergence of modern rates of sea-level rise. Here the authors estimate that global sea-level rise emerged by 1863 and find spatial variability of emergence at sites within the North Atlantic.

T he time of emergence (ToE) identifies when a climate change signal emerges above background variability, reflecting the onset of significant periods of change [1][2][3] . ToE has been evaluated for measures such as surface air temperatures, atmospheric variables, sea surface temperatures (SST), ocean biogeochemistry, and sea level 1,2,4 . However, these studies are often restricted to comparisons of 21st-century projections to a background reference period of observations over decades to the last century or simulated pre-industrial climate rather than using available paleoclimate data 5 . With regard to sea level, highresolution proxy reconstructions provide a record of preindustrial variability extending through the Common Era (0-2000 CE) 6,7 , which allows for a more detailed analysis of ToE.
The goal of identifying the timing of modern rates of sea-level rise has been pursued using several different techniques and datasets. Although there is agreement that rates of sea-level rise globally and in many locations exceed Common Era background rates by the late 19th to early 20th century [8][9][10] , global and regional estimates differ and the spatial variability in the timing among locations is unclear. Using a global tide-gauge compilation beginning in 1870 CE, Church and White 8,11 found that rates of global-mean sea-level rise were already accelerating in the late 19th century. A global sea-level reconstruction using Monte Carlo Singular Spectrum Analysis of tide gauge records suggests an acceleration starting at the end of the 18th century 12 , but only a limited number of tide gauge records, restricted to northwestern Europe, extend back through to the 18th and early 19th centuries. Kopp et al. 9 estimated global sea-level change through the Common Era by applying a spatiotemporal hierarchical model to a global database of relative sea-level (RSL) reconstructions and found a significant global sea-level acceleration that began in the 19th century. From the combination of tide gauge observations and proxy reconstructions, the Sixth Assessment Report (AR6) of the Intergovernmental Panel on Climate Change (IPCC) concluded that a sustained increase of global mean sea level began between 1820 and 1860 3 . At a regional scale, Gehrels and Woodworth 10 used tide gauge and proxy reconstructions from the North Atlantic and Australasia to suggest an increase in the rate of rise between 1895 and 1945 CE; however, this estimate was simply made by visual inspection of sea-level trends. Previous studies on the U.S. Atlantic coast have used change point analysis to examine the timing of increased rates of rise at individual locations, quantifying common timing among western North Atlantic proxy records to 1865-1873 CE 13 . At individual sites, change point analysis and linear regression have identified increases in rate varying in timing from the early 19th to the early 20th centuries 7,14-17 . However, the change-point model is limited in that it assumes acceleration is instantaneous.
The range of suggested timings both globally and regionally may reflect the method of analysis or the influence of the longterm rate of change and the amplitude of pre-industrial variability among locations or may arise from regional variability in sea-level change, due to processes such as the gravitational, rotational, and deformational fingerprints of mass loss from ice sheet and glacier melt or changes in ocean and atmosphere circulation 3,18 . Concurrent analysis of sea-level change, both globally and at many individual sites, is needed to identify the source of discrepancies among estimates of the timing of increased rates of sea-level rise and to determine when sea-level rates emerged to be clearly distinguishable from background variability.
Here, we use a previously published global database (Supplementary Fig. 1, Source Data) of instrumental and proxy (e.g., foraminifera, diatoms, testate amoebae, coral microatolls, archeological evidence, sediment geochemistry) sea-level records of the Common Era 19 to examine the ToE of modern rates of sea-level rise above pre-industrial background variability. These records are incorporated into a spatiotemporal empirical hierarchical model 9,20 to examine magnitudes and rates of past RSL and global sea-level change with associated uncertainty, using shorter time periods (60-year rates) than the previous analysis 9, 19 while still minimizing the effects of interdecadal fluctuations (Supplementary Table 4). We analyze the timing of the onset of modern rates of sea-level rise using the entire reconstructed spatiotemporal sea-level field from the global database. The ToE of modern rates of RSL rise is defined as when it is very likely (P ≥ 0.90) that the rate of RSL change during a 60-year industrial-era period (from 1700 to 2000 CE) and all subsequent periods is greater than that of a random pre-industrial period (from 0 to 1700 CE). We improve upon previous analyses by assessing both global ToE, using the common global sea-level signal from proxy records, and local ToE at individual sites probabilistically and simultaneously, as opposed to separate studies of individual sites or an analysis of only the global signal. In addition, we evaluate the spatial variability of the ToE among individual locations in the North Atlantic, where the highest density of reconstructions is located, and we only analyze those sites with the highest resolution reconstructions. By decomposing RSL change into different spatiotemporal scales, we can identify regional anomalies in RSL, suggesting potential underlying processes contributing to the variability in ToE.

Results
Global ToE. The rate of global sea-level rise emerged above preindustrial rates by 1863 CE (P = 0.9; range of 1825 [P = 0.66] to 1873 CE [P = 0.95]), when the global rate of rise was 0.4 ± 0.2 mm/yr from 1840 to 1900 CE (Fig. 1). Over the pre- Global sea-level rise is largely driven by thermal expansion of warming ocean water and increases in ocean mass due to the melting of land-based glaciers and ice sheets 22,23 . The global ToE of modern sea-level rise (1863 [1825-1873] CE) is in concordance with the onset of warming oceans. While a global synthesis of sea surface temperatures identified a cooling trend from 1 to 1800 CE 24 , instrumental records of sea surface and subsurface temperatures suggest warming from the 1870s to present 25 . Using paleoclimate records since 1500 CE, Abram et al. 5 found that sustained, significant industrial-era warming trends commenced during the 19th century, preceding the global sea level ToE. Specifically, an onset of warming of sea surface temperatures occurred in the early to mid-19th century in the Western Atlantic (1828 CE), Western Pacific (1834 CE), and Indian oceans (1827 CE). While intermediate water temperatures in some regions of the ocean interior or global ocean heat content may lag surface warming trends 26,27 , the onset of increased surface temperatures may signify the initial change in ocean warming, which precedes the global sea level ToE. The global ToE of modern sea-level rise is also similar in timing to global glacier retreat. Glaciers expanded globally beginning in the 13th century and reached a maximum between the mid-16th to 19th centuries, but subsequently began to retreat in the middle to late 19th century [28][29][30] . It is unlikely that the Greenland and Antarctic Ice Sheets had large positive contributions to global sea-level rise prior to the 20th century 31,32 . The Greenland Ice Sheet had minimal variability in ice-mass loss over the Common Era 32 but has experienced an acceleration in the rate of ice loss in the last century 33 . The mass balance of the Antarctic Ice Sheet is not as well constrained until the last several decades when there has been an acceleration in mass loss over time 34 . Additionally, the ToE occurs after the end of the Little Ice Age (~1850 CE), which was marked by the recovery from naturally forced climate cooling due to a series of volcanic eruptions and solar variability 35,36 .
ToE in the North Atlantic. Local records in the North Atlantic exhibit more variability than the common global signal; therefore, the ToE of modern rates of global sea level (1863 [1825-1873] CE) is earlier than individual locations in the North Atlantic. At sites in the North Atlantic, the ToE exhibits spatial variability and ranges from the late 19th century to the mid-20th century ( Table 1, Fig. 2). The spatial pattern in the ToE is unrelated to pre-industrial rate variability or proxy data resolution (Supplementary Figs To highlight the underlying nonlinear regional to local signal at each site responsible for the spatial pattern of the ToE, we use the spatiotemporal model decomposition (Fig. 3) to remove the global signal (which is common to all sites) (Fig. 3b) and linear signal (which is consistent with the effects of glacial isostatic  (Fig. 3c). There is a clear spatial variability of regional (global and linear signals removed) trends that correlates to the spatial pattern in the ToE (Supplementary Fig. 5). In the mid-Atlantic, it is very likely (P > 0.9) that the regional rates are positive for the period from 1700 to 2000 CE, with an increase beginning around 1400 CE, reaching up to 0.7 ± 0.4 mm/yr at the end of the 20th century in southern New Jersey ( Supplementary  Fig. 5a, Supplementary Table 3). The regional increase in rate is muted or absent in the northeastern and southeastern U.S. (e.g., 0.2 ± 0.4 mm/yr in Massachusetts and 0 ± 0.4 mm/yr in northern Florida at the end of the 20th century), where the probability of a positive regional rate from 1700 to 2000 CE ranges from 0.7 to 0.9 ( Supplementary Fig. 5b, Supplementary Table 3). In Canada and Europe, there is no continuous pre-industrial regional rate increase to the present. In fact, it is likely (P > 0.66) that most sites in Canada and Europe had a negative regional rate from 1700 to 2000 CE (e.g., −0.6 ± 0.4 mm/yr in Denmark and Scotland by the end of the 20th century), with a decline beginning around 1800 CE ( Supplementary Fig. 5c, Supplementary Table 3). The spatial variability in the ToE among sites suggests different centennial-timescale, regional spatial scale physical mechanisms have influenced North Atlantic sea level over the last millennium. Site-specific local mechanisms, such as increasing groundwater withdrawal causing land subsidence from growing populations at mid-Atlantic sites, could contribute to faster rates of RSL rise. For example, coastal New Jersey experienced up to~0.7 mm yr −1 of subsidence in the 20th century 37 and other sites on the U.S. Atlantic coast had subsidence rates in the 21st century up to double the long-term geologic rates due to groundwater withdrawal 38 . However, these are relatively recent processes and the consistency of the pattern of the ToE among and within regions suggests common regional scale driving mechanisms rather than local site-specific processes.
Along the U.S. Atlantic coast, the spatial pattern of regional rate anomalies does not suggest mass loss of the Greenland Ice Sheet, where we would expect to see an amplified RSL response to the south along the U.S. Atlantic coast 39 , and therefore also a south to north pattern of the earliest to latest ToE. Previous analysis has also not detected a fingerprint of Greenland Ice Sheet melt on the U.S. Atlantic coast over the Common Era 40 . An ice sheet fingerprint may be too small or overprinted by other processes to be detectable over this time period; however, it is unlikely that the Greenland Ice Sheet significantly contributed to sea-level rise prior to the 20th century 32,33 .
We thus hypothesize that the regional increase in RSL rates on the U.S. Atlantic coast is related to regional changes in the ocean and/or atmosphere. Ocean circulation changes in the North Atlantic through Atlantic Meridional Overturning Circulation (AMOC) and the strength and/or position of the Gulf Stream can affect regional sea level on the U.S. Atlantic coast [41][42][43] . For example, many theoretical and modeling studies have found a scaling coefficient between AMOC transport and U.S. Atlantic coast sea level on the order of −1 to −2 cm/Sv 44 . Specific climate changes coincident with the initiation of elevated regional rates of RSL rise in the mid-Atlantic have been documented in the proxy record. In particular, proxy foraminifera have been used to infer a~3 Sv reduction in the Gulf Stream through the strength of the Florida Current during the Little Ice Age (at~1350-1750 CE) 45 . However, sea-level rise associated with a Florida Current weakening is expected to be largest and most coherent south of Cape Hatteras; 44 therefore, the coupling of changes across the North Atlantic involving the Florida Current, Gulf Stream, and U.S. Atlantic coast sea level requires further analysis. The timing of the circulation changes and the regional increase in RSL rates in the mid-Atlantic is also roughly coincident with a shift from a sustained positive phase of the North Atlantic Oscillation (NAO) to a sustained negative phase beginning around 1400 CE 46 . Changing patterns of atmospheric winds, buoyancy fluxes, and pressure, such as those connected with the NAO 47 , can affect the Gulf Stream and AMOC 48,49 , which could manifest in centimeter-scale regional sea-level changes 50 .
At sites in Canada and Europe, regional rates are relatively stable until~1800 CE; however, there are several short periods of increased pre-industrial rates at several sites. Specifically, there is an increase in the regional rate between 1300 and 1500 CE, reaching 0.7-0.8 mm/yr in both sites in Greenland and one of the two sites in Newfoundland, but the rates then subsequently decrease into the 20th century ( Supplementary Fig. 5). The increase in the rate in Greenland may be related to a period of growth of the Greenland Ice Sheet and crustal subsidence during this time 51 . The site in Newfoundland may experience more varied RSL trends due to its local-scale geomorphology situated on a peninsula indirectly connected to the ocean 40 . These elevated pre-industrial rates could contribute to the later ToE at these sites; however, elevated pre-industrial rates are absent at other sites within the same region that have a similar ToE, suggesting that the early elevated rates are locally driven and the ToE is substantially driven by regional processes. We thus conclude that the later ToE in Canadian and European sites is driven by a regional rate decline beginning around 1800 CE, which masks the positive rate contribution from the increasing global signal occurring during the same time. Whether the decline is related to climatic changes in the North Atlantic such as a proposed coincident AMOC weakening 52,53 is uncertain. Model evidence and theory suggest that sea-level changes associated with AMOC are smaller along the eastern boundary of the Atlantic 44,54 . There is a potential explanation outside ocean dynamics, at least for the last century: a regional sea-level fall in eastern Canada and northwestern Europe is consistent with the sea-level fingerprint of Greenland ice mass loss 55 . Overall, the ToE of modern rates of sea-level rise is later along the eastern North Atlantic margin (European coast) compared to the western North Atlantic margin (North American coast), supporting previous interpretations of differences in RSL histories between the North American and European coastlines over the last several centuries 56,57 . Therefore, while modern rates of RSL emerge at all sites by the mid-20th century due to the large magnitude common global signal (1.4 ± 0.2 mm/yr by the end of the 20th century), the spatial variability in the ToE among sites is driven by asynchronous regional centennial-timescale trends unrelated to long-term linear rates of change associated with glacial isostatic adjustment (Fig. 3,  Supplementary Fig. 5). Specifically, regional and local sea-level changes with an increase on the U.S. Atlantic coast since~1400 CE and a decrease in Canada and Europe since~1800 CE drive the spatial pattern in ToE, suggesting multiple processes underlie North Atlantic RSL variability over the Common Era. Further analysis of climate proxies and the relationship of elements of the North Atlantic circulation with RSL is needed to fully determine the cause of the regional RSL changes occurring over the last millennium, but the timing of these trends suggests that the spatial patterns in ToE are unrelated to anthropogenic climate forcing.

Methods
Sea-level database. The global sea-level database 19 (Source Data) comprises Common Era RSL proxy records with high-resolution chronologies from 36 regions around the world (Supplementary Fig. 1) and includes 2274 sea-level data points from proxies such as foraminifera, diatoms, testate amoebae, coral microatolls, archeological evidence, and sediment geochemistry. In addition, decadalaverage values from instrumental tide gauge records in the Permanent Service for Mean Sea Level (PSMSL; 58 were included, provided they were either (1) longer than 150 years, (2) within 5°distance of a proxy site and longer than 70 years, or (3) the nearest tide gauge to a proxy site that is longer than 20 years 9,40 . Also included are multicentury records from Amsterdam (1700-1925 CE) 59 , Kronstadt (1773-1993 CE) 60 , and Stockholm (1774-2000 CE) 61 , as compiled by PSMSL. As in Kopp et al. 9 and Kemp et al. 40 , the input data also include 1880-2010 global mean sea-level reconstruction of Hay et al. 62 from tide-gauge records.
Spatiotemporal model. The data is used with a spatiotemporal empirical hierarchical model 9,20 . A process-level characterizes RSL over space and time and a data level links RSL observations (reconstructions) to the RSL process.
At the process level, the RSL field f x; t ð Þ is modeled as the sum of seven components 40 where x represent geographic location, t represents time and t 0 is a reference time point (2000 CE). The components are fast and slow common global terms (g f t ð Þ and g s t ð Þ), a regional linear term (m x ð Þ t À t 0 À Á ), fast and slow regional non-linear terms (r f x; t ð Þ and r s x; t ð Þ), and fast and slow local terms (l f ðx; tÞ and l s x; t ð Þ). Predictions from the ICE5G-VM2-90 Earth-ice model 63 act as prior means for the regional linear term.
The data level includes the RSL reconstructions with observations, y i , where where f x i ; t i À Á is the true RSL value at location x and time t, ε i is the vertical uncertainty, w x i ; t i À Á is supplemental white noise, and y 0 ðx i Þ is a site-specific vertical datum correction to ensure that the RSL reconstructions are directly comparable to one another. The true age of an RSL observation (t i ) is the mean estimate (t i ) and its error (δ i ).
The hyperparameters are set through a maximum-likelihood optimization to characterize prior expectations of spatial and temporal scales, as well as amplitudes, of RSL variability (Supplementary Table 1). The non-linear terms were characterized by three spatial scales (global, regional, and local) and two temporal scales (fast and slow), which allow RSL to be decomposed into global, regional temporally linear, regional non-linear, and local components. Here, we remove a constraint on the model that was applied in Kopp et al. 9 , Kemp et al. 40 , and Walker et al. 19 that mean global sea level over −100 to 100 CE is equal to mean global sea level over 1600-1800 CE because, with the current, expanded database, the results of the analysis with and without the constraint are nearly the same; a constant± 0.1 mm/yr ambiguity in the long-term trend is no longer apparent.
Time of emergence. We use the reconstructed spatiotemporal field to determine the timing at which the rates in the last three centuries of the records emerge above the spread of the previous variability over the Common Era (Fig. 1). To minimize the effects of interdecadal fluctuations 12,64 and limited reconstruction resolutions, we focus on 60-year average rates. We define the pre-industrial background rates by the distribution of 60-year averages from 0 to 1700 CE at 20-year increments (e.g., 0-60 and 20-80 CE). Industrial-era rates are defined by 60-year averages from 1700 to 2000 CE at 20-year increments (e.g., 1700-1760 and 1720-1780 CE). The pre-industrial 60-year periods are enumerated from 1 to m, while industrial-era periods are enumerated from m + 1 to m + n. The ToE of modern rates of RSL rise is defined as when it is very likely (P ≥ 0.90) that the rate of RSL change during industrial-era period l and all subsequent periods is greater than that of a random pre-industrial period. We define rates as x = [x 1 , x 2 , …, x m , …, x m+n ], where x follows a multivariate normal distribution estimated by the statistical model. We define U as the (m + n) × (m + n) matrix with elements u ij = x i − x j . Accounting for the full covariance among rates, we compare the industrial-era rates to the background distribution using a Monte Carlo approach, taking 10,000 samples of U, enumerated U 1 ; ;Û 10000 , with elements enumeratedû ijk . For each Monte Carlo sample k, we randomly select one pre-industrial period r k 2 ½1; m to serve as the reference period. For each industrial-era period l, we calculatev lk ¼ min i2½l;mþnûir k k and estimate the distribution of v l from the samplesv lk . The ToE of modern rates of RSL rise is defined as the period l when P(v l > 0) ≥ 0.90. In addition, representing each 60-year period by its central year, we interpolate the probability curve to identify the year in which the probability reaches 0.90. We provide an uncertainty range for the ToE with a lower bound of when the probability reaches 0.66 and an upper bound of when the probability reaches 0.95.
Both the highest density of RSL reconstructions from the Common Era database and the highest resolution reconstructions are in the North Atlantic. We, therefore, focus our ToE analysis on individual sites in the North Atlantic where the spatial resolution allows an examination of variability in timing. To increase the likelihood that any variability in the ToE is due to process and not the proxy data resolution, we also only analyze the ToE at proxy reconstructions in the North Atlantic that fit the following criteria: (1) the proxy record is at least~1000 years in length to provide sufficient background information; and (2) the proxy record has data from 1700 to 2000 CE. Using these criteria, we analyze the ToE at six sites in Europe and fifteen sites along the eastern coast of North America (Table 1). In addition, we establish a "null hypothesis" by predicting RSL and the ToE at an indicative generic site in northeast Asia (38°N, 127°E) with no instrumental data and far from available proxy data to compare with the RSL and ToE at individual locations where there is available instrumental and proxy data ( Supplementary  Fig. 1). At this site, predicted RSL is equal to global sea level plus additional uncertainty associated with regional variability (Supplementary Fig. 2). When the RSL predictions and ToE at individual records differ from the null hypothesis, it reflects the influence of meaningful local information at individual sites.
We also conducted a series of sensitivity analyses to evaluate our ToE method. When comparing the ToE when using 40, 60, or 80-year average rates, the global ToE only varies by 15 years between the ToE using 40-year rates and 80-year rates (Supplementary Table 4). Using 40-year rates does shift the ToE at individual sites to a later date as the effects of interdecadal fluctuations become more pronounced, but the overall spatial pattern in the ToE remains broadly the same with earlier ToE in the western North Atlantic compared to the eastern North Atlantic. Using 80-year rates shifts the ToE at individual sites to slightly earlier dates, but again, the spatial pattern is largely consistent with the ToE when using 60-year average rates. Supplementary Table 5 shows the consistency of ToE when using a 0-1400 CE background reference period compared to a 0-1700 CE reference period, where the ToE globally and at individual sites remains the same. Additionally, we compared the global ToE using different datasets, including earlier versions of the proxy sea-level database from Kopp et al. 9 and Kemp et al. 40 (Supplementary Table 6). Since there is an abundance of data in the North Atlantic and western North Atlantic compared to the rest of the globe, we also compare the global ToE using only North Atlantic sites, and then all the data excluding North Atlantic sites, and all the data excluding western North Atlantic sites. In all scenarios, the global ToE still occurs in the mid to late 19th century and varies by 29 years. Finally, we tested the influence of proxy data RSL and chronological uncertainties on ToE by reducing uncertainties at European sites to be comparable to the uncertainties of proxy data on the North American coast. In this case, the ToE at European sites remained broadly the same (Supplementary Table 7).

Data availability
Data related to this study are provided in the Supplementary Information and Source Data file. Source data are provided with this paper.

Code availability
Code for the spatiotemporal model results that are reported in the paper 65