Evolving patterns of sterodynamic sea-level rise under mitigation scenarios and insights from linear system theory

Long-term behaviour of sea-level rise is an important factor in assessing the impact of climate change on multi-century timescales. Under the stabilisation scenario RCP4.5, Sterodynamic Sea-Level (SdynSL) and ocean density change in the CMIP5 models exhibit distinct patterns over the periods before and after Radiative Forcing (RF) stabilisation (2000–2070 vs. 2100–2300). The stabilisation pattern is more geographically uniform and involves deeper penetration of density change than the transient pattern. In RCP2.6, 4.5 and 8.5, the spatiotemporal evolution of SdynSL change can be approximated as a linear combination of the transient and stabilisation patterns. Specifically, SdynSL change is dominated by the transient pattern when RF increases rapidly, but it is increasingly affected by the stabilisation pattern once RF starts to stabilise. The growth of the stabilisation pattern could persist for centuries after RF ceases increasing. The evolving patterns of SdynSL change can also be approximated as a linear system's responses (characterised by its Green’s function) to time-dependent boundary conditions. By examining SdynSL change simulated in linear system models with different estimates of Green's functions, we find that both the climatological ocean circulation and the ocean's dynamical response to RF play a role in shaping the patterns of SdynSL change. The linear system model is more accurate than the univariate pattern scaling in emulating the CMIP5 SdynSL change beyond 2100. The emergence of the stabilisation pattern leads to a 1–10% decrease in the ocean's expansion efficiency of heat over 2000–2300 in RCP2.6 and 4.5.


Introduction
Sea-level rise is a major consequence of anthropogenic climate change and poses severe threats to coastal communities. Satellite altimetry reveals more than 7 cm of global mean sea-level rise since 1993 at a rate about 3 mm year −1 (Church and White 2011;Nerem et al. 2018). Thermal expansion of the ocean (termed the thermosteric component) accounts for about one third of that. The rest (termed the barystatic component) is due to transfer of mass to the ocean from ice sheets, glaciers and terrestrial water storage (Church et al. 2013;Chambers et al. 2017). Models contributing to the Coupled Model Intercomparison Project phase 5 (CMIP5; Taylor et al. 2012) have been used to project sea-level rise under various Representative Concentration Pathway (RCP) scenarios (van Vuuren et al. 2011). The resulting likely global mean sea-level change ranges from 0.26-0.55 m in RCP2.6 (strong mitigation scenario) to 0.45-0.82 m in RCP8.5 (high emission scenario) during 2081-2100 relative to 1986-2005, and with ongoing sealevel rise over subsequent centuries (Church et al. 2013).
The ocean responds to climate change on timescales from decades to millennia. For a doubling of CO 2 concentration in the atmosphere, most of the upper ocean (depth < 1 km) reaches its full response within several hundred years, while it takes thousands of years for the deep ocean and the Southern Ocean to reach their full response (Stouffer 2004). Stabilising Radiative Forcing (RF) therefore has little effect on stabilising sea level within decades or centuries (Meehl et al. 2005). In RCP4.5 for example, Global Mean Thermosteric Sea-level (GMTSL) rise continues long after RF is 1 3 stabilised at 2100, making sea-level rise by 2300 about three times more than that by 2100 (0.14-0.23 m vs. 0.38-0.66 m; Church et al. 2013). Li et al. (2013) showed that GMTSL rise could continue for several thousand years after RF stabilisation, due to slow adjustments in the deep ocean. Mass loss from ice sheets and glaciers also continues after RF stabilisation, giving an ongoing barystatic contribution to sea-level rise. For 2000-2300, its magnitude could be as much as twice the GMTSL contribution as suggested by a simple climate model (Nauels et al. 2017).
Sea-level rise is not spatially uniform. In the altimetry record, regional trends of sea-level rise can deviate from the global mean trend by a factor of 2-3 over two decades (e.g. Merrifield 2011;Zhang and Church 2012;Fasullo and Nerem 2018). These regional deviations are strongly associated with internal variability of the climate system (referred to as unforced variability) and do not necessarily reflect a response to long-term climate change (Lyu et al. 2014). Statistical methods such as regression and principal component analysis can be used to reduce the impact of unforced variability on regional sea-level trends seen in the altimetry record (e.g. Zhang and Church 2012;Hamlington et al. 2014Hamlington et al. , 2019. Sea-level rise in response to anthropogenic climate change (termed the forced response) exhibits regional deviations due to various physical processes (Slangen et al. 2014. A major part of it, called Sterodynamic Sea Level (SdynSL) change, arises from changes in density and circulation of the ocean . SdynSL change has both global mean and regional signals; the former is GMTSL change, and the latter is termed Dynamic Sea Level (DSL) change. SdynSL change can be alternatively partitioned into steric sea level change, due to local ocean density change, and regional deviations of manometric sea level change with respect to its global mean, due to redistribution of ocean mass (Gill and Niller 1973;Gregory et al. 2019). Steric sea level change is a good approximation of DSL change in regions away from the continental shelf and without strong barotropic currents (Lowe and Gregory 2006).
Atmosphere-Ocean General Circulation Models (AOGCMs) are often used to study the forced response of SdynSL change under prescribed greenhouse gas forcing. In the twenty-first century projections, AOGCMs simulate a smaller sea-level rise than the global mean in the Southern Ocean (south of 50° S) and a greater sea-level rise than the global mean in the Arctic. Regionally AOGCMs produce DSL change with complex ridge-and-trough structures in middle and high latitudes (e.g. Yin et al. 2010;Yin 2012;Church et al. 2013;Lyu et al. 2020). Some of these features are associated with changes of oceanic gyres and boundary currents (Zhang et al. , 2017Terada and Minobe 2018). The spatial Standard Deviation (STD) of DSL change is roughly 15-30% of GMTSL change by the end of the twenty-first century in Gregory et al. (2001). Bilbao et al. (2015) argued that the twenty-first century DSL change in the CMIP5 models has a nearly constant pattern with its magnitude related to the GMTSL change. In addition to DSL change, regional sea-level change is also affected by gravitational, rotational and solid-Earth deformational effects . These latter processes are not considered in this study.
Understanding long-term behaviour of sea-level rise is important in assessing the impact of anthropogenic climate change on multi-century timescales. Efforts have been made to study the long-term response of global mean sea-level rise under various mitigation scenarios (e.g. Schaeffer et al. 2012;Mengel et al. 2018), and how such information can be used to guide climate policy (e.g. Clark et al. 2016Clark et al. , 2018Nauels et al. 2019). However, much less is known about how regional sea level will change under mitigation scenarios, as existing studies generally focus on its shorter-term transient response instead . What is the pattern of regional sea-level change after RF stabilisation? How is it similar to or different from the pattern in the twenty-first century projections? Will the regional differences in sealevel rise further increase after RF stabilisation? Questions like these remain unanswered, limiting our confidence in multi-century sea-level projections on regional scales.
Another motivation for studying this issue is the need to reduce the computational cost of multi-century sea-level projections, so that its impact over a wide range of emission scenarios can be assessed. To address this need, simple climate models for quantities averaged over the globe or hemispheres have been developed (e.g. Meinshausen et al. 2011a;Palmer et al. 2018), and their outputs used to project DSL change by scaling a constant pattern derived from the CMIP5 simulations (e.g. Perrette et al. 2013;Palmer et al. 2020). The scaling method assumes that global and local responses of a system are linearly related in time. That however is not necessarily the case when the response unfolds over a range of timescales. For example, after RF stabilisation one would expect an emergence of the long-term response, which gradually replaces the shorter-term response that dominates during the increase of RF. An alternative and more general approach is offered by linear system theory, which in our application assumes that sea-level change at any given time is the sum of responses to all the forcings that occurred in the past. It makes no assumption about how the response to any forcing develops in time, except that the responses are independent from one another, which justifies their being summed. This theory forms a simple conceptual model for understanding sea-level change under any given emission scenarios.
Linear system theory has been applied in several climate-related studies. For example, Good et al. (2011Good et al. ( , 2013Good et al. ( , 2015 showed that abrupt CO 2 experiments can be used to reconstruct and understand climate responses in the 1% per year CO 2 increase experiment and RCP scenarios. Long et al. (2020) decomposed the RCP2.6 response into a rescaled RCP4.5 response and a response to a RF decrease. Marčelja (2010) showed that global ocean heat content change and GMTSL change in the CMIP5 models can be approximated using a linear system model based on a single column model of the upwelling-diffusion process. Nonetheless, to what extent linear system theory applies to regional sea-level change in climate models remain unclear.
In this paper, we study the time-evolving patterns of SdynSL change under mitigation scenarios (in which RF stabilises or decreases after an initial increase) and investigate to what extent they can be approximated as a linear system's response to changing RF. Details of the CMIP5 model outputs and how they are processed are included in Sect. 2. We first demonstrate the transient and stabilisation patterns of SdynSL change from RCP4.5 in Sect. 3. In Sect. 4, we apply the RCP4.5 patterns to explain SdynSL change in RCP2.6, 4.5 and 8.5 via a two-pattern regression model, and demonstrate a method to construct a two-pattern scaling model for DSL projections. In Sect. 5, ocean density change associated with the SdynSL change is investigated. In Sect. 6, we introduce the method of emulating SdynSL change using linear system models based on the Green's function of AOGCMs. The following two Sects. (7 and 8) then compare the CMIP5 SdynSL change with that from two configurations of the linear system models. Discussion and conclusion are presented in Sects. 9 and 10, respectively. In order to gain a comprehensive view of sea-level rise, both global mean and regional deviations of SdynSL rise are studied. We focus more on the latter aspect as it has seen less progress in recent years. Note that linearity has different meanings in linear system models and pattern scaling. For linear system models the "linearity" means that the response to a sum of forcings equals the sum of responses to individual forcings, while for pattern scaling it means that regional changes are linearly related to the global mean change.

Data and variables
The CMIP5 model outputs under Representative Concentration Pathway (RCP) scenarios and their extensions (to 2300) are used in this study. The RCP scenarios are a common set of scenarios for assessing climate change with climate models (Meinshausen et al. 2011b;van Vuuren et al. 2011). The three emission scenarios used here are RCP2.6, 4.5 and 8.5. Their RF profiles are shown in Fig. 1a. In RCP4.5, RF increases to 4.5 W m −2 near 2100 and then stabilises. In RCP2.6, RF peaks at 3 W m −2 near 2040 and then declines to 2.6 W m −2 by 2100 and to 1.7 W m −2 by 2300. In RCP8.5, RF keeps increasing throughout the twenty-first and twenty-second centuries (at a decreasing rate after 2050), and gradually stabilises by 2250 at about 12.5 W m −2 . Among the three RCPs, RCP4.5 is most suitable for investigating the ocean's response to a stabilisation of RF as it features a 200 year period (2100-2300) with a steady RF. We therefore focus on RCP4.5 when analysing sea-level change of different RF profiles in Sects. 3, 5 and 7.
The values shown in Fig. 1 are from integrated assessment models, and the exact RF might vary in different CMIP5 models. In addition, RF in the CMIP5 models also has a complex and likely time-variable spatial pattern, which is not explicitly considered in this study; we are concerned only with the time-evolution of the global mean RF.
Eleven CMIP5 models that have run to 2300under RCPs (2006-2100 and their extensions (2100-2300) are used in this study (Table 1). Their sea surface height "zos", ocean temperature "thetao" and ocean salinity "so" are yearly averaged and interpolated onto a 2 × 2° longitude-latitude grid (using the ocean mask shown in Fig. S1) before further processing. Model drifts are removed by subtracting a 3rd order polynomial fit of the parallel control run at each grid point (Sen Gupta et al. 2013). After de-drifting, linear trends of DSL change during 1860-2300 of the piControl experiment are within ± 0.5 cm per 100 years in most regions Fig. 1 Time evolution of total anthropogenic and natural Radiative Forcing (RF) relative to 1765. a Global mean RF in RCP2.6, 4.5 and 8.5. b Rates of global mean RF change in the three RCPs. Data are available from http:// www. pik-postd am. de/ ~mmalte/ rcps/ for individual model (drifts would be further reduced in the multi-model mean). The cubic fit was used in Lyu et al. (2014) and Zhang et al. (2017), while a linear or quadratic fit was used by other studies (e.g. Meyssignac et al. 2017;Lyu et al. 2020). For four models that do not have control runs extending to 2300, linear fits of available control runs are first extrapolated to 2300 and then subtracted.
DSL is calculated from sea surface height "zos" with its global mean removed at each year. Steric Sea Level (SSL) Anomaly (SSLA) is calculated from in-situ density anomalies ′ following: where ′ is relative to the 1861-1880 level of the historical experiment, 0 = 1029 kg m −3 is a constant reference density, and H is the local ocean depth. In-situ density is calculated from ocean potential temperature "thetao" and ocean salinity "so" using the 1980 Equation of State. GMTSL change is calculated as the global average of SSLA using ocean salinity averaged over 1861-1880. Calculating GMTSL change from density change is also to address the issue that most of the CMIP5 models are based on the Boussinesq approximation and hence conserve volume (rather than mass). The sum of DSL change and GMTSL change gives SdynSL change. The sum of SSLA and manometric sea level change gives DSL change (after the global mean is removed). A 20-year running mean is applied to all diagnostics used in this study to reduce the unforced variability.
Throughout this paper, the "transient" and "stabilisation" periods are referred to when RF is changing rapidly (e.g. 2006(e.g. -2070Fig. 1b) and when RF is staying fairly constant (e.g. 2100-2300 in RCP4.5), respectively. We also use transient and stabilisation patterns to denote changes of a quantity during the transient and stabilisation periods, respectively. Note that the nomenclature only refers to the temporal behaviour of the RF, not that of ocean states.

Transient and stabilisation patterns of SdynSL change in RCP4.5
We start by demonstrating that SdynSL change has different patterns for the transient and stabilisation periods of RCP4.5. These patterns are derived as changes of SdynSL over the corresponding periods. We focus on the multi-model mean results in the following because the forced response is our main concern here. The multi-model mean patterns shown here are robust features if the ensemble size ≥ 5. For convenience, we use "Δ" to denote a change of a quantity over a time period, and ΔDSL to describe how local ΔSdynSL deviates from the global mean change hereafter. All changes are calculated as difference between the last and first 20 year averages over the given time interval. Using the last and first 10 year averages or calculating the leading empirical orthogonal function lead to very similar results. Due to a lack of deep water formation in the North Pacific, zonal averages of ΔSdynSL are calculated for the Indo-Pacific and Atlantic Oceans separately. The same basin average is applied to Sects. 5 and 7. Definitions of ocean basins are shown in Fig. S1. The transient pattern of ΔSdynSL exhibits rich regional structures, for example meridional dipoles in the North Pacific and the North Atlantic and zonal bands in the Southern Ocean (Fig. 2a). ΔDSL is 10-50% of ΔGMTSL (about 10 cm; Fig. 2e) in some subtropical and subpolar Cross marks indicate models from which sea surface height "zos", ocean temperature "thetao" and ocean salinity "so" are available under given experiments. The time periods during which the pre-industrial Control (piControl)  oceans and greater than 50% of ΔGMTSL in the North Atlantic, the Arctic and the Southern Ocean ( Fig. 2a black and red contours). In the tropics, ΔDSL is less than 10% of ΔGMTSL. The spread of ΔSdynSL across the CMIP5 models is lower in the tropics than in higher latitudes, and greater in the Atlantic Ocean than in the Indo-Pacific Ocean (Fig. 2c, d black lines). These features are consistent with previous studies of ΔDSL in the twenty-first century projection (Yin 2012;Bilbao et al. 2015;Lyu et al. 2020). The stabilisation pattern of ΔSdynSL is different from the transient pattern in both structure and magnitude (Fig. 2b). In the North Atlantic, the meridional dipole seen in the transient pattern is replaced by a negative ΔDSL of about 5 cm (compare red and black lines in Fig. 2c). In contrast, a positive ΔDSL of 1-3 cm develops in most of the Indo-Pacific  (Fig. 2d red lines). The contrast between the Pacific and the Atlantic was also found in Krasting et al. (2016) but with the opposite sign during high carbon emission rates (related to a rapid increase of RF). They suggested that a difference in ventilation timescale controls the difference in sea-level response between the two basins. In the Arctic, ΔDSL reverses sign from positive in the transient pattern to negative in the stabilisation pattern (compare red and black lines in Fig. 2c), consistent with the change in the North Atlantic (see also Lyu et al. 2020). The pattern of ΔSdynSL is least modified in the Southern Ocean after RF stabilisation (Fig. 2b). The zonal bands in the transient pattern are still visible in the stabilisation pattern, but become less coherent spatially. Regional deviations of ΔSdynSL (i.e. ΔDSL) during the stabilisation period (2100-2300) are mostly less than 10% of ΔGMTSL (Fig. 2b), suggesting SdynSL rise becomes more uniform after RF stabilisation. Even in the Southern Ocean, where the pattern of ΔSdynSL is relatively stable, local ΔSdynSL relative to the global mean is greatly weakened after RF stabilisation, indicated by the shrinking of red contours from Fig. 2a to Fig. 2b. Consequently, ΔDSL over 2000-2300 is similar to 2000-2070, with a smaller contribution from 2100 to 2300 (compare blue and black lines in Fig. 2c, d).
The temporal evolution of ΔDSL exhibits a slowdown around 2070 and a stabilisation after that to the north of 45° S. This is demonstrated in Fig. 2f where the spatial STD is used as a global measure of ΔDSL (relative to [1986][1987][1988][1989][1990][1991][1992][1993][1994][1995][1996][1997][1998][1999][2000][2001][2002][2003][2004][2005]. The slowdown is less pronounced for the global ocean, because DSL continues to change in the Southern Ocean (compare red and black lines in Fig. 2f). The stabilisation of ΔDSL can be identified in all eleven models, but their patterns stabilise at different magnitudes. This is caused by differences in models' forced responses as well as their unforced variability. Because the unforced variability of DSL is reduced in the multi-model mean, the STD of multimodel mean ΔDSL tends to be lower than that from individual models (Fig. 2f).

Formulation of the two-pattern regression model
The transient and stabilisation patterns of SdynSL change in RCP4.5 represent the short-and long-term responses of SdynSL change over a 300 year period (i.e. 2000-2300). In this section, we examine whether DSL change at any given time can be explained as a linear combination of those two patterns (after their global means are removed). This is equivalent to assuming that DSL change is governed by a system of two timescales, which can be formulated as a multi-variable linear regression model: 1 ( ) and 2 ( ) are column vectors of the transient and stabilisation patterns ( is the position vector) as in Fig. 2a, b. is a column vector of SdynSL change relative to 1986-2005 from the CMIP5 multi-model mean, and is a column vector of regression errors. Angle brackets denote deviations from the spatial mean, which converts SdynSL change to DSL change. x 1 and x 2 are regression coefficients for corresponding patterns, solved through the ordinary least-squares method for each individual year ( t = t 1 ⋯ t n ). Equation (2) is referred to as the two-pattern regression model hereafter.
We calculate x 1 and x 2 for from RCP4.5, the high-emission scenario RCP8.5 and the strong mitigation scenario RCP2.6. A similar method was used by Ceppi et al. (2018) to study patterns of sea surface temperature change.
The two-pattern regression model is applied to DSL change (denoted by the use of " < > " in Eq. (2)) so that x 1 and x 2 are solely determined by regional structures, which are our focus. If the two-pattern regression model were applied to SdynSL change (i.e. with global mean change included), the algorithm would aim to scale up 1 ( ) and 2 ( ) to fit the global mean change ( where "[]" denote the spatial mean) at the expense of worsening the fit of regional changes ( ⟨ ⟩ ), because is much larger than ⟨ ⟩ (Fig. 2c, d).
DSL change from the two-pattern regression model is denoted as ⟨̂ ⟩ and calculated as: Contributions from the transient and stabilisation patterns to ⟨̂ ⟩ can be quantified using the terms associated with x 1 and x 2 , respectively. To evaluate how well DSL change in the CMIP5 models can be explained by the two-pattern regression model, we focus on two global metrics: (1) spatial STD of DSL change and (2) Root Mean Square Error (RMSE) of DSL change. The STD is calculated for both ⟨ ⟩ and ⟨̂ ⟩ , while the RMSE is calculated for differences between ⟨ ⟩ and ⟨̂ ⟩ . One caveat of our regression model is that ⟨ 1 ⟩ and ⟨ 2 ⟩ are correlated at 0.7, which might hinder a clear attribution of ⟨ ⟩ to ⟨ 1 ⟩ and ⟨ 2 ⟩ . To test that, we repeat the regression by excluding the Southern Ocean (90° S-45° S), which reduces the correlation between ⟨ 1 ⟩ and ⟨ 2 ⟩ to -0.3. The resulting x 1 and x 2 are not significantly different from when the Southern Ocean is included (not shown). (2) We also compare DSL change from the univariate pattern scaling with that from the CMIP5 models. Following Bilbao et al. (2015) the univariate pattern scaling can be written as: ( ) is a column vector of the time series of GMTSL change. c(r) is the regression coefficient between the GMTSL change ( ) and DSL change ⟨ (r, )⟩ for a given location r . We first solve for c(r) based on ( ) and ⟨ (r, )⟩ from 2000-2100 in RCP4.5 using the ordinary (4) ( )c(r) = ⟨ (r, )⟩ + (r, ), r = r 1 … r n , least-squares method (Eq. (4)), then use the resulting c(r) to estimate DSL change for a given emission scenario from the GMTSL change in that scenario (Eq. (5)).

Skill of the two-pattern regression and behaviours of individual patterns
We start by regressing RCP4.5 DSL change against the transient and stabilisation patterns. Because the patterns are derived from RCP4.5, we expect that DSL change is mainly explained by an increase of the transient pattern during 2000-2070 and by an increase of the stabilisation pattern after that. The regression result is in excellent agreement In comparison, the univariate pattern scaling overestimates the STD of DSL change in CMIP5 after 2100 ( Fig. 3a red line), which leads to a RMSE of 5 cm by 2300 ( Fig. 3d red line). The overestimate is due to the emergence of the stabilisation pattern of SdynSL change, which is more uniform than its transient counterpart as discussed in Sect. 3. A major part of DSL change in RCP8.5 and 2.6 can also be approximated as a linear combination of the transient and stabilisation patterns of RCP4.5. The two-pattern regression explains 21 out of 25 cm and 3.0 out of 3.2 cm of the STD of DSL change by 2300 in RCP8.5 and 2.6, respectively (Fig. 3b, c black lines). In both RCPs, the behaviours of individual patterns are associated with their RF profiles. Specifically, the amplitude of the transient pattern tends to follow that of the RF. For example in RCP2.6, the transient pattern declines after 2070 ( Fig. 3i black line) following the decline of its RF (Fig. 1a). The amplitude of the stabilisation pattern is small while RF is increasing rapidly, but it starts to increase after RF starts to stabilise. For example in RCP8.5, the stabilisation pattern starts to grow at 2100 (Fig. 3h black line) after the increase of RF starts to slow down at 2050 (Fig. 1b). The RMSE from the two-pattern regression in RCP8.5 (2.5-10.0 cm) is higher than that in RCP4.5 (about 1.0 cm) over 2100-2300, and increases with time (Fig. 3e). This implies that additional patterns are required to better explain DSL change in RCP8.5. Similar to RCP4.5, the univariate pattern scaling overestimates the STD of DSL change in RCP2.6 after 2100 ( Fig. 3c red line), which leads to 4.0 cm of RMSE by 2300 compared to 1.5 cm in the twopattern regression (Fig. 3f). The overestimate is also evident in RCP8.5 but by a lesser extent, which leads to 15 cm of RMSE by 2300 in the univariate pattern scaling compared to 10 cm in the two-pattern regression (Fig. 3e).

A simple two-pattern scaling model
Given the two-pattern regression model is more accurate than the univariate pattern scaling in emulating the CMIP5 DSL change, a potential application is to convert it to a twopattern scaling model. Note that the two-pattern regression model itself cannot be used for DSL projections directly because the regression coefficients x 1 and x 2 are unknown for general emission scenarios. In the following, we demonstrate a simple two-pattern scaling model that is able to capture the two-timescale behaviour seen in the CMIP5 DSL change.
The two-timescale behaviour can be generated by a simple linear system model with its response to an impulse of RF prescribed as: where 1∕ 1 and 1∕ 2 are the e-folding timescales of the system, and a 1 (r) and a 2 (r) represent the patterns and magnitudes of the responses on each timescale. For a given RF profile, the evolution of the simple linear system model can be estimated by convolving G in Eq. (6) with RF: where t is elapsed time since a reference state at t 0 , lead time relative to t , and RF ′ is RF anomalies relative to RF at t 0 . s i (t) indicates the time-evolution of each response pattern, which corresponds to x 1 and x 2 in the two-pattern regression model (Eq. (3)). For a positive step change of RF and 1∕ 1 ≪ 1∕ 2 , s 2 would continue to rise long after s 1 has stabilised, consistent with the behaviours of x 1 and x 2 in the two-pattern regression model ( Fig. 3 bottom row). To make DSL projections, one can first solve for a 1 (r) and a 2 (r) based on DSL change and RF from the CMIP5 models using the ordinary least-squares method (Eqs. (7) and (8)), then use the resulting a 1 (r) and a 2 (r) to estimate DSL change for a given emission scenario from its RF profile (Eqs. (7) and (9)).
As an example, we calculate 1 ( ) and 2 ( ) by choosing 1∕ 1 = 10 years and 1∕ 2 = 1000 years, setting the year 2000 as t 0 , and solving for a 1 (r) and a 2 (r) based on DSL change and RF ′ (Fig. 1a) from 2000 to 2300 in RCP4.5. The resulting a 1 (r) and a 2 (r) are very similar to the transient and stabilisation patterns in Fig. 2a, b (not shown). Note that by setting t 0 to the year 2000, we neglect the ongoing DSL response to RF imposed before 2000. The two-pattern scaling model based on above estimates of parameters is also used to project DSL change in RCP2.6 and 8.5 from their RF profiles (Fig. 1a). These estimates of DSL change are compared with that from the CMIP5 multi-model mean to assess the performance of the simple two-pattern scaling model (Fig. 3).
The two-pattern scaling model captures the two-timescale behaviour of DSL change as seen in the two-pattern regression model in all three RCPs (Fig. 3 bottom row). It even has a lower RMSE (< 1 cm) than the two-pattern regression in explaining DSL change in RCP4.5 (Fig. 3d blue line). The performance of the two-pattern scaling model degrades (6) G(r, t) = a 1 (r)exp − 1 t + a 2 (r)exp − 2 t , when applying a 1 (r) and a 2 (r) derived from RCP4.5 to project DSL change in RCP2.6 and 8.5. In RCP2.6, the twopattern scaling underestimates the magnitude of the CMIP5 DSL change (Fig. 3c blue line), but it still features a lower RMSE than the univariate pattern scaling (2 cm vs. 4 cm by 2300, Fig. 3f blue line). This implies that the two-pattern scaling model is an alternative to the univariate pattern scaling in strong mitigation scenarios (e.g. RCP2.6 and 4.5). In RCP8.5, the two-pattern scaling has a similar performance compared to the univariate pattern scaling (Fig. 3b, e blue lines). It has a larger RMSE than the two-pattern regression model mostly because it overestimates the magnitude of the stabilisation pattern (Fig. 3h green line). The above illustration of the two-pattern scaling model is not intended to be an optimal solution. Instead, the focus is to present a "minimal" model that can explain the distinct behaviours of the transient and stabilisation patterns, and a framework for the two-pattern scaling that can be further developed. For example, the two-pattern scaling could be further improved by optimising the e-folding timescales, or using a different function than an exponential decay to represent the impulse response of DSL to RF. A similar two-pattern scaling model is developed by Yuan and Kopp (2021), in which they use upper and deep ocean temperatures derived from a two-layer energy balance model as 1 ( ) and 2 ( ) , instead of estimating them from RF. In Sects. 6-8, we will present a method to estimate G in Eq. (6) from AOGCMs directly, which forms a more comprehensive linear system model and avoids arbitrary choices of formulation.

SSL and ocean density change
The transient and stabilisation patterns of SdynSL change in RCP4.5 are primarily determined by SSL change (Eq. (1)) integrated over the 0-1000 m depth range (compare Fig. 4a, b with c, d). The 0-1000 m SSL change is equivalent to the baroclinic component of DSL change in Lowe and Gregory (2006) assuming the layer of no motion at 1000 m depth. Differences between SdynSL change and SSL change are evident in the Arctic, the Southern Ocean and marginal seas, suggesting that density change below 1000 m and redistribution of ocean mass are non-negligible in those regions. To shed light on how density changes at different layers affect the patterns of SdynSL change, SSL change is further decomposed into to the 0-200, 200-500 and 500-1000 m components. In the transient pattern, the 0-500 and 0-200 m components dominate the meridional dipoles in the North Pacific and the North Atlantic, respectively, while all the three components are important in shaping the meridional contrast in the Southern Ocean (Fig. 4e, g). In the stabilisation pattern, a uniform change of SSL across latitudes is evident in all the three components (to a lesser extent in the Atlantic), except in the Southern Ocean due mostly to the 500-1000 m layer (Fig. 4f, h). The uniform change of SSL sees stronger contributions from deeper layers in both basins, which is the opposite of the transient SSL change in the Indo-Pacific (Fig. 4e-h).
Ocean potential density change exhibits different depth distributions during the two periods in RCP4.5. It is strongly surface-intensified during the transient period (Fig. 5a, d, and c, f black lines), with the strongest change found in the tropics and northern subtropics (0.4-0.6 kg m −3 century −1 ). Density change tends to be more surface-intensified in the tropics than higher latitudes where subduction and convection quickly carry the surface signal deeper into the ocean. Compared to the transient pattern, the stabilisation pattern is more uniform across depths, with a subsurface maximum within the 200-1000 m depths (Fig. 5b, e). The subsurface maximum of density change has a magnitude of 0.04-0.08 kg m −3 century −1 , and is stronger in the Indo-Pacific than the Atlantic. It corresponds to the depths of the minimum salinity layer in the Southern Ocean (indicated by the 34.5 PSU contour in Fig. 5b, e), a signature of the Antarctic Intermediate Water (Sallée et al. 2013). In addition, the stabilisation pattern has a much weaker surface signature than the transient pattern (by 0.2-0.3 kg m −3 century −1 ). The main features of both the transient and stabilisation patterns are consistent across models, but the exact shapes and magnitudes of those features are model-dependent, likely due to the spread of forced response across models (Fig. 5c, f).

Method of emulating SdynSL change in linear system models
In this section, we introduce the general procedures of calculating SdynSL change in a Linear System Model (LSM). This is to lay a foundation for Sects. 7 and 8, in which we examine how accurately the evolving patterns of the CMIP5 SdynSL change can be emulated using LSMs based on the Green's function of AOGCMs. In contrast to Sect. 4 where we prescribe the Green's function of a LSM (i.e. Equation (6)), here we construct Green's functions from AOGCMs directly. In a LSM, change of an ocean state X under the influence of some Boundary Conditions (BCs) can be written as a convolution of a constant kernel G with the BCs: where r is the position vector of ocean interior, r s is the position vector of ocean surface, t elapsed time since an equilibrium state and lead time relative to t . G represents the system's response to an impulse of BC (i.e. the system's Green's functions). The integral is conducted over all surface locations in r s and all time before t (i.e. = 0 → t ).
For example, one can calculate X at 2300 by integrating BCs from 1870 to 2300 using Eq. (10), assuming that the ocean is at equilibrium at 1870 and earlier (i.e. BC = 0 prior to 1870). The BC could be either the Dirichlet BC or the Neumann BC. In this study, the X(r, t) we are interested in is the forced response of SdynSL and ocean density change to a change of RF in the CMIP5 models, and we investigate this using two configurations of LSMs. The first configuration investigates the role of the pre-industrial ocean circulation in shaping the patterns of SdynSL change, and is referred to as the piControl LSM. The information of RF change is added to the system by prescribing Sea Surface Temperature (SST) and Sea Surface Salinity (SSS) anomalies (relative to an equilibrium state) from the CMIP5 models as the Dirichlet BC. The corresponding G is constructed to be a boundary propagator (Haine and Hall 2002) which depicts how an impulse of the SST or SSS anomaly spreads into the ocean interior under the influence of pre-industrial ocean circulation. The temperature and salinity anomalies from the piControl LSM are used to calculate SSL change, which is used to approximate SdynSL change. In the second configuration, we directly use the RF as the BC (i.e. the Neumann BC), and use the ocean's response (e.g. DSL change) to an impulse of CO 2 quadrupling as G. This model is referred to as the abrupt4xCO 2 LSM. One important difference between the piControl LSM and the abrupt4xCO 2 LSM is that only the latter accounts for the impact of change in ocean circulation. Details on how to derive the piControl LSM and the abrupt4xCO 2 LSM, and how their results are compared to the CMIP5 models are presented in the following two sections.

Deriving the piControl LSM
Deriving G for the piControl LSM requires conducting passive tracer release experiment for different surface regions ( r s ) under the pre-industrial condition. Because this is not a standard experiment of CMIP5, we conduct the tracer experiment using the AOGCM HadCM3 (Gordon et al. 2000) and use its G as a substitute for the G of CMIP5 models. The HadCM3 ocean is based on the Cox (1984) model with Horizontal eddy mixing in HadCM3 is parameterised using the Gent and McWilliams (1990) and Redi (1982) schemes. The near surface vertical mixing is parameterised via the Kraus and Turner (1967) mixed layer model. Gregory et al. (2016) showed that under greenhouse gas forcing, HadCM3 generates DSL change and ocean heat uptake that are similar to the CMIP5 models.
To derive G, the surface domain of HadCM3 is separated into 27 patches (Fig. S2) based on its climatological surface density, which turns r s in Eq. (10) into the surface patch index (e.g. Khatiwala et al. 2012). Passive tracer release experiments are carried out for each patch ( r s = 1...27 ) separately in a 300 year integration of HadCM3 under the preindustrial condition (i.e. a control experiment with regard to climate change). In each experiment, the tracer boundary condition is set to unity for the corresponding patch and zero elsewhere throughout the experiment, which gives X(r, t|r s ) = ∫ d 2 r s ∫ t 0 G(r, |r s ) d from Eq. (10). We then derive G(r, |r s ) (the kernel for each patch) by calculating the time derivative of X(r, t|r s ) (the tracer field from each patch experiment). A similar method was used to study anthropogenic CO 2 in Khatiwala et al. (2009) and ocean heat content in Zanna et al. (2019).
The adequacy of using the G of HadCM3 to approximate the G of CMIP5 models is examined by comparing estimated and diagnosed passive heat uptake in their FAF-passiveheat experiments (Gregory et al. 2016), in which a truly passive heat tracer Θ passive (representing heat) is initialised to zero everywhere in the ocean and a common surface flux of Θ passive is applied to each model for 70 years. The FAFpassiveheat experiment does not simulate G in each model directly. Rather, we compare G across models by comparing their Θ passive at the surface (i.e. BC of the piControl LSM) and Θ passive in the interior (i.e. X of the piControl LSM) in FAF-passiveheat. We assume that if both are similar across models, then G is similar across models too.
In the Indo-Pacific, change of Θ passive in HadCM3 after 70 years is in good agreement with four CMIP5 models used here for both surface values (Fig. 6a) and depth-integrals (Fig. 6c), except in 60° S-40° S where the spread across models is more pronounced. In the Atlantic, the agreement Fig. 6 Comparing passive tracer uptake due to the pre-industrial control circulation between HadCM3 and four CMIP5 models in the FAF-passiveheat experiment. Results shown are passive heat tracer Θ passive averaged in year 51-70. a, b Zonal averages of Θ passive at the surface. c, d Zonal and depth integrals of Θ passive . Averages and integrals are done for the Indo-Pacific Ocean (in a, c) and the Atlantic Ocean (in b, d) separately. Results from different models are denoted by different colours between HadCM3 and the CMIP5 models is encouraging south of 30°N (Fig. 6b, d). North of that, a significant spread of Θ passive is found, implying a spread of G. Based on results from the FAF-passiveheat experiment, we conclude that the G of HadCM3 is adequate for approximating the G of CMIP5 models in the Indo-Pacific and the South Atlantic, but less so in the North Atlantic and the Southern Ocean due to the spread of G across four CMIP5 models.
Another caveat of the piControl LSM is that its G is derived for large surface patches, instead of each of surface grid cells, to reduce the computational cost. This choice has the effect of neglecting the covariance between G and BCs within each patch. As a result, the patterns of tracer uptake estimated from the piControl LSM tend to be smoother than that simulated from AOGCMs which fully account for the variability of G and BCs within grid cells. In Ito and Wang (2017), a similar piControl LSM is developed using the G of the global patch derived from an ocean model. Their results showed that the latitude-depth distribution of CFC-11 estimated from the G of the global patch correlates with that from the ocean model at > 0.87 in space. This implies that the dominant patterns of tracers estimated from the piControl LSM are not sensitive to the choice of patches. However, the exact values of tracer fields from the piControl LSM are likely to be different if different setups of patches are used.

Transient and stabilisation patterns in the piControl LSM
The transient and stabilisation patterns of SSL and ocean density change in RCP4.5 are calculated using the piControl LSM and compared with those of the CMIP5 multi-model mean. Specifically, this is done by convolving G derived in Sect. 7.1 with the CMIP5 SST and SSS anomalies in RCP4.5. These anomalies are relative to the 1861-1880 level, and averaged into the 27 patches (Fig. S2). To minimise the unforced variability, the CMIP5 multi-model mean of the SST and SSS anomalies, rather than those from individual models, are used as the BC of the piControl LSM. The use of SST and SSS anomalies as BCs has the caveat that they are not entirely due to changes in surface fluxes. They also contain a signal due to redistribution of climatological states resulting from climate change. However, a separation of the redistributed component is not possible with available CMIP5 outputs. In the FAF-all experiment (Gregory et al. 2016) where the added and redistributed SST and SSS changes under a global warming scenario can be separated, we find that the forced component is an order of magnitude larger than the redistributed component in regions outside of the Atlantic and Arctic (not shown). This justifies neglecting the redistributed SST and SSS anomalies in the Indo-Pacific and the Southern Ocean when calculating the BCs.
The piControl LSM captures the dominant features of CMIP5 ocean density change in the Indo-Pacific during both the transient and stabilisation periods of RCP4.5 (Fig. 7a-c). One difference between density change in the piControl LSM and the CMIP5 models is that the former penetrates to greater depths than the latter. This is evident in most latitudes in the transient pattern and south of 30° S in the stabilisation pattern (Fig. 7a, b). Qualitatively, such a difference is expected because the effect of upper ocean warming, which tends to reduce vertical tracer uptake (by stabilising water columns), is neglected in the piControl LSM. Part of the difference in the Southern Ocean is also related to the caveats regarding the piControl LSM, namely that G is derived from HadCM3 and based on large patches.
Differences in density change between the piControl LSM and the CMIP5 models are more pronounced in the Atlantic. In the North Atlantic, the piControl LSM features density changes greater than 0.2 kg m −3 century −1 from the surface down to 3000 m during the transient period ( Fig. 7d black contours), but that depth range is only 0-300 m in the CMIP5 models ( Fig. 7d green contours), probably due to the weakening of deep convection and the slowdown of the Atlantic Meridional Overturning Circulation (AMOC) in the CMIP5 models, both of which transport surface tracers to depth. In the stabilisation pattern, the deep density change seen in the piControl LSM moves further south, with its core reaching 0-40°N and 500-2500 m (Fig. 7e shading). In comparison, the density change in CMIP5 is much weaker and peaks at a shallower depth (Fig. 7e green contours). The piControl LSM simulates greater density change than the CMIP5 models across all depths in the Atlantic (compare solid and dashed lines in Fig. 7f), suggesting a higher uptake efficiency for passive tracers than active heat, consistent with Romanou et al. (2017). Again, part of the differences discussed here might be related to the caveats of the piControl LSM.
The transient and stabilisation patterns of < SSL > change (0-1 km integrated) simulated by the piControl LSM exhibit some resemblances to those of the CMIP5 models in the Indo-Pacific (" < > " denotes global means are removed). But, there are also marked differences between the two estimates, evident in both structures and magnitudes, especially in the Atlantic Ocean. In the transient pattern, the piControl LSM simulates a meridional dipole and zonal bands of < SSL > change in the Indo-Pacific, similar to those in the CMIP5 models (compare Fig. 8a, c). However, those from the piControl LSM have a wider structure and a greater magnitude than their CMIP5 counterparts. The difference is qualitatively similar to the effect of circulation change discussed in Fig. 7b of Lowe and Gregory (2006): a poleward redistribution of density anomalies.
In the Atlantic, the piControl LSM features a hemispheric contrast in the transient pattern of < SSL > change ( Fig. 8c), which is absent in CMIP5 (Fig. 8a). Similar patterns of < SSL > change are also found in  when the effect of ocean circulation change is excluded in a climate model. The absence of the hemispheric contrast in the CMIP5 models is likely due to the effect of the AMOC slowdown, which tends to oppose the contrast in the piControl LSM [see Fig. 7b in Lowe and Gregory (2006) and Fig. 9d, h in ]. Moving to the stabilisation pattern of < SSL > change, the piControl LSM simulates a more uniform pattern compared to its transient counterpart (Fig. 8 second  row), which agrees with the CMIP5 results ( Fig. 8 top  row). Apart from that, the piControl LSM simulates weaker < SSL > change in the Indo-Pacific and the Southern Ocean (by 1-2 cm century −1 ), and different patterns of < SSL > change in the Atlantic and the Arctic, compared to those from the CMIP5 models. Nonetheless, the finding that the piControl LSM is able to capture some transient and stabilisation features of the CMIP5 < SSL > change suggests that the climatological ocean circulation plays a role in shaping the patterns of SdynSL and ocean density change, especially in the Indo-Pacific.

Timescale decomposition
One benefit of approximating the CMIP5 response using a LSM is that it allows a clear separation of response into different timescales. Such a separation is often difficult to conduct in conventional AOGCM simulations. Here, we use the piControl LSM to identify how responses on various timescales contribute to ocean density and SSL change over the transient and stabilisation periods of RCP4.5. Responses on the timescale t 1 to t 2 are derived by integrating Eq. (10) over from t 1 to t 2 : For example, the response on the 0-20 year timescale is calculated as the response to BCs from current time up to 20 years ago. The timescale decomposition is applied to ocean density and SSL change averaged over the Indo-Pacific Ocean in RCP4.5, for which we find that the piControl LSM is a good approximation of the CMIP5 models ( Fig. 7 top row, discussed in Sect. 7.2).
In the Indo-Pacific, the transient vertical profile of ocean density change in RCP4.5 is dominated by the response on timescale < 20 years, which peaks at the surface, decaying exponentially from − 0.25 kg m −3 to almost 0 at 1000 m ( Fig. 9a). The response on timescale of 20-100 years is of secondary importance and an order of magnitude weaker than the 0-20 year component. It features a subsurface maximum at 500 m depth and is negligible at the surface and below 1500 m (Fig. 9a). The response on longer timescale (> 100 years) is negligible for the transient profile of ocean density change in the Indo-Pacific. Similarly, the basin-averaged SSL change over the transient period of RCP4.5 is also dominated by the response of 0-20 year timescale. The piControl LSM simulates SSL change of 15 cm (13 cm in CMIP5), with over 70% (10 cm) coming SSL and ocean density change over the stabilisation period of RCP4.5 are determined by responses over a wider range of timescales compared to the transient one. In density change, the 0-20 year component remains surface-intensified during the stabilisation period, but is weaker by 50% compared to its transient counterpart (Fig. 9b). In contrast, the responses on longer timescales (20-100, 100-200 and 200-300 years) see an increase during the stabilisation period relative to their transient signals, most pronounced for the 100-200 year component (Fig. 9b). These long-term responses all peak at subsurface layers with deep-reaching signals penetrating to well below 1500 m. As the timescale increases, the density change peaks at a greater depth (with a weaker magnitude) and exhibits a longer tail into the deep oceans. The piControl LSM simulates 25 cm of SSL change averaged over the Indo-Pacific during the stabilisation period of RCP4.5, which is close to the value of 28 cm in CMIP5 (Fig. 9d). The responses on timescales of 0-20, 20-100, 100-200 and 200-300 years account for 22%, 27%, 34% and 17%, respectively of the total in the piControl LSM (Fig. 9d).
From the perspective of response timescale, the transition from the transient pattern to the stabilisation pattern is caused by a reduction of the fast response (timescale < 20 years) and an increase of the slow (and deeper) response (timescale > 20 years). The reduction of the fast response is a consequence of RF stabilisation, which slows the rate of surface change and the fast response associated with it. The fast response is not zero in the stabilisation pattern because surface change does not stabilise with RF, but continues to evolve in high latitudes (Long et al. 2014). On the other hand, the increase of the slow response is in part due to stronger surface change prior to 2100 than 2000, therefore a stronger lagged response during 2100-2300 than 2000-2070. In addition, a longer period used to derive the stabilisation pattern might also contribute to a stronger slow response, as it allows more time for the slow response to develop. The vertical profiles in Fig. 9b indicate that the longer-term responses mostly occur deeper in the ocean where water masses are replenished on decadal and longer timescales (England 1995), consistent with the ocean dynamics. The important role of the 100-300 year component in the SSL change during the stabilisation period explains why stabilising RF or surface change has a limited impact on mitigating SSL change on centennial timescales.

Linear system model based on the abrupt CO 2 quadrupling experiment
The comparison in Sect. 7.2 suggests an important role of ocean circulation change in determining the pattern of SdynSL change, particularly in the North Atlantic. In this section, we compare SdynSL change simulated from the abrupt4xCO 2 LSM, which includes the effect of changes in ocean circulation, with that from the piControl LSM and the CMIP5 models.

Deriving the abrupt4xCO 2 LSM
To derive the abrupt4xCO 2 LSM G, which quantifies evolution of an ocean state after an impulse of RF, we use simulations from the CMIP5 abrupt CO 2 quadrupling (abrupt4xCO 2 ) experiment. Because RF in the abrupt4xCO 2 experiment is a constant, one can evaluate G for a given ocean state by calculating the time derivative of its evolution in the abrupt4xCO 2 experiment relative to the piControl experiment, and normalising it using the RF of CO 2 quadrupling [7 W m −2 based on Andrews et al. (2012)]. For example, the G of DSL change can be estimated as: where DSL 4× and DSL piControl are DSL in the abrupt4xCO 2 and piControl experiments, respectively, and RF 4× is the RF of CO 2 quadrupling. The G is derived for 9 CMIP5 models (Table 1) for DSL and ocean temperature and salinity, and their multi-model mean is used as the G of the abrupt4xCO 2 LSM to reduce the impact of unforced variability. Evolution of ocean states under different RCPs can then be calculated by convolving the abrupt4xCO 2 LSM G with corresponding RFs (Fig. 1a) using Eq. (10). Unlike with the piControl LSM, we cannot distinguish the responses due to RF imposed in different regions using the abrupt4xCO 2 experiment. Therefore, we use the global mean RF as the BC for the abrupt4xCO 2 LSM in Eq. (10) (which removes the integral across r s in it). This approximation holds if the spatial pattern of RF is mostly constant, or the ocean responses are not significantly affected by the time-dependence of the pattern of RF. The abrupt4xCO 2 LSM G is derived directly from the CMIP5 models, which avoids the error due to using the HadCM3 G as a substitute.
A limitation of the abrupt4xCO 2 LSM is that the CMIP5 abrupt4xCO 2 experiment only runs for 140 years, therefore responses on timescale > 140 years cannot be resolved. We expect this limitation would affect the simulation of the abrupt4xCO 2 LSM after 2100, but have little impact prior to that, based on the finding of Sect. 7.3.

SSL and DSL change in the abrupt4xCO 2 LSM
The abrupt4xCO 2 LSM has better performance than the piControl LSM in emulating the patterns of < SSL > change (0-1000 m) in the CMIP5 models. For the transient pattern of RCP4.5, the abrupt4xCO 2 LSM reproduces the meridional dipoles and zonal bands of < SSL > change in CMIP5 accurately (compare Fig. 8 left column). Differences of < SSL > change between the CMIP5 models and the abrupt4xCO 2 LSM are evident in the Arctic, which implies a non-negligible role of the non-linear response there. In the stabilisation pattern of RCP4.5, the abrupt4xCO 2 LSM captures the positive < SSL > change in the Indo-Pacific, the negative < SSL > change in the North Atlantic, and the zonal bands in the Southern Ocean, but all with a weaker magnitude compared to those in CMIP5 (compare Fig. 8 right column). This difference could be partially due to the response on timescale > 140 years, which is not included in the abrupt4xCO 2 LSM. We test this hypothesis by examining whether the agreement of < SSL > change between the two estimates improves for the 2070-2140 period, since in that case a smaller fraction of RF history is neglected in the abrupt4xCO 2 LSM (compared to the case of 2100-2300).
There is indeed an improvement in the agreement of < SSL > change for the 2070-2140 period, which supports our hypothesis (not shown). The abrupt4xCO 2 LSM is also accurate in emulating the CMIP5 transient pattern of DSL change in RCP4.5, but features a weaker stabilisation pattern of DSL change compared to that of CMIP5 (compare Fig. 10 top and bottom rows). The abrupt4xCO 2 LSM has a similar RMSE compared to the univariate pattern scaling in emulating the CMIP5 DSL change during 2000-2100, but it outperforms the univariate pattern scaling during 2100-2300 in all three RCPs ( Fig. 3d-f, green lines). The RMSE of the abrupt4xCO 2 LSM has a slower growth rate than that of the univariate pattern scaling. For example, it increases from 1 to 3 cm in the former, but 1 to 5 cm in the latter during 2100-2300 in RCP4.5. Compared to the two-pattern regression, the abrupt4xCO 2 LSM has a similar RMSE in RCP4.5 and 2.6, and even a lower RMSE in RCP8.5 ( Fig. 3d-f, green lines). The STD of DSL change from the abrupt4xCO 2 LSM starts to diverge from that of the CMIP5 models after 2150 (Fig. 3a-c), associated with an increase of RMSE ( Fig. 3d-f). This is consistent with the caveat that the abrupt4xCO 2 LSM only accounts for 140 years of RF 1 3 history, which leads to an underestimate of long-term DSL change. Despite the caveat, our results confirm that the LSM is a better formula for sea-level projections than the univariate pattern scaling. We expect that the RMSE of DSL change from the abrupt4xCO 2 LSM would be further reduced if a longer simulation of the abrupt4xCO 2 experiment were available.

Discussion
The transient and stabilisation patterns of the ocean's response to RF were also studied in Long et al. (2014Long et al. ( , 2020 for ocean temperature. It is useful to clarify a difference between Long et al. (2014Long et al. ( , 2020) and our study. The two patterns are referred to as "fast" and "slow" responses in Long et al. (2014Long et al. ( , 2020, based on the time lag between a change in RF and a change in ocean temperature. This definition of time lag is however different from in Eqs. (10) and (11) for the piControl LSM, where measures the time difference between a change in SST and a change in ocean interior temperature. Therefore, all changes that occur after RF stabilisation are "slow" responses in Long et al. (2014Long et al. ( , 2020, while that is not the case in Fig. 9d. The time-dependence of the pattern of SdynSL change in RCP4.5 has implications for the univariate pattern scaling. First, it suggests that the constant pattern assumption (e.g. Mitchell 2003;Perrette et al. 2013;Bilbao et al. 2015) is not valid once RF starts to stabilise. Second (and more important), the ratio between DSL change and GMTSL change is not a constant, but is reduced after RF stabilisation (as sealevel rise becomes increasingly uniform). Therefore, applying the transient ratio to stabilisation scenarios overestimates regional sea-level rise (see Fig. 3a-c red lines). To improve on the univariate pattern scaling, the behaviour of RF should be taken into account, for example by using separate patterns and scaling parameters under different RF profiles.
The evolving patterns of ocean density change have implications for the ocean's expansion efficiency of heat ( ), which measures GMTSL rise for a given amount of global ocean heat uptake (Russell et al. 2000;Kuhlbrodt and Gregory 2012). The transition from a surface-intensified response to a subsurface response shown in Fig. 5 indicates that a higher fraction of ocean warming occurs in regions of lower thermal expansion coefficient ( ) after RF stabilisation (this holds for 0-2000 m, below that increases with pressure/ depth). As a result, we expect a reduction of in mitigation scenarios. To verify this hypothesis, we calculate for the twenty-first century and the twenty-third century in the three RCPs, by regressing GMTSL change against global ocean heat content change (Table 2). In RCP2.6, is reduced from 0.127 (the median value) (0.125 -0.128) m YJ −1 during the twenty-first century to 0.102 (0.097-0.113) m YJ −1 during the twenty-third century (uncertainties denote the 25th-75th percentiles), a reduction of about 20%. In RCP4.5, where RF stabilises at a later time than in RCP2.6, the reduction of is also smaller, at about 10% (from 0.130 to 0.117 m YJ −1 ). The reduction of was also reported in Russell et al. (2000) and Körper et al. (2013) under other mitigation scenarios. Interestingly, in RCP8.5 exhibits a slight increase from 0.136 (0.133 -0.137) m YJ −1 over the twenty-first century to 0.137 (0.136-0.144) m YJ −1 over the twenty-third century. This suggests that the effect of ocean warming, which increases and hence (e.g. Hallberg et al. 2013), overcomes the effect of subsurface warming, which reduces , in RCP8.5. We expect that the ocean warming also increases in RCP2.6 and RCP4.5, but it is not strong enough to overcome the effect of subsurface warming due to much weaker RF in them. In other words, the relative importance of the transient and stabilisation patterns determines the sign of the change of . Beyond 2300, we would expect to decrease in RCP8.5 too as the stabilisation pattern becomes dominant.
DSL/SSL change is decomposed into responses to changes in momentum flux and buoyancy flux in many previous studies (e.g. Fukumori et al. 1998;Lowe and Gregory 2006;Köhl and Stammer 2008;Forget and Ponte 2015;Roberts et al. 2016;Gregory et al. 2016). It is helpful to explain differences between our LSM approach and flux perturbation experiments. The piControl LSM does not distinguish among SSL changes due to changes in different fluxes. Instead, the effect of changes of all surface fluxes is combined in the piControl LSM through their imprints on the SST and SSS anomalies (BCs of the piControl LSM). These surface anomalies are carried into the ocean interior differently in the piControl LSM and in the flux perturbation experiments. In the former, the ocean transport is not affected by changes in the surface fluxes, while that effect is fully accounted for in the latter. Similar to the piControl LSM, the abrupt4xCO 2 LSM does not distinguish among DSL/SSL changes due to different forcings. However, it does fully account for the ocean's responses to changes in surface fluxes as in the flux perturbation experiments.

Conclusion
Sterodynamic Sea Level (SdynSL) rise continues for centuries after stabilisation of Radiative Forcing (RF) in the CMIP5 models. However, its regional pattern exhibits different characteristics before and after RF stabilisation. The transient features of SdynSL change, a subject of many studies in the literature, are replaced by a more uniform basinscale pattern after RF stabilisation, although Global Mean Thermosteric Sea Level (GMTSL) rise shows little sign of slowdown. A linear combination of the transient and stabilisation patterns explains a substantial part of Dynamic Sea Level (DSL) change in RCP2.6, 4.5 and 8.5. In all cases, DSL change is dominated by the transient pattern while RF increases rapidly, but it is increasingly affected by the stabilisation pattern after RF starts to stabilise. As a result, the univariate pattern scaling overestimates the magnitude of DSL change after 2100 in RCP2.6 and 4.5. The growth of the stabilisation pattern could last for centuries after RF ceases increasing. Similarly, ocean density change, which accounts for most of DSL change, exhibits different patterns before and after RF stabilisation, with its maximum shifting from surface to subsurface. The two-pattern regression model can be converted to a two-pattern scaling model with its scaling parameters estimated from RF and a two-timescale linear system model. This approach has a better skill than the univariate pattern scaling in emulating the CMIP5 DSL change in RCP2.6 and 4.5 after 2100.
To explain the evolving patterns of SdynSL change, Linear System Models (LSM) are constructed based on the AOGCMs' responses to an impulse of boundary conditions (i.e. their Green's functions). In the piControl LSM, in which the CMIP5 transient ocean circulation is replaced by the HadCM3 pre-industrial ocean circulation, some transient and stabilisation features of the CMIP5 Steric Sea Level (SSL) and ocean density change are still present, particularly in the Indo-Pacific and Southern Ocean. This suggests that the climatological ocean circulation plays a role in determining the pattern of SdynSL change in these areas. The piControl LSM further reveals that in the Indo-Pacific, the basinaveraged ocean density and SSL change is dominated by the decadal (0-20 years) response during the transient period of RCP4.5, while they are shaped by the decadal to centennial responses jointly during the stabilisation period. Nonetheless, differences in SSL change between the piControl LSM and the CMIP5 models are evident in many regions (particularly the North Atlantic), which implies that the effect of ocean circulation change is non-negligible.
To address this issue, the abrupt4xCO 2 LSM is developed, which uses the CMIP5 models' responses to an impulse of RF as its kernel. The abrupt4xCO 2 LSM is more Table 2 The ocean's expansion efficiency of heat during the twenty-first century and the twenty-third century and the difference between the two Results shown are the median value and the 25th-75th percentiles in unit of m YJ −1 (1 YJ = 10 24 J) accurate than the piControl LSM in emulating the CMIP5 SSL change, which confirms the important role of ocean circulation change in shaping the patterns of SdynSL change. The abrupt4xCO 2 LSM also reproduces the transient pattern of DSL change in CMIP5 to a large extent (RMSE ≈ 1 cm), but its stabilisation pattern of DSL change has a weaker magnitude than CMIP5, because it only accounts for 140 years of RF history (due to data availability). Importantly, the abrupt4xCO 2 LSM has better performance than the univariate pattern scaling in emulating the CMIP5 DSL change beyond 2100 in three RCPs examined here.
Our results have implications for the ocean's expansion efficiency of heat ( ). The growth of the transient response tends to increase , while that of the stabilisation response tends to reduce . Their relative importance determines the sign of the change of . In three RCPs examined here, the change of over 2000-2300 is on the order of 1-10% relative its twenty-first century value. Given the performance of the abrupt4xCO 2 LSM in emulating the CMIP5 DSL change, it would be useful if modelling groups could extend the abrupt4xCO 2 experiment to 500 years, so that most of RF history can be accounted for in the abrupt4xCO 2 LSM for simulations up to 2300.
Some configurations of the CMIP5 models are based on the understanding of the current climate system. At high levels of warming in the future (e.g. in RCP4.5 and 8.5), it is unclear whether the approximations or parameterisations of the models will hold if the mean states of the ocean and atmosphere undergo a major change. This might affect the realism of CMIP5 simulations after 2100. How realistic the CMIP5 projections are on multi-century timescales is a challenging problem that is well beyond the scope of this study. Instead, our effort focuses on gaining a better understanding of the CMIP5 models, and using that knowledge to provide a computationally efficient method to emulate their behaviours, such as linear system models.

Availability of data and material
CMIP5 model data are available from the Earth System Grid Federation data portal at https:// esgf-node. llnl. gov/ proje cts/ cmip5/. Green's functions simulated in HadCM3 are available upon request.