Extreme sea levels at different global warming levels

1Joint Global Change Research Institute, Pacific Northwest National Laboratory, College Park, MD, USA. 2Department of Coastal and Urban Risk & Resilience, IHE Delft institute for Water Education, Delft, the Netherlands. 3Harbour, Coastal and Offshore Engineering, Deltares, Delft, the Netherlands. 4Water Engineering and Management, Faculty of Engineering Technology, University of Twente, Enschede, the Netherlands. 5European Commission, Joint Research Centre (JRC), Ispra, Italy. 6Princeton School of Public and International Affairs, Princeton University, Princeton, NJ, USA. 7Department of Atmospheric Sciences, University of Illinois, Urbana-Champaign, IL, USA. 8Department of Infrastructure Engineering, University of Melbourne, Melbourne, Victoria, Australia. 9Rutgers Institute of Earth, Ocean, and Atmospheric Sciences and Department of Earth and Planetary Sciences, Rutgers University, New Brunswick, NJ, USA. 10Present address: Department of Physics and Astronomy Augusto Righi, University of Bologna, Bologna, Italy. ✉e-mail: claudia.tebaldi@pnnl.gov Extreme sea levels (ESLs)1 are triggered by the combination of storm surges, tides and waves. At vulnerable locations, high ESLs can constitute severe hazards, causing extensive damages to both human settlements and coastal ecosystems when natural and engineered defenses are overtopped or breached 2–6. Inevitably, increasing global temperatures will continue to cause global mean sea-level rise7. Even without the potential effects of climate change on waves and storm surge, sea-level rise alone is expected to lead to increases in coastal flooding and/or erosion8–11. Several recent studies or assessments that address coastal flooding at the global scale have characterized future evolution of ESLs, particularly as driven by relative sea-level change (RSLC)1,12. Some of these adopted a scenario perspective, analysing projected changes under future emission trajectories for various time horizons3,5,13–18. Others have used idealized projections of RSLC11. Only a few studies, some regional in scale, have characterized the effects of reaching alternative global warming levels (GWLs) independently of specific scenarios of future emissions19–22. These focused on the low end of the GWL range, in accordance with what motivated the Paris Agreement, that is, the measurable differential impacts for a world at 1.5 °C compared to one at 2 °C above pre-industrial levels23,24. Here, we take this GWL perspective but model RSLC for a wider range of values, from strong mitigation (1.5 °C) to high-end global warming (5 °C) and analyse consequences for ESL future frequencies. We note, however, that while the GWL framing is a natural fit for variables for which (at least on a century timescale) climatological changes with warming are largely time independent, for sea-level change the effects of different GWLs cannot be examined in a time-independent manner, as global mean sea level is more closely related to time-integrated temperature change than to the value of temperature at a given point in time25. Thus, sea level at a time when a given degree of warming is reached can be very different from sea level after the system has had time to equilibrate at that same warming level. We therefore also choose 2100 as a the time horizon by which different warming levels are achieved. We strongly emphasize epistemic uncertainty in the modelling of RSLC, using two methods for deriving its projections. One of them builds on aggregating existing results obtained by representative concentration pathway (RCP)-driven Earth System Model simulations according to their behaviour at 2100 (the GWL they reach, independently of the emission scenario followed, and not necessarily stabilizing at it, in most cases), using an ensemble of projections to account for uncertainties associated with parameters such as climate sensitivity, ocean heat uptake, glacier response and dynamic ocean response19. The second uses a simple model and combines information from Coupled Model Intercomparison Project Phase 5 (CMIP5) RCP-based projections of global temperature together with other IPCC-assessed and/or observationally constrained ranges of those same critical processes in its comprehensive uncertainty exploration26–28. We also include the effects of deeply uncertain factors contributing to ice-sheet changes29 for the 2 and 5 °C GWLs using results from a Structured Expert Judgment (SEJ) study30 (Methods). For present-day ESLs at individual locations, our analysis also adopts a multimethod approach and gathers up to three estimates that used either observational records19 or hydrodynamic modelling16,17. The latter enables us to produce results for a large number of locations covering most of the world’s coastlines, and in one case16 also models the effects of tropical cyclones potentially making landfall at exposed locations. We pair the alternative estimates of current ESLs with the two alternative RSLC projections, producing six alternative fully probabilistic projections of future ESLs for a subset of Extreme sea levels at different global warming levels

xtreme sea levels (ESLs) 1 are triggered by the combination of storm surges, tides and waves. At vulnerable locations, high ESLs can constitute severe hazards, causing extensive damages to both human settlements and coastal ecosystems when natural and engineered defenses are overtopped or breached [2][3][4][5][6] . Inevitably, increasing global temperatures will continue to cause global mean sea-level rise 7 . Even without the potential effects of climate change on waves and storm surge, sea-level rise alone is expected to lead to increases in coastal flooding and/or erosion [8][9][10][11] .
Several recent studies or assessments that address coastal flooding at the global scale have characterized future evolution of ESLs, particularly as driven by relative sea-level change (RSLC) 1,12 . Some of these adopted a scenario perspective, analysing projected changes under future emission trajectories for various time horizons 3,5,[13][14][15][16][17][18] . Others have used idealized projections of RSLC 11 .
Only a few studies, some regional in scale, have characterized the effects of reaching alternative global warming levels (GWLs) independently of specific scenarios of future emissions [19][20][21][22] . These focused on the low end of the GWL range, in accordance with what motivated the Paris Agreement, that is, the measurable differential impacts for a world at 1.5 °C compared to one at 2 °C above pre-industrial levels 23,24 . Here, we take this GWL perspective but model RSLC for a wider range of values, from strong mitigation (1.5 °C) to high-end global warming (5 °C) and analyse consequences for ESL future frequencies. We note, however, that while the GWL framing is a natural fit for variables for which (at least on a century timescale) climatological changes with warming are largely time independent, for sea-level change the effects of different GWLs cannot be examined in a time-independent manner, as global mean sea level is more closely related to time-integrated temperature change than to the value of temperature at a given point in time 25 . Thus, sea level at a time when a given degree of warming is reached can be very different from sea level after the system has had time to equilibrate at that same warming level. We therefore also choose 2100 as a the time horizon by which different warming levels are achieved. We strongly emphasize epistemic uncertainty in the modelling of RSLC, using two methods for deriving its projections. One of them builds on aggregating existing results obtained by representative concentration pathway (RCP)-driven Earth System Model simulations according to their behaviour at 2100 (the GWL they reach, independently of the emission scenario followed, and not necessarily stabilizing at it, in most cases), using an ensemble of projections to account for uncertainties associated with parameters such as climate sensitivity, ocean heat uptake, glacier response and dynamic ocean response 19 . The second uses a simple model and combines information from Coupled Model Intercomparison Project Phase 5 (CMIP5) RCP-based projections of global temperature together with other IPCC-assessed and/or observationally constrained ranges of those same critical processes in its comprehensive uncertainty exploration [26][27][28] . We also include the effects of deeply uncertain factors contributing to ice-sheet changes 29 for the 2 and 5 °C GWLs using results from a Structured Expert Judgment (SEJ) study 30 (Methods).
For present-day ESLs at individual locations, our analysis also adopts a multimethod approach and gathers up to three estimates that used either observational records 19 or hydrodynamic modelling 16,17 . The latter enables us to produce results for a large number of locations covering most of the world's coastlines, and in one case 16 also models the effects of tropical cyclones potentially making landfall at exposed locations. We pair the alternative estimates of current ESLs with the two alternative RSLC projections, producing six alternative fully probabilistic projections of future ESLs for a subset of Extreme sea levels at different global warming levels Claudia Tebaldi 1 ✉ , Roshanka Ranasinghe 2,3,4 , Michalis Vousdoukas 5 , D. J. Rasmussen 6 , Ben Vega-Westhoff 7 , Ebru Kirezci 8 , Robert E. Kopp 9 , Ryan Sriver 7 and Lorenzo Mentaschi 5,10 The Paris agreement focused global climate mitigation policy on limiting global warming to 1.5 or 2 °C above pre-industrial levels. Consequently, projections of hazards and risk are increasingly framed in terms of global warming levels rather than emission scenarios. Here, we use a multimethod approach to describe changes in extreme sea levels driven by changes in mean sea level associated with a wide range of global warming levels, from 1.5 to 5 °C, and for a large number of locations, providing uniform coverage over most of the world's coastlines. We estimate that by 2100 ~50% of the 7,000+ locations considered will experience the present-day 100-yr extreme-sea-level event at least once a year, even under 1.5 °C of warming, and often well before the end of the century. The tropics appear more sensitive than the Northern high latitudes, where some locations do not see this frequency change even for the highest global warming levels.
179 locations and four alternative projections at a larger set of 7,283 locations (when not using the observationally based current ESL estimates). We introduce a voting system, treating the alternative projections as individual experts and producing a central estimate, derived by asking for a majority vote, a lower-pessimistic-bound, akin to taking the most pessimistic projection, and an upper bound, which requires unanimity (Methods). The voting system yields a broad distribution comparable to that produced by 'p-box'-based possibilistic approaches (for example, ref. 31 ) and would appear better suited at respecting the disparities in the alternative estimates we seek to merge. In addition, however, we compare the results of this synthesis approach to a more conventional full convolution of the four, or six, alternatives into a single distribution.
As in several previous studies 10,14,19,32,33 , future ESL changes are here taken to be solely a function of changes in mean sea level at a location, translating essentially into a uniform shift of the present-day ESL distribution. One of the methods considered here for the estimation of ESLs 16 also modelled future changes in storm surge and waves but found that these contributions did not substantially affect the temporal dynamics of ESLs, especially at the global scale adopted herein, compared to the overwhelming effect of RSLC. Similarly, the focus on long-term changes and the global perspective made us choose not to address explicitly natural variability of ESLs, even if localized, mostly present-day analyses have shown its importance 10,34 . The expectation is that the size of RSLC will overwhelm such variability by the end of the century. Of course, even if limited by the length of the historical period used for estimating present-day ESLs (≥30 yr depending on the method), a portion of that internal variability is necessarily reflected in the ESL central estimates and their uncertainty bound.
Our focus here is on large-scale global assessment and we position our study as eminently relevant for scientists and policy-makers who are interested in the spatial dynamics of sea-level rise and its implications. Local stakeholders have to contend with additional uncertainties and components of impact risk, like exposure and vulnerability linked to topography and future socioeconomic conditions, for example. We do not address those here, as we focus our analysis on changes in a standard benchmark metric of hazards, the 100-yr ESL event (the extreme total water level expected to be experienced on average once in 100 yr or with 0.01 probability every year), recognizing that it may or may not be of immediate salience to individual localities.

gWLs triggering changes in frequency of ESLs
The study is framed to answer a specific question: what is the lowest warming level to trigger at least a 100-fold change in frequency of the present-day 100-yr ESL event, making it at least an annual occurrence by 2100?
We start from the subset of 179 locations at which we have six alternative estimates. For all components, RSLCs and ESLs, full probability distributions are available (Methods); in Supplementary Figs. 1-3 we show the three alternative mean estimates of the present-day 100-yr ESL and the associated standard deviations. Supplementary Figs. 4-11 show the six GWL-based RSLC median projections available from the two alternative methods, plus the two projections that include ice-sheet contributions based on SEJ available for one of the methods. In Extended Data Fig. 1, we also show the mean difference across the three ESL datasets between the present-day 100-yr and 1-yr ESL events, together with the standard deviation around it. The latter provides important context against which to evaluate the salience of the results we present below, since at some locations the difference is much smaller than at others. By design, our analysis does not extend to the evaluation of actual, location-specific impacts. Therefore, in the following, we do not highlight how large a change in the magnitude of the annual event at 2100 the enhanced frequency implies.
The central estimate of the GWL triggering the change in ESL frequency obtained by the voting system (akin to a majority vote, as described in the Methods) is shown in Fig. 1a. Already at 1.5 °C of warming by 2100 a large number of locations (54%; also Table 1, first row) will experience their present-day 100-yr ESL events annually (or more often). This is overwhelmingly true for locations in the southern hemisphere and the subtropics (Hawaii, the Maritime continent and the Caribbean) and the southern half of North America's Pacific coast. However, a consideration of the uncertainties (Methods) shows that this frequency change would occur at up to 99% of locations in the most pessimistic lower bound produced by the voting system ( Fig. 1c and Table 1, second row), while occurring at only 2% of locations for the optimistic upper bound ( Fig. 1e and Table 1, third row). These bounding outcomes, akin to only asking for the vote of the single lowest GWL estimate, in the case of the pessimistic bound, or to requiring a unanimous vote in the case of the optimistic bound (as explained in the Methods) highlight the substantial level of disagreement among the six estimates, stemming from the wide ranges of the estimated present-day ESLs and the RSLC projections. Further considering the central estimate ( Fig. 1a and Table 1, first row), an additional 14% of locations undergoes the change in frequency if global warming reaches 2 °C by 2100. But up to 20% of locations do not reach that annual frequency even for a warming of 5 °C and even if the SEJ-derived effects of ice-sheet dynamics are included (as all purple dots signify). Most of this last type of location are along the coastlines of Alaska and Northern Europe but are also found in regions exposed to tropical cyclone activity. Exposure to tropical cyclone activity makes the present-day ESL distributions wider (Extended Data Fig. 1) and therefore less sensitive to uniform shifts from RSLC. Further consideration of the upper bound shows that up to 80% of locations do not see the 100-fold increase in frequency for any of the considered GWLs, suggesting that at least one of the six estimates produces a more moderate frequency change even at the highest GWL considered.
Results for a much larger number of locations (7,283) at which four alternative projections are available and are combined through our voting system are presented in Fig. 1b,d,f and Table 1 (fourth to sixth rows). These have almost uniform global coverage, with an along-coast spacing of ~1° (~100 km at the equator; Methods). The differences between 100-yr and 1-yr events at this larger set of locations are shown in Extended Data Fig. 2.
The spatial dynamics shown by the analysis at the sparse subset of the tide-gauge locations (Fig. 1a,c,e and Table 1, first through third rows) appears to be representative of that along the global coastline ( Fig. 1b,d,f and Table 1, fourth through sixth rows). Most of the tropics see the change in frequency from 100-yr to 1-yr (or more frequent) ESL event already at 1.5 °C of warming also for the larger set of locations, while the Northern Hemisphere high latitudes are again projected not to experience such a dramatic ESL frequency change, even for the highest level of warming, and even with SEJ-derived ice-sheet contributions (as all purple dots signify). The greater coverage indicates that parts of the Mediterranean coasts and the Arabian Peninsula may also become ESL hotspots. A more complex picture emerges for some of the coastal areas, with a high level of along-coast variability; for example, along the Atlantic coast of North America, where the whole range of GWLs can result in the 100-fold change in frequency at nearby locations. The West coast of North America appears less sensitive to GWLs up to 5 °C and so does the Pacific coast of the Asian continent. Japan appears to experience a dichotomy between a sensitive West coast and its East coast, showing to endure the whole range of GWLs. Most of these features are reflected in the mean and standard deviation of the differences between 100-yr and 1-yr events shown in Supplementary  Fig. 13. All this detail, however, does not translate in radically different distributions of the proportion of sites under the triggering GWLs, as a comparison of the corresponding rows in Table 1 shows.
For the larger set of locations, the central estimate projects that 43% of locations will experience the present-day 100-yr ESL event at an annual or higher frequency even at 1.5 °C and that such a large frequency change will not occur at about 23% of the locations, even at the highest warming level (and even when including the effects of ice-sheet melt). The pessimistic lower bound estimate shows 99% of the locations experiencing the frequency change at 1.5 °C. The optimistic upper bound has only 7% experiencing this change at 1.5 °C and 68% avoiding the change altogether.
For both sets of locations, a common feature of the two central estimates of the distributions of sites (the relative magnitude of the numbers along the tables' first and fourth rows), is their U-shape. Most of the sites either see the 100-fold increase in frequency at the lowest levels of warming of 1.5 and 2 °C, or do so only-if at all-for the highest level of warming, 5 °C and only when that includes the effects of ice-sheet dynamics. This is visible in the histograms in Extended Data Fig. 3, left-hand side, where the height of the bars correspond to the content of the two sets of rows in Table 1. Thus, the voting system outcomes appear to be less sensitive to the range of GWL above the lower ones.
We test the robustness of this result by considering the outcomes of an alternative synthesis method, a full convolution that constructs a unique distribution of future ESLs from the six (or four) distributions that we have treated as individual expert voters in our approach so far. The full convolution could be seen as a more traditional approach at merging different probabilistic estimates but given the substantial differences in the individual estimates we propose our voting system as better suited at respecting such disagreements. Focusing on central estimates and lower and upper bounds available from both methods, a comparison of the histograms in Extended Data Fig. 3, left versus right column, indicates that the choice of synthesis method does not impact the lower bound (pessimistic) estimates (histograms in the third and fourth row) but changes the central estimate (top two rows) and the upper bound (bottom two rows) in two ways: by slightly increasing the fraction of sites for which the increase in ESL event frequency occurs at the lowest GWL (for the central estimate) or does not happen at all (for the upper bound) and by spreading more evenly across the other GWLs the remaining fraction of sites. Thus, a full-convolution approach erases the U-shape of the voting approach and creates a more uniformly distributed behaviour between the two extremes of the GWL range. (The tractability of a single probability distribution rather than the need of reconciling four or six of them through our voting system, allows us to gain additional insights, by filling in intermediate quantiles to better represent the range of probabilistic outcomes. The distribution of GWLs including two intermediate quantiles, bounding 66% of probability, is shown in Supplementary Tables 4 and 5. As expected, the outcomes across GWLs are in these intermediate cases more evenly spread.) Spatially, the differences emerging between the two approaches do not change the large-scale geographical patterns but make the coastal regions, other than the hotspots, show a more gradual sensitivity to increasing GWLs in the convolution approach. Nonetheless, this method confirms that the Atlantic coast of North America and the coasts of Northeast Asia are affected by a high degree of along-coast variability. These results are mapped in Extended Data Fig. 4.

Timing of the ESL frequency changes
Last, even if our analyses focussed on what happens at the time when the discrete set of warming levels are reached (2100), the availability of the corresponding RSLC timeseries over the twenty-first century (2020-2100) allows us to answer a further question, about when 100-yr to 1-yr ESL frequency changes are first observed. We stress, however, that the answer needs to be always conditional on the trajectories identified as being consistent with the individual GWLs. Supplementary Figs. 12 (for 1.5 °C) through 19 (for 5.0 °C) and Table 2 show the decade when the 100-fold frequency change is first observed, for the larger set of locations, and using the full-convolution approach. For the locations that are sensitive to even the lowest warming level of 1.5 °C-identified in our previous analysis mostly with the tropics and subtropics, that is, all locations not coloured in purple in Supplementary Fig. 12

Discussion
We use a voting system, which we also compare to a more traditional full-convolution approach to synthesize alternative projections of ESL frequency changes (from 100-yr event at present to annual or more frequent by 2100) at a range of GWLs, for a large number of locations all along the world's coastlines.
According to our central estimate, by the end of the century close to half of the locations considered will experience the present-day 100-yr ESL event at least once a year, even for a trajectory of global temperature that limits warming to 1.5 °C or at most 2 °C. The lower bound estimate sees practically all sites (98 or 99%) experiencing that dramatic change already for 1.5 °C. On the opposite end of the spectrum, an optimistic estimate sees about 70-80% of the locations avoiding such an increase in frequency altogether. Locations around the world can largely be distinguished according to a bimodal behaviour with respect to this change in frequency, either showing large sensitivities already at the lowest warming levels of 1.5 or 2 °C or being insensitive even up to the highest warming considered (5 °C), as about 20% of locations do not experience the frequency change, even when the possibility of rapid ice-sheet loss is included in the estimates of RSLC. We tested the sensitivity of this result to the approach taken to determine our 'consensus' by applying a more traditional convolution of the alternative probability distributions. If the bimodality disappears, the overall message does not change substantially. According to this last approach, we were also able to estimate the time at which the 100-fold increase in ESL frequency happens and we find that most of these locations will experience such change earlier than the end of the century, that is, in the decade between 2070 and 2080 according to our central estimate under the lowest GWL of 1.5 °C, one decade earlier for the intermediate GWLs and as early as 2060 for the highest GWLs of 4.0 and 5.0 °C. These last results are found to be very consistent with the analogous results shown in Fig. SPM4 of ref. 35 , if one considers RCP 2.6 as a 2.0 °C end-of-century scenario and RCP 8.5 as a 4.0-5.0 °C end-of-century scenario 36 .
Coastlines in the tropics and parts of the Mediterranean and the Arabian Peninsula appear to be the places where these 100-fold changes in frequency of the 100-yr ESL events will take place even Percentage of the 179 locations depicted in Fig. 1a,c,e (first to third rows of the table) and of the 7,283 locations depicted in Fig. 1b,d,f (fourth to sixth rows of the table) at which the frequency of the present-day 100-yr ESL event changes to at least an annual event by 2100. Indicated are central estimates, lower (pessimistic) and upper (optimistic) bounds for the GWL required. Percentages along each row may not add up to 100 exactly because of rounding errors. Each percentage value under GWLs of 2 °C or higher is to be interpreted as the additional fraction of sites experiencing the change (that is, in addition to the total along the row to its left). The + sign associated with 2 and 5 °C indicates projections that include SEJ-derived estimates of ice-sheet contribution to RSLC.
if 2100 warming is limited to 1.5 or 2 °C, consistent with previous studies that focused specifically on these lower levels of warming 37 . The coastal regions of the highest latitudes of both hemispheres are those where even higher levels of warming will not produce such frequency increases. While there may be localized exceptions, this means that, in terms of a 100-fold increase in the 100-yr ESL frequency, microtidal areas are highly sensitive to even the smaller GWLs considered here, while meso-and macro-tidal areas are not. The same applies for coastal areas that are either protected or exposed to meteorological extremes, with the former being more vulnerable to changing mean sea levels that expose them to unprecedented ESLs. The RSLC maps in the Supplementary Information also show that sea level is projected to increase above-average in many of these regions. The Atlantic coast of North America and the coast of Northeast Asia appear to host a high degree of spatial variability in the results, with locations adjacent to one another either experiencing a large change in frequency at very low warming levels, or not experiencing it even for the highest warming levels considered here, calling for more indepth analysis and detailed modelling of the local dynamics. These geographically differentiated results are consistent with the analyses of refs. 11,16 , which also highlighted similar differential sensitivities for return period changes due to RSLC. Topography and the history of extreme events experienced in the record at these locations are probably the source of such variations.
Our findings have important policy and practical implications as they highlight that even if the Paris Agreement goals will be achieved, extreme events potentially conducive to coastal flooding will be experienced at unprecedented frequencies in many parts of the world's coasts. Only a very strict unanimity rule required of the alternative estimates produces substantially more optimistic projections. Our analysis was cast in terms of return period changes as return periods are typically the basis for coastal protection infrastructure design and hazard communication and therefore underline the potential need for global-scale adaptation efforts, even if our global assessment cannot match the detail and accuracy of local/regional studies. However, apart from allowing an indepth analysis of the spatial dynamics of sea-level rise using a globally homogeneous framework, our results allow the identification of ESL change 'hotspots' for which detailed dynamical modelling and a more indepth interpretation of the impact of changes in ESL hazard on human and natural systems [38][39][40][41] are needed to serve local decisions.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41558-021-01127-1.  We underline the first fraction to become larger than 0.5 to highlight the decade by which more than half the sites considered experience such change. The last column highlights the fraction of the sites that do not experience the 100-fold change by 2100 for that GWL.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Methods
Our analysis focuses on frequency changes of present-day 100-yr ESL events (events having 0.01 probability in any given year) at locations along the world's ice-free coastline (Greenland, the Russian Arctic coast and Antarctica excluded).
At each location, future probabilistic estimates of RSLC are added to present-day probability distribution of ESLs. The addition has the first-order effect of shifting the distribution of ESLs uniformly and therefore changing the height of the 100-year event (and any other) by an offset equal to the amount of RSLC. The addition of an offset also changes the return period (the expected frequency) of the present-day 100-yr event, by making it more frequent if RSLC is positive or less frequent in the few locations and time periods where RSLC projections are negative. While the change in height of the future 100-yr event is by construction in our analysis the magnitude of RSLC, the change in frequency will depend on the interplay between the RSLC magnitude and the shape of the present-day ESL distribution. Importantly, all uncertainties in the parameters of the ESL distribution and the RSLC projections are taken into account.
All three studies used here for the ESL component fit peak-over-threshold Poisson-Generalized Pareto Distribution (GPD) models 42 to present-day ESLs at each location. One of the original studies looked at alternative functional forms for the fit and concluded that the GPD had optimal properties 17 . Depending on the study, these ESLs are either observed from tide gauges or generated by hydrodynamic models that are validated with observations. Each Poisson-GPD model is characterized by four parameters: a threshold, determining the location of the distribution; a scale parameter, determining its width; a shape parameter, determining how fast the tail declines; and a rate parameter, determining the expected frequency of threshold exceedances. The addition of a positive amount of RSLC can be equivalently characterized as a decrease in the threshold by the same amount and therefore an increase in the likelihood of exceedances of a given height. The magnitude of the increase in the likelihood of occurrence of a given ESL caused by a given shift is determined by the scale and shape parameters.
If a quantity-in our case the extreme total water level 1 , which we here call ESL-follows the GPD, its probability distribution can be expressed as a function of three parameters: a threshold μ, a scale parameter σ and a shape parameter ξ.
is the GPD density function. The additional parameter of interest that needs to be estimated to compute return periods is λ the Poisson rate of exceedance, indicating the expected annual number of events exceeding the threshold μ.
According to the GPD definition, the expected number of events exceeding the threshold in a year is then In our study, when we consider RSLC corresponding to year 2100 GWLs of 1.5, 2, 2.5, 3, 4 or 5°C above pre-industrial levels, we assume that those levels are reached by a gradual, non-decreasing trajectory of temperature change over the course of the twenty-first century.
Present-day ESL estimates from three previous studies. The three previously published studies that we use provide GPD parameters at locations around the world using the following approaches. . This study 16 uses the baseline period of 1980-2014 to generate realizations of the additive components of the total water levels (tidal elevation, wave setup and storm surges) by forcing a storm surge model and a wave model (Delft3D-FM (refs. 43,44 ) and WW3 (refs. 45,46 ) respectively) with ERA-INTERIM wind and pressure fields. Tropical cyclones effects are included in this baseline estimates by using satellite data (for wave heights) and simulating all historical tropical cyclones recorded during the 1980-2014 period (for storm surge) on the basis of the datasets available through the Hurricane Research Division of the National Oceanic and Atmospheric Administration (United States), the Joint Typhoon Warning Center and the UNISYS database (Data availability). These weather-driven components are added to the high-tide level at each location. A GPD is then fitted to these data. While here we use only the baseline period GPD parameters at each location, ref. 16 also estimate future changes in all the ESL components, by using CMIP5 model output and two emission scenarios to force the hydrodynamic models, but find that at most locations the main driver of change remains RSLC. The study relies on the preceding validation of the individual models used. Kirezci et al. (2020). This study 17 combines global models of tide (FES2014, ref. 47 ), storm surge (GTSR, ref. 43 ) and waves (GOW2, ref. 48 ) to reconstruct the historical total water level at the DIVA coastal locations 49,50 for the baseline period 1979-2014. Total water levels are defined as above, as the linear summation of tide, surge and wave setup. ESLs are determined by fitting ten different extreme value analysis methods and the global best fit is determined to be a GPD using the 98th-percentile threshold. Both the historical total water levels and the ESLs are rigorously validated against the quasi-global GESLA-2 tide-gauge observations 51 . The study does not consider possible future changes in storm surge heights and wave heights. Tropical cyclones effects are not included in these ESL estimates. Rasmussen et al. (2018). This study 19 , in contrast to the two previously described, relies only on observed hourly records of still water height to fit GPDs at a global network of about 200 tide gauges (University of Hawaii Sea Level Center (Data availability)). Only tide-gauge records of length >30 consecutive years and for which each year has >80% of observations available are considered. Unlike the above model-based approaches, this would neglect wave setup and swash contributions, both of which can be important contributors to ESLs 52,53 , unless the gauges are located nearshore. For each day in a given tide-gauge record with >12 h of data, the daily maximum sea level is estimated. Following ref. 32 , the variation in ESLs is isolated by subtracting the annual mean sea level change from each daily maximum value (that is, values are detrended). The detrended daily maximum tide values are then referenced to local mean sea level. Daily maximum sea levels that are (1) above the 99th percentile and (2) within 3 d of each other are declustered to meet the statistical independence assumption of the GPD. The GPD is fitted to these daily maxima. Here, as in the previous two studies, future changes in storm frequency 54 , intensity 55 or track 56 , which could all modify the GPD parameters, are not considered. Also not considered are changes in the tide-surge interaction 53,57,58 or future changes in geomorphology, which can both impact the return periods of ESLs 59,60 .
Future RSLC estimates based on two previous studies. Our projections of future RSLC, geographically detailed, are obtained according to two methods developed in previously published studies but here applied to a wider range of GWLs. Rasmussen et al. (2018). This approach 19 uses a set of local, probabilistic, RSLC projections conditional on each temperature target. Projections are made at a global tide-gauge network (Permanent Service for Mean Sea Level (Data availability)) as well as at the centre points of a 1° × 1° grid covering the coastlines of nearly all global landmasses, including several islands. The methodology of the projections follows the 'bottom-up' , component-based framework of ref. 61 , with modifications to accommodate temperature targets 19,62 and estimates of dynamic ice-sheet melt from an SEJ exercise 30 .
Following ref. 19 , temperature targets are accommodated by constructing alternative ensembles for each temperature scenario using Atmosphere-Ocean Global Climate Model (AOGCM) output filtered according to each AOGCM's 2100 global mean surface air temperature. Specifically, AOGCM outputs from the CMIP5 archive 63 are used for global mean thermal expansion, local ocean dynamics and as a driver of a surface mass balance model of glaciers and ice caps 64 . We detail the GCMs/RCPs used in Supplementary Tables 1-3. We create ensembles for 1.5, 2, 2.5, 3, 4 and 5 °C with AOGCMs that have a twenty-first-century global mean surface air temperature increase (19-yr running average) of 1.5 ± 0.25 °C, 2 ± 0.25 °C, 2.5 ± 0.25 °C, 3 ± 0.5 °C, 4 ± 0.5 °C and 5 ± 0.5 °C. For consistency with the ref. 61 framework, which models 19-yr running averages of SLC relative to 2000, anomalies of global mean surface air temperature from 1991 to 2009 are computed and then shifted upward by 0.72 °C to account for warming since 1875-1900 65 . The global mean surface air temperature trajectories and global mean sea-level contributions from thermal expansion and glacial ice from selected CMIP5 models that are binned into each temperature target category are shown in Supplementary Figs. 20 and 21 respectively. Supplementary Tables 1-3 list the AOGCMs used in each temperature target ensemble, the RCP and the sea-level components derived.
Two different approaches are used for modelling contributions from ice-sheet melt. In the first approach, Antarctic Ice Sheet (AIS) and Greenland Ice Sheet (GIS) contributions are estimated by combining IPCC's Fifth Assessment Report (AR5) 'likely range' projections of ice-sheet dynamics and surface mass balance (Table  13.5 in ref. 66 ) with information about distribution tail shape from an SEJ study of ice-sheet mass loss from ref. 67 . Total ice-sheet contributions are computed by this approach for each temperature target considered in this study, by randomly sampling the AIS and GIS ice-sheet distribution for each RCP (Table 13.5 in ref. 66 ) in proportion to the representation of each RCP in the groups of CMIP5 models selected for each temperature target. For example, the 1.5 °C target uses 12 CMIP5 models from RCP 2.6 and two from RCP 4.5, so 86% of the samples are drawn from the RCP 2.6 distribution and 14% are drawn from the RCP 4.5 distribution (Supplementary Table 1).
In the second approach, which we apply only to projections for the 2 and 5 °C temperature targets, contributions from ice-sheet melt are obtained from the SEJ study of ref. 30 . Projections of RSLC after mid-twenty-first century are highly dependent on ice-sheet melt because of its potential for substantial contributions to global mean sea-level rise 5,68 . However, incomplete understanding of the physical processes that govern ice-sheet melt inhibits realistic representations in process-based models. In such cases of incomplete scientific understanding, SEJ using calibrated expert responses is one approach for estimating such uncertain quantities (as used here). Note that we treat these two as additional levels of the GWL range and pair them with the same 2 and 5 °C estimates from the alternative method (not specifically containing SEJ ice-sheet melt estimates) when synthesizing through the voting or the convolution methods. -Westhoff et al. (2019). RSLC projections from this approach 27 are from a recent perturbed physics ensemble using the reduced complexity climate model Hector-BRICK 26,27,69 . The model includes a one-dimensional diffusive heat and energy balance model, combined with a sea-level change module that represents contributions from thermal expansion, glaciers and small ice caps and polar ice sheets. The model also includes a simple parameterization of fast ice-sheet dynamic disintegration in the Antarctic ice-sheet component 26 . We estimate 39 uncertain parameters using a Bayesian calibration method (adaptive Markov Chain Monte Carlo), with observational constraints from global surface temperature, global ocean thermal expansion and polar land ice. The ensemble and applications are discussed in more detail in refs. 27,28 .
Hector-BRICK simulates global sea-level rise, so gridded and local tide-gauge estimates (Permanent Service for Mean Sea Level (Data availability)) are constructed from Hector-BRICK results using a fingerprinting technique that accounts for regional variability due to mass redistribution from melting land ice 28,70 . For this fingerprinting, we assume land water storage and ocean thermal expansion are globally uniform and we neglect changes in local ocean dynamics. We also assume that contributions from Antarctica and Greenland are uniform over their respective ice-sheet surfaces. We also neglect the possibility of different paces of mass change for different glacier regions over the twenty-first century.
Both of the methods produce RSLC projections at the same set of tide gauges used in ref. 19 and in addition cover a regular 1° × 1° grid over the world's coastlines. Both methods produce projections in the form of an empirical distribution of values generated through Monte Carlo sampling.
The methods can produce substantially different projections, especially in the high latitudes . We attribute these large differences to two reasons. First, the methods use two different approaches for accounting for non-climatic vertical land motion. Following ref. 61 , ref. 19 uses a spatiotemporal Gaussian process model to account for background rates of vertical land motion. This includes deltaic processes, tectonic uplift, anthropogenic subsidence processes (for example, ground water and fossil fuel withdrawal) and glacial isostatic adjustment. On the other hand, ref. 27 only considers glacial isostatic adjustment estimates from ref. 70 . Second, the two methods consider different estimates of GIS melt and AIS melt. For example, ref. 27 considers a median estimate for AIS melt that is nearly 0.5 m larger than for ref. 19 . Greater ice-sheet melt can amplify differences in RSLC due to gravitational-rotational and deformational effects 1 .
Matching the two components and computing future frequencies. The GPDs estimated by each of the three studies include measures of standard errors for the parameters and the 100-yr event. To sample from the distributions fully accounting for the uncertainty in the parameters, we use an asymptotic approximation to their covariance matrix based on the Fisher information matrix, as presented in refs. 71 . The Fisher information matrix, when calculated at the maximum likelihood estimates of the parameters, can be interpreted as the inverse of the covariance matrix estimate and in our notation takes the simple form where N is the expected total number of exceedances, estimated by λ multiplied by the number of years used in the estimation. The approximation is valid as long as ξ > −0.5 which is true for all but a handful of locations, which we discard.
We therefore resample at each location pairs of the shape and scale parameters according to a normal distribution with mean their maximum likelihood estimates and covariance matrix (equation (3)), and for each pair we derive the corresponding value of the 100-yr event. We then match them with as many random samples from the RSLC distributions.
For each location and each sample i, given the magnitude of the present-day 100-year event, x i , the four parameters of the GPD, μ, σ i , ξ i and λ and the magnitude of RSLC by 2100 for a given warming target, s i , the new frequency RP(x i ) (RP, return period) can be computed as if xi ≤ μ ′ else: where μ ′ = μ + si. A population of samples large enough to characterize future statistics robustly (we use 1,000 after verifying that larger samples would not change substantially our results) is computed separately for each of the six combinations of the three sets of ESL estimates and the two sets of RSLC estimates.
Note that the use of this formula caps the highest obtainable frequency to the 1-year return period. Therefore, by construction that is going to be the highest future frequency projected by our approach, even if some locations could experience present-day 100-year ESL events more often than once a year under different GWLs 19 . Consequently, when throughout the paper we use the expressions 1-yr, or annual, ESL event, we have added wording to the effect that this is a lower estimate of the possible change in frequency.
Synthesizing. For a subset of 179 locations, and for a given warming level, we have six full probability distributions of future (2100) return periods of the current 100-yr ESL event, one for each combination of the three ESL estimates with the two RSLC estimates. For a much larger set of 7,283 locations we have four, as the observation-based ESL estimates from ref. 19 are only available for the subset of 179 common locations.
By considering the evolution of these distributions across the discrete set of GWLs at 2100, it is possible to pinpoint the warming level, if any, at which the current 100-yr event becomes an annual (or more frequent) event. Estimates derived are by construction not deterministic but hinge on choosing a level of confidence, for example, a quantile of the six distributions and some form of 'voting' among them.
We use the following three ways to synthesize this complex set of projections at each location where all six probabilistic estimates are available.
1. For each 2100 warming level, we consider the median value (50th quantile) of each of the six distributions of future return periods and we will take the lowest warming level for which most of the six medians predict a change from 100-yr to 1-yr (that is, at least four of the medians are equal to 1-yr RP). This will be our 'central estimate' of the warming level at which the frequency of the current 100-yr event becomes annual (or higher). This estimate could be characterized as a majority vote (using the medians of the distributions as the individual votes). 2. For each 2100 warming level, we will consider the 5th quantile of each of the the six distributions. This choice represents a 'worst case' because it favours low values in the distribution, translating to a short return period, and therefore a high frequency. Our voting system in this case will pick the most pessimistic of the six estimates; that is, we will choose the lowest warming level at which the minimum of the six 5th quantiles is equal to 1-yr RP. This will be the pessimistic (lower) bound of our estimates. This estimate can be characterized as taking the most extreme (in a negative sense) of the alternative estimates; that is, using a low value of the most sensitive distribution. 3. For each 2100 warming level, we will consider the 95th quantile of each of the six distributions. Mirroring the pessimistic case above, this choice represents a 'best case' because it favours large values in the distribution, translating to a long return period and therefore a lower frequency. Our voting system in this case will pick the lowest warming level at which the maximum of the six 95th quantiles is 1-yr RP. This will constitute our optimistic upper bound estimate. Note that by definition requiring the maximum to be 1 is equal to requiring that all six 95th quantiles are equal to 1; that is, we require unanimity among the individual votes.
We apply this synthesis to the 179 locations at which all three ESL estimates are available (we define one available if it produces an estimate at a location within a distance less than a (spatial) degree, which is also the criterion for matching RSLC and ESLs estimates). We then show the same synthesis (with the required modification in defining the majority vote out of four) applied to the 7,283 locations at which both ESL estimates based on hydrodynamic modelling are available but not the observation-based one.
An alternative, more conventional, approach to synthesizing these multiple distributions is to consider them equally likely and combine them by Monte Carlo sampling into a single full-convolution distribution, for which median and quantiles of interest can be extracted. We perform this alternative approach as well, by sampling 1,000 values from the individual distribution with equal weights and discuss differences in the outcomes of the two synthesis approaches.