Increase in sea level variability with ocean warming associated with the nonlinear thermal expansion of seawater

Sea level variability increasingly contributes to coastal flooding and erosion as global sea levels rise, partly due to the thermal expansion of seawater, which accelerates with increasing temperature. Climate model simulations with increasing greenhouse gas emissions suggest that future sea level variability, such as the annual and interannual oscillations that alter local astronomical tidal cycles and contribute to coastal impacts, will also increase in many regions. Here, we present an analysis of the CMIP5 climate model projections of future sea level to show that there is a tendency for a near-global increase in sea level variability with continued warming that is robust across models, regardless of whether ocean temperature variability increases. Specifically, for an upper-ocean warming by 2 °C, which is likely to be reached by the end of this century, sea level variability increases by 4 to 10% globally on seasonal-to-interannual timescales because of the nonlinear thermal expansion of seawater. As the oceans continue to warm, future ocean temperature oscillations will cause increasingly larger buoyancy-related sea level fluctuations that may alter coastal risks. Sea level variability will increase by up to 10% globally on seasonal-to-interannual timescales for future upper-ocean warming of about 2 °C, because thermal expansion accelerates at higher temperatures, suggest analyses of climate model simulations.

G lobal sea levels are rising 1 in part due to the expansion of seawater with increasing temperatures 2,3 . Since the rate, or coefficient, of thermal expansion increases with greenhouse warming 4,5 , seawater density becomes increasingly sensitive to higher temperatures; thereby contributing to observed and future projections of accelerating sea level rise 1,6,7 . So far unexplored is how the nonlinear thermal expansion property of seawater will affect the variability of future higher sea levels.
Variations in coastal sea levels are already causing more frequent flooding and erosion due to increasing sea level rise [8][9][10] . Regionally, and on seasonal-to-interannual timescales, the sea level variability is mostly determined by the ocean temperature structure, and hence the densities in the seawater column below 11,12 . Large regional sea level variations, such as those associated with the El Niño-Southern Oscillation (ENSO; refs. [13][14][15], are linked to wind-driven shifts of the thermocline 16 as well as the oceanic mixed-layer heat content 17 . Many climate models project increased future ocean temperature variability related to more extreme and frequent ENSO events [18][19][20] . As a consequence, temperature-driven sea level variability (i.e., the thermosteric component 21 ) also increases in the tropical Pacific Ocean 22,23 .
Given the increase in the rate of thermal expansion with temperature 4 and the often dominant role of the thermosteric component in explaining sea level variability 11,12,24 , we hypothesize that sea level variability must increase relative to temperature variability in a warming ocean. Combining the effect of nonlinear thermal expansion with increasing temperature variability, e.g., projected in ENSO-affected areas, the increase in sea level variability must then be even larger than that in temperature variability.
Here, we use an ensemble of climate models to show that there is a near-global tendency for the sea level variability to increase with greenhouse warming. We quantify across these models the components of sea level variability change related to either changes in mean temperature (i.e., the change in the rate of thermal expansion) or changes in temperature variability. Thereby, we assess the inter-model uncertainty of future sea level variability associated with these two components. To explain our methodology, we first discuss the oceanic changes near a sample of coastal cities, then expand our analysis globally. The interpretation of results is supported by analytical solutions to a reduced-gravity ocean model prescribed with thermal expansion characteristics from the climate models. Lastly, we discuss implications for describing future coastal flooding and erosion risks.

Results
Changing sea level variability. Observed sea levels vary seasonally and interannually everywhere, although there are pronounced gradients between regions of larger and smaller variability in sea surface height (SSH; Fig. 1a, b). Some of the largest annual ranges of sea level (10 cm to greater than 30 cm from the minimum to maximum of monthly averages) occur near the continental margins of the northwestern Atlantic and Pacific Oceans, as well as in the northern Indian and the tropical Pacific Oceans (Fig. 1a). For example, in the Florida Strait near the City of Miami, the sea level annual range averages 13 cm ( Fig. 2a; 22 cm at the local tide gauge). The annual maximum of sea level, which typically occurs during September-November in Miami, often determines the season of greatest risk for coastal flooding 25 . Whereas the sea level annual cycle is incorporated into most tidal predictions 26 , interannual variability (Figs. 1b and 2a) that either amplifies or dampens the annual cycle 27 , is not. Should interannual high sea level anomalies ( Fig. 1b; e.g., standard deviations greater than 5 cm) occur during the season of highest sea levels, then coastal flooding is more likely to occur; especially if accompanied by large astronomical tides or exacerbating meteorological events 28 (e.g., storm surges or runoffs from heavy rainfall).
Global climate models are able to simulate many of the salient features of the observed sea level annual range ( Supplementary  Fig. 1) and interannual standard deviation ( Supplementary Fig. 2 Table 2) are 22 and 26 for the tropics and midlatitudes, respectively. There are stark regional differences in the projected increases of the annual range and interannual standard deviation (Fig. 1c, d), with reduced variability also being projected in some areas. Considered regionally, the inter-model uncertainty of future sea level variability changes is substantial, even concerning the sign of changes (stippling in Fig. 1c, d).
Physical processes. One large region where many models do agree on increasing sea level variability is in the tropical Pacific (interannual standard deviation changes of 20%; Fig. 1d), which has been shown to be related to more frequent occurrences of strong El Niño and La Niña events in the future 23 . Increasing ENSO variability [18][19][20] would intensify primarily wind-driven fluctuations of the tropical Pacific thermocline, upper-ocean temperatures, and sea level 12,13,29 . Yet, the future sea level is projected to become more variable interannually in many regions that are not directly affected by ENSO ( Fig. 1d; e.g., in the mid-latitudes). Furthermore, the annual cycle of sea level is also projected to become larger in more areas than not (Fig. 1c), such as near Miami where the range during the 21st century future simulation is 12% larger than during the 20th century historical simulation (Figs. 1c and 2b). For seasonal-tointerannual variability in many regions, ocean temperature fluctuations will directly affect the sea level due to the sensitivity of seawater buoyancy to temperature (i.e., thermal expansion or contraction of the water column), which is the dominant mode of variability 11,12 (i.e., sea level can often be substituted by thermosteric sea level; Methods). At Miami, for example, the annual maximum of sea level usually occurs close to when the upper-100 m ocean temperature (T100, which is a proxy for most of the column-integrated temperature variability that typically occurs above the thermocline at this location 27 ) is the warmest (Fig. 2a, c; October or September peaks). Similar coherence between temperature and sea level cycles 11 are observed for most coastal regions 30,31 , with wind-driven mass redistributions and other processes explaining differences, especially in shallow seas 32 .
For Miami, the CMIP5 average annual range of T100 is projected to increase by 7% (Fig. 2d). If the rate of thermal Projected T100 annual cycle (°C) +7% expansion was constant (and assuming no changes in mass or salinity), the increased temperature variability would translate into a 7% increase in sea level variability (Fig. 3a). However, according to the equation of state (EOS; Methods) of seawater 33 , which describes empirically the ocean density as a nonlinear function of temperature, salinity, and pressure 34 , the thermal expansion coefficient increases with ocean temperature. While the EOS is nonlinear in several ways, it is the relation between the thermal expansion coefficient and temperature that is of most interest here. According to this relation, the density (and therefore sea level) range associated with a temperature fluctuation of the same amplitude will always be larger for warmer (e.g., 21st century) than cooler (e.g., 20th century) mean temperatures, as illustrated in Fig. 3a. Assuming small changes in the thermal expansion coefficient and temperature variability, changes in thermosteric sea level variability can be expanded linearly and are approximately given by the sum of the two contributions (i.e., associated with the EOS and temperature variability; Methods). For example, the future Miami T100 warms by about 2°C (Fig. 2d), increasing the coefficient of thermal expansion by about 5% (Fig. 3a). Together, the change in annual temperature range (7%) and rate of thermal expansion add up to a 12% increase in thermosteric sea level variability, which matches the local CMIP5 increase in sea level variability (Figs. 1c and 2b). All 29 CMIP5 model projections that we assessed for the Miami example location have larger total future annual ranges of seawater density compared to if only temperature variability changes are considered (Fig. 3b). Consistent results are found near the six largest coastal cities in the world ( Supplementary  Fig. 3), which together sample seawater density changes through parts of the tropics and mid-latitudes in the Pacific, Indian, and Atlantic Oceans. Whereas the contribution of the EOS is positive and similar in magnitude across models everywhere (because each model projects future warming seawater near Miami as well as Tokyo, Jakarta, Manila, Mumbai, Lagos, and New York City), the inter-model spread of how T100 variability will change is −4% to 20% for Miami ( Fig. 3b; 6% standard deviation) and even larger for most other locations (Supplementary Fig. 3; e.g., −9% to 40% spread with a 10% standard deviation for New York City). For near these seven example cities ( Fig. 3b and Supplementary  Fig. 3), the EOS always contributes to increasing sea level variability, which is in contrast to the CMIP5 projection of decreasing T100 variability for three of the cities (Tokyo, Mumbai, and Lagos). Hence, the temperature variability contribution to future sea level variability is much more uncertain compared to the EOS effect, both across models and among the example locations. Yet, there are places where the future sea level variability decreases (Fig. 1c, d), which presumably must match where decreases in temperature variability are larger than the EOS contribution (e.g., near Mumbai and Lagos in the tropical Indian and Atlantic Oceans, respectively; Supplementary Fig. 3).
We globally computed the relative contributions of changes in the rate of thermal expansion (i.e., the EOS) and temperature variability (e.g., associated with the annual cycle or ENSO changes) to thermosteric sea level (Fig. 4). This analysis follows that for Miami and the other large coastal cities used as examples, except that we consider ocean temperature over the full water column (Methods). We also assessed the veracity of the assumption that thermosteric sea level variability is a good proxy for sea level variability by comparing the CMIP5 inferred annual range and interannual standard deviation with the direct model output of SSH (Supplementary Figs. 1, 2), as well as the observed monthly anomalies of thermosteric sea level and SSH (Supplementary Fig. 4). As has been shown previously 21,35,36 , thermosteric sea level and SSH correlate well, although amplitudes deviate in some regions (e.g., the North Atlantic for interannual variability). The mostly larger amplitude of thermosteric sea level variability compared to SSH in CMIP5 (Supplementary Figs. 1, 2) is indicative of other processes, such as salinity 37 , also contributing to sea level variability. Since the impact of salinity on ocean stratification tends to increase with latitude 38 , we have restricted the sea level analysis to between 60°S and 60°N. Here, the similarity between thermosteric sea level and SSH, and especially the future changes of each (Supplementary Figs. 1, 2), gives confidence for our approach to determining the components of sea level variability changes.
The future sea level variability changes associated with the EOS (i.e., holding constant ocean temperature variability) are shown in Fig. 4a, b (annual range and interannual standard deviation, respectively). As expected for future near-global warming, the change in variability is positive almost everywhere for all models (parts of the equatorial Pacific and North Atlantic are exceptions, which we discuss in the next section). The patterns of variability changes are primarily determined by two factors that are illustrated in Fig. 3a: (i) future warming amount (for greater warming, there is a larger difference in the amount of thermal expansion per temperature fluctuation; i.e., the EOS slope increases with increasing temperature) and (ii) historical mean temperature (the rate that thermal expansion increases is less as seawater warms; i.e., the decreasing rate of slope change of the EOS with increasing temperature). Combining these factors, for both annual and interannual timescales of variability (Fig. 4a, b), the future sea level changes are largest at higher latitudes (>12%, except in parts of the North Atlantic and the South Pacific) and smallest in the tropics (0-4%). The greater changes in the higher latitudes compared to the tropics for variability associated with the EOS (e.g., New York City versus Jakarta; Supplementary  Fig. 3) are expected considering that, in the former region, the future warming is projected to extend to greater depths (factor i; Supplementary Figs. 5, 6) and historical mean temperatures are relatively cooler (factor ii).
Inter-model consensus and uncertainty. Unlike the almost certain likelihood that future ocean mean temperatures will be warmer nearly everywhere, for many regions there is inter-model uncertainty whether temperature variability will increase or decrease 5,39,40 . Considering the sea level response to changes in the annual range of ocean temperatures (Fig. 4c), there are only a few regions of multi-model agreement that the future variability will increase (29% area between 60°N/S, which is not stippled) and be of much larger magnitude than the change expected from the EOS (e.g., part of the tropical Pacific south of Hawaii). Around Miami, the effect on sea level of increasing temperature annual range (5%) is comparable to that expected from the EOS. For the interannual variability, inter-model consensus is also weak that future ocean temperature fluctuations will increase in most locations ( Fig. 4d; only 14% of the global area). Although, the tropical southwestern Pacific is a notable exception where more than 2/3 models agree on increased temperature-sea level variability, which is associated with future ENSO changes 23 . In general, where there is inter-model consensus of increasing ocean temperature annual and interannual variability (Fig. 4c, d), such as in the tropical Pacific, the CMIP5 projected change in SSH variability is especially large and consistent across models (Fig. 1c, d).
In a future warming environment, there is robust inter-model agreement of increasing sea level variability associated with the EOS that will almost certainly affect the annual cycle range and interannual standard deviation throughout most of the tropics and mid-latitudes (Fig. 4a, b and Supplementary Tables 1, 2). The only regions where some models project smaller future sea level variability related to the EOS are in parts of the equatorial Pacific and the high-latitude North Atlantic (stippling in Fig. 4a, b). Depth profiles of the ocean temperature variability and future warming for regions of reduced future sea level variability associated with the EOS reveal that some models do project future cooling at certain ocean levels, including where temperature variability is large ( Supplementary Figs. 5, 6). As the total change in thermal expansion rate is a vertical average weighted by the amount of temperature variability (Methods), the EOS component of sea level variability change can be negative (i.e., an opposite thermodynamic change to that illustrated in Fig. 3a) under such conditions even if the (unweighted) depth-averaged temperature increases. One such example is future cooling near the thermocline of the equatorial western Pacific, which increases the regional upper-ocean stratification ( Supplementary Fig. 7) more than if only near-surface warming were to occur 38 .
We have considered so far the components of future sea level variability changes under the inference that sea level variations are fully explained by the ocean temperature and density characteristics (i.e., sea levels determined solely by the EOS and temperature variability). It is informative to compare these results with the sea level variability changes projected by CMIP5; the multi-model mean of which (Fig. 1c, d) we discussed previously. Figure 5 relates, model-by-model, the CMIP5 output of SSH variability change (y-axes) to the component of change associated with either the EOS or temperature variability (x-axes). For both annual and interannual timescales, the multi-model mean changes associated with the EOS are closely aligned with the SSH multi-model mean changes in the tropics and mid-latitudes (Fig. 5a, b). Whereas the inter-model spreads of the temperature variability components explain much of the uncertainty of SSH projections (Fig. 5c, d), in a multi-model mean sense, less than half of the SSH changes are explained by temperature variability. In fact, the inter-model spread of the change of interannual temperature variability (Fig. 5d; e.g., −7.2% to 11.7% in the tropics; see Supplementary Tables 1, 2 for model-by-model statistics) suggests that this is the most uncertain of the sea level change components that we considered.
The proportion of increasing SSH variability associated with the EOS depends on the amount of future ocean heat uptake. For upper-ocean warming that is likely to occur in most of the tropics and mid-latitudes by the end of this century (e.g., regionalaverage T100 warms by 2°C; Fig. 6), the projected increases of the SSH annual cycles associated with the EOS are 4% (tropics; Fig. 6a) and 10% (mid-latitudes; Fig. 6b). Since we are assessing the projected changes with respect to the amount of future warming, these increases are relative to the latter part of the historical simulation (Methods), rather than the entire 20th century (e.g., as in Fig. 5). EOS contributions to increases of the SSH interannual variability, for 2°C warming, are similar (4% and 8% in the tropics and mid-latitudes, respectively; Fig. 6c, d). Even larger SSH annual cycle and interannual variability changes associated with the EOS are expected for greater warming (e.g., 2.5°C; Fig. 6). In contrast, the temperature variability changes are smaller and much more uncertain (especially interannual variability) for all future warming amounts expected anytime this century (Fig. 6). In fact, for the majority of models that simulate at least 2°C of T100 warming in the tropics and midlatitudes, there is no multi-model consensus that interannual temperature variability will increase; yet, the SSH variability increases in most of these same models (Fig. 6c, d).
Interpretation. Overall, for the CMIP5 climate models and RCP8.5 greenhouse warming scenario that we considered, it is perceivable that the annual cycle and interannual variability of sea level will both increase at many coastal locations. Such a tendency exists, at least in part, because the thermal expansion rate of seawater increases with warming (i.e., the EOS nonlinearity), which inherently causes density-related sea level variability to also increase even if the temperature variability itself does not change. Furthermore, the CMIP5 SSH projection is characterized by large inter-model spread in many regions (Fig. 1c, d) that is primarily related to uncertainty across models in how future ocean temperature variability will respond to continued greenhouse warming (Figs. 4c, d, 5c, d and 6).
Whereas our CMIP5 analysis supports the hypothesis that seasonal-to-interannual sea level variability will increase relative to changes in ocean temperature variability due to nonlinearity of the EOS, the EOS alone does not constrain future sea level variability to increase in a warming ocean. From this perspective, future sea level variability could also remain constant with reduced temperature variability. Indeed, CMIP5 models diverge substantially in the projected amount of sea level variability increase because of the uncertainty in ocean temperature variability changes (respectively, stippling in Figs. 1c, d and 4c, d; see also shading in Fig. 6). Certainly also contributing to future sea level and temperature variability changes in CMIP5 are changes in the variability of atmospheric forcing. Because such forcing impacts variability in the oceans locally as well as globally, the effects of atmospheric changes are difficult to separate from those of increased oceanic thermal expansion and stratification in coupled climate model simulations. For this reason, we consider how thermal expansion and stratification impact sea level and thermocline variability (a proxy for temperature variability) in an analytic, reduced-gravity ocean model prescribed with future warming but otherwise unchanged atmospheric forcing. We obtain analytic model solutions to periodic heat or wind forcings (Cases 1 and 2, respectively; Methods) in an idealized domain representing the subtropical gyre of the North Pacific during historical (cool) and future (warm) climate conditions. Overall, for the CMIP5 projected thermal expansion and ocean stratification increases in the subtropical North Pacific (Fig. 4a and Supplementary Fig. 7, respectively), sea level variability tends to increase, and thermocline variability to decrease, in spite of no change in atmospheric forcing between the historical and future solutions (Fig. 7). For the heat-forcing experiment (i.e., Case 1; Fig. 7a), sea level variability increases proportionally to the prescribed thermal expansion coefficient change, which   determines how much the surface-layer thickness varies in response to the unaltered variability of heating. For the windforcing experiment (i.e., Case 2; Fig. 7b), assuming that the forcing period is sufficiently long such that Rossby wave adjustments redistribute mass and heat throughout the basin (i.e., the circulation remains near equilibrium 41 ), sea level variability is insensitive to changing ocean stratification. (Note that increasing stratification is only partly due to the EOS effect on buoyancy, but also to the ocean warming and thus expanding faster near the surface than at depth 38 , e.g., Supplementary  Fig. 5.) Since Sverdrup transports are not affected by stratification 41 , and with no periodic heating to affect the surface-layer thickness, the sea level gradients in the gyre will only adjust to the wind forcing, which in this experiment vary the same in both climate conditions. Even though the periodic wind forcing is unchanged, thermocline variability is reduced with increasing stratification (Fig. 7b), because of decreasing thermocline slopes ( Supplementary Fig. 8), implying reduced redistribution of heat (i.e., less temperature variability) in the future. Larger sea level variability due to increased thermal expansion, combined with unchanged sea level variability due to increased ocean stratification (despite reduced temperature variability), in the analytic model (Fig. 7a, b and Supplementary Fig. 8) suggests that the EOS will on average contribute to increasing sea level variability in the future, both relative to any change in temperature variability and in absolute terms. The same thermodynamic constraints imposed by the EOS presumably hold in a cold climate as well (i.e., decreased thermal expansion and stratification; Fig. 7a, b). Consequently, reduced thermosteric sea level variability and enhanced thermocline variability would be expected during glacial periods.
Although future thermocline variability is reduced in the analytic model when ocean stratification increases (Fig. 7b), there is no wide-spread projected decrease of temperature variability in CMIP5 in the subtropical North Pacific or elsewhere (Fig. 4c, d). This is suggestive of stronger atmospheric forcing in the future (e.g., associated with increasing ENSO variability [18][19][20] balancing such tendencies on the thermocline, as well as contributing to projected sea level variability increases beyond that explained by increasing thermal expansion (Fig. 7a).

Discussion
In summary, the thermal expansion property of seawater, which is described by the nonlinear characteristic of the EOS, contributes to a global-scale tendency for the sea level to become more variable relative to temperature fluctuations with continued greenhouse warming. Such a tendency, which would alter the risks of coastal flooding and erosion beyond changes associated with anticipated accelerating sea level rise 10 , exists regardless of whether or not the ocean temperature variability increases in the future. An inference from these results is that the inter-model uncertainty of CMIP5 projected changes in SSH variability is explained in most regions by the uncertainty of future ocean temperature variability. A way forward to better describing coastal risks related to sea level variability is therefore either to reduce the uncertainty of how ocean temperature variability will change or, if that is not possible, to make assumptions based on the tendency for increasing sea level variability caused by the nonlinear thermal expansion of seawater.

Methods
Observations and reanalysis products. To describe the observed SSH and threedimensional ocean temperature, we used the ECMWF Ocean Reanalysis-System 5 (ORA-S5; ref. 42 ). We performed analyses on a globally uniform 1°latitude × 1°l ongitude grid between 60°S-60°N to encompass a large portion of the world's oceans, yet limiting the area covered by sea ice. During the 1993-2018 period that we analyzed, ORA-S5 assimilated satellite altimetry measurements of SSH along with observed depth profiles of ocean temperature and salinity. We compared the SSH from ORA-S5 for an example region around Miami, Florida with the nearby Virginia Key shore-based tide gauge record (1994-2018; Fig. 2a), which we acquired from the Quality Assessment of Sea Level Data archive 43 .
For all variables, we calculated the mean annual cycle and monthly mean anomalies with respect to the observed period. We also subtracted the locationspecific linear trend for the same period; thus the contribution of recent sea level rise is removed from our assessment of sea level variability (i.e., annual cycle range and interannual standard deviation). For analyses of interannual variability, we lastly high-pass filtered the monthly anomalies to remove any oscillations with periods longer than 11 years. This final step is necessary to distinguish the effect of accelerating ocean heat uptake associated with greenhouse warming from the interannual variability, which is especially important when considering future climate projections on centennial timescales.  Table 1. We assessed one experiment from each model, covering the period 1906-2005 using historical anthropogenic and natural forcings and then the future emission scenario (RCP8.5) for 2006-2100, which ignores volcanic and other natural aerosols. For each model, we first interpolated the dynamic SSH and three-dimensional ocean temperature to the 1°latitude × 1°longitude grid using bilinear interpolation. We calculated the mean annual cycle for the historical and future periods and then monthly mean anomalies with respect to each period. Following ref. 23 , we derived changes in the variability of SSH and ocean temperature by comparing the first 95 years (1906-2000, historical period) to the later 95 years (2006-2100, future period); thus, there was a large ratio between the climate change signal and any higher-frequency variability internal to the models. To assess how the ocean condition varies as a function of the future warming amount (Fig. 6), we compared 30-year running windows starting every 5 years from 1906-2100, with the 30-year climatology from the end of the historical experiment .
In addition to the multi-model mean, we considered the inter-model uncertainty of the projected changes by indicating regions where less than 2/3 of models (i.e., 19 out of 29 models) agree on the future change sign (stippling in Figs. 1c, d and 4c, d) or showing the inter-model spread (i.e., ±1 standard deviation; shading in Fig. 6). We also assessed the projections model-by-model for specific locations, such as near Miami and six of the largest coastal cities in the world, and globally (Figs. 3 and 5; Supplementary Tables 1 and 2; Supplementary Fig. 3). Furthermore, for each model and region (tropics and mid-latitudes), we calculated the percentage area where the sign of projected change disagrees from the multimodel average (Supplementary Tables 1 and 2; i.e., listing the area of decreasing variability).
Empirical EOS to describe seawater density. Seawater density is fully related to its state of temperature, salinity, and pressure (i.e., the EOS). From combining the First and Second Laws of Thermodynamics into the Gibbs function (G), volume (v), and therefore density (ρ), is given by the pressure (p) derivative of G that describes the molecular structure of the substance 45 (in this case, seawater): Since no analytical solution to the EOS exists, empirical expressions based on the measured properties of seawater are used to approximate it in various polynomial forms. We evaluate the density property of seawater using the International Thermodynamic Equation of Seawater-2010 (TEOS-10; ref. 33 ) and, in particular, the Gibbs Sea Water (GSW) Oceanographic Toolbox of TEOS-10, which uses an expression with 75 coefficients to calculate in situ density as a function of temperature, salinity, and sea pressure 46 . For our calculation of seawater density shown in Fig. 3a, we simplify the EOS by using a reference value for salinity (S, 35.2 g kg −1 ) and considering the sea pressure at 50-m depth (60.4647 dbar). The sea level nonlinearity 4 that we explore is explained by the temperature dependence of the thermal expansion coefficient (α), which is the temperature (T) derivative of the EOS with respect to the independent state variables (salinity and pressure) that we hold constant 45 : Since the thermal expansion coefficient increases with temperature, the change in seawater density and, hence, sea level also increases with temperature. This nonlinear behavior of seawater density is indicated by the changing slope of the EOS (Fig. 3a; i.e., the illustrated response of density to variations in T100). If instead the thermal expansion coefficient did not increase with temperature (e.g., under greenhouse warming conditions), then there would be no acceleration of sea level rise through thermosteric processes 4 . There would also be no increase of sea level variability without a change in temperature variability (i.e., the linearized EOS illustration in Fig. 3a) or some other mechanism that affects sea level (e.g., changes in seawater mass). We note that the other derivatives of the EOS (i.e., the coefficients of haline contraction and isothermal compressibility; ref. 45 ), are not considered in our calculations of seawater density since we use constant reference values of salinity and pressure.
Inference calculation of the components of sea level variability change. The sea level annual cycle and interannual variability are largely explained by ocean temperature variability (i.e., thermosteric sea level, H T ) and the resulting seawater buoyancy changes 11,12 . We therefor infer the thermosteric sea level variability (δH T ) based on the ocean density response to temperature changes, which we calculate for each latitude, longitude, and depth between 60°S-60°N using the rho_mwjf function 47 of the NCAR Command Language. Similar to our use of the empirical EOS to illustrate the seawater density response to changes in T100 around Miami, Florida (Fig. 3a) and the other example cities ( Supplementary  Fig. 3), the effect of changing seawater salinity is not considered. Although in our assessments of the full ocean (e.g., Fig. 4a, b), we do adjust for changing depth and pressure (i.e., isothermal compressibility) unlike in (2). Throughout, the effect of mass changes on sea level variability are not considered.
Since changes in the curvature of the EOS are small for ocean temperature fluctuations that are typically observed on seasonal-to-interannual timescales, we approximate the thermosteric sea level variability (e.g., annual cycle range or interannual standard deviation) as where α is the thermal expansion coefficient and δT is the temperature variability (i.e., range or standard deviation) that are each calculated at all ocean depths (z) and integrated from the ocean bottom (−D) to the surface. As in (2), α is representative of the time-mean temperature (i.e., climatology). For future changes of the mean temperature (hence also α) and δT that are small relative to their historical values, we decompose the future change (denoted by Δ) in the thermosteric sea level variability using a linear expansion of (3): On the right-hand side of (4), the first term represents the future change of α that is related to the nonlinearity of the EOS (e.g., Fig. 4a, b), and the second term is due to the future change in temperature variability (e.g., Fig. 4c, d). Specifically, ΔδH T is the future change relative to the historical condition of thermosteric sea level variability (either the annual cycle range or interannual standard deviation; Supplementary Figs. 1, 2) that we infer by combining the terms on the right-hand side of (4): α historical and Δα are respectively the historical mean and future change of the thermal expansion coefficient, and the temporal distinction is similar for temperature variability (δT historical and ΔδT).
Through comparisons of the observed SSH (i.e., ORA-S5) and inferred sea level (i.e., Eq. 4) characteristics, using calculations of the anomaly correlation coefficients at each grid location ( Supplementary Fig. 4), we determined that the inference method resolves most of the annual cycle and interannual variability (in 83% and 77% of the area, respectively). The comparison is especially close for the annual cycle range almost globally and the interannual standard deviation in the tropical Pacific. The patterns of root mean square errors mostly mirrors the correlations ( Supplementary Fig. 4).
As we noted, the inferred sea level variability that we calculate is larger than the direct CMIP5 output of SSH almost everywhere ( Supplementary Figs. 1, 2 show maps of the comparison of annual cycle range and interannual standard deviation, respectively, for the historical and future periods). Such a result is to be expected, as there are contributions to sea level variability other than temperature (e.g., salinity) that may compensate the effect of temperature variability on the ocean density 37 . Yet, comparing the future change of sea level variability from CMIP5 and using the inference method shows similar patterns for both the annual cycle and interannual standard deviation (Supplementary Figs. 1, 2). Since we consider only future changes in our decomposition of the sea level variability associated with the EOS and temperature variability (Figs. 4-6, Supplementary Tables 1, 2), the robustness the total future change patterns gives confidence in the method. Analytical determination of the sea level variability response to future warming. To illustrate sensitivities of the ocean response to changing mean temperature when the atmospheric forcing (i.e., variability of heat or momentum fluxes) is not changed, we present two solutions to a reduced-gravity ocean model 48,49 . In each case, the analytic model is driven by periodic surface fluxes. Specifically, we analyze how the SSH and thermocline variability (a proxy for temperature variability) each change as a function of changes in the mean thermal expansion coefficient and ocean stratification. Solutions are obtained in an idealized oceanic domain resembling the subtropical gyre in the North Pacific (Supplementary Fig. 8); however, the model configuration and experiments are generic enough to make generalizations from the results about the broader oceanic response to future warming (or past cooling) of the climate.
Our analytic model has one active layer with uniformly varying density ρ 1 at the surface, and an infinitely thick layer below with constant density ρ 2 = 1,030 kg m −3 .
Solutions are obtained in a closed ocean basin extending from latitude y s = 15°N in the south to y n = 35°N in the north. Longitudinally, the domain extends from 140°E to 120°W, which, respectively, are the western (x w ) and eastern boundaries (x e ). The model is directly forced by Ekman pumping of the form w ek ¼ W ek sinð y À y n y s À y n Þ, where the amplitude W ek may depend on time. In steady state, the interior ocean circulation is in Sverdrup balance 48,49 , and the surface-layer thickness is described by where f is the Coriolis frequency, β its meridional derivative, and g 0 ¼ gðρ 2 À ρ 1 Þ=ρ 2 is the reduced gravity with the gravitational acceleration being g = 9.81 m 2 s −1 . The eastern boundary surface-layer thickness h e is determined by requiring that the total mass in layer 1, M ¼ R R ρ 1 h dx dy = 1028 kg m −3 × 200 m × area of the domain, is conserved. Note that infinitesimally thin boundary layers are required to close the circulation, which are assumed to have no further impact on the interior ocean solution 49 . The SSH (η) is described by relative to a constant level η o . Note that the model, and the sea level Eq. (6) in particular, imply that barotropic adjustments, which eliminate horizontal pressure gradients in the deep layer, occur instantaneously. Analytical solutions are compared for a historical and future ocean, with the difference being that the future ocean is more stratified than the historical one due to warming from above and an increased thermal expansion coefficient, which enhances the buoyancy of the upper layer. Specifically, we increase the mean difference between deep-and surface-layer densities (i.e., ρ 2 À ρ 1 ) from 2 kg m −3 in historical solutions by a range from −10% to 20% in the future simulations, with 13.3% being the average CMIP5 stratification increase in the subtropical North Pacific region that we are modeling ( Supplementary Fig. 8). Stratification changes in CMIP5 models are calculated following the methodology of ref. 38 as the density difference between 200 m and the surface (Supplementary Fig. 7).
Case 1, oscillating heat forcing: first, we consider a solution without wind forcing (W ek ¼ 0) and a heat flux that is uniform over the domain, but oscillating in time t. Assuming constant salinity, the EOS for the surface density is ρ 1 ¼ ρ 1 À αT 0 , where ρ 1 is the mean density, α is the thermal expansion coefficient appropriate at the mean temperature (i.e., a linearization is performed), and T′ is the upper-layer temperature anomaly (e.g., associated with the annual cycle).
It follows directly from (5) that the thermocline depth remains constant across the domain at all times in this solution (i.e., h = h e ), as shown in Fig. 7a. According to the linearized EOS and (6), the SSH anomaly associated with T′ is given by Since temperature anomalies are zero initially, the uniform heating that we prescribe ensures that T′ is also spatially uniform at all times. Assuming the same thermocline depth (h e ) and heat flux (i.e., T′) in historical and future solutions, and neglecting changes in deep-ocean density (ρ 2 ), it follows from (7) that changes in η′ (e.g., the annual cycle range of sea level) are equal to changes in the α (i.e., the future change in sea level variability and the thermal expansion coefficient are equal). Figure 7a shows results for thermal expansion changes ranging from −10% to 20%, with 4.0% being the CMIP5 multi-model average increase in the subtropical North Pacific (Fig. 4a; domain average between 15°N-35°N, 120°E-120°W). Case 2, oscillating wind forcing: densities remain constant in the second set of solutions (i.e., no oscillating heat forcing), however, the wind forcing is modulated in time according to W ek ¼ W o þ W 0 cos t, where W o and W′ are the curl of the mean wind and the y-dependent amplitude of the time-varying wind, respectively. The frequency of the forcing (ω) is assumed to be sufficiently slow, such that the circulation remains near the Sverdrup balance; i.e., (5) holds at all times. The only effect of thermal expansion coefficient changes is that associated with the stratification change in g′. Since g′ occurs in the denominator of the forcing term in (5), it follows that for the same W ek , increased stratification is associated with a smaller slope of the future thermocline as well as reduced thermocline variability. Supplementary Fig. 8 shows, for the historical and future climate conditions, the thermocline depth (i.e., surface-layer thickness, h) and SSH solutions when the wind forcing reaches its maximum strength (at t = 0) and the difference between that time and when minimum strength of the wind forcing occurs (at t ¼ π=ω). We set the wind forcing as follows: W o ¼ 0:5 10 À6 m s À1 and W 0 ¼ 0:1 W o . In all solutions, wind forcing causes the thermocline to be deep in the west and shallow in the east. In the future, with increased stratification and the associated reduction of the thermocline slope, and since the domain-integrated mass is conserved, the thermocline deepens in the east and shoals in the west. This weaker wind-forced mean response of the future thermocline (i.e., smaller slope) also translates to a reduced thermocline depth range (δh) in the future relative to the historical solution (i.e., Δδh < 0). Despite reduced thermocline variability, the historical and future SSH variability (i.e., range between strong and weak wind forcing) are nearly identical. This follows from the Sverdrup transports and mean thermocline depth being unchanged between historical and future solutions. Figure 7b shows the domain-average relative changes between the future and historical solutions for the variability of thermocline depth and SSH (i.e., the effect of CMIP5 projected stratification increase), as well as for a range of other stratification changes (−10% to 20%) perceivable in past or future climates. For thermocline depth, the change is defined as 100 Δ ∂h j j= δh j j historical , which is the absolute value of the surface-layer thickness difference between strong and weak wind forcing (i.e., the thickness range) for the future minus historical solution, relative to the historical variability of such ( Supplementary Fig. 8). The relative change for SSH is defined analogously, with surface-layer thickness replaced by SSH.

Data availability
The data that support the findings of this study are available from the Asia-Pacific Data-