Atmosphere-ocean general circulation models (GCMs) are the primary tool used to project ocean dynamic sea level (DSL) change, but the computational demands of these models limit the utility of GCM ensembles for estimating the likelihood of different levels of future sea-level change. Ensembles such as the Coupled Model Intercomparison Project Phase 5 (CMIP5, Landerer et al., 2014;Taylor et al., 2012) are composed of models contributed based on voluntary effort, not the product of systematic experimental design; as such, they are an "ensemble of opportunity" rather than a probabilistic ensemble (Tebaldi & Knutti, 2007). The CMIP future projection experiments are driven by a small number of forcing scenarios-Representative Concentration Pathways (RCPs) in the case of CMIP5-and model simulations are of different lengths; some simulations run the RCPs to the year 2100, while others extend these to 2300.

The computationally intensive nature of GCMs makes it challenging to produce large perturbed-physics ensembles that represent uncertainties in key feedback parameters, as well as to simulate forcing conditions intermediate between the RCPs. Simple climate models (SCMs) provide an alternative tool for estimating the uncertainties of future projections at the global scale, as they can capture the overall physics of climate evolution and can be run very fast even on a personal computer (Held et al., 2010;Meinshausen et al., 2011;Millar et al., 2017;Perrette et al., 2013). However, SCMs represent the climate at a highly aggregated (e.g., global or hemispheric) scale, and thus cannot produce spatial patterns of climate change at each time step.

Pattern scaling approaches are often used to translate the global mean surface air temperature (GSAT) change into regional-scale changes for impact analysis (Mitchell, 2003;Rasmussen et al., 2016;Santer, 1990;Tebaldi & Arblaster, 2014;Tebaldi et al., 2011). Generally speaking, pattern scaling uses a simple statistical model (often, linear regression) to relate local climatic changes to a variable such as GSAT change, assuming the patterns of local response to external forcing remain constant under increased forcing (Tebaldi & Arblaster, 2014). Some previous studies use the pattern scaling approach to estimate the uncertainty in DSL projections (Bilbao et al., 2015;M. D. Palmer et al., 2020;Perrette et al., 2013). For example, Perrette et al. (2013) regressed DSL change on GSAT. At New York City, they found that r 2 values across models vary between 0.02 and 0.85, and also that the linear relationship between DSL and GSAT becomes weaker after the 21st century. Bilbao et al. (2015) examined the relationship between DSL and several variables, including GSAT, global-mean sea-surface temperature, ocean volume mean temperature, and global-mean thermosteric sea-level rise (GMTSLR). They found that GSAT performed best in predicting 21st-century DSL change in a high emissions scenario (RCP 8.5), while ocean-volume mean temperature and GMTSLR performed better in lower emissions scenarios (RCP 2.6 and 4.5). They speculated that this difference reflects a more important role for surface warming relative to deep warming in a more strongly forced scenario. They found that, across models and scenarios, area-weighted average root mean square error in pattern-scaled 2081-2100 DSL change ranged from ∼1 to 3 cm. Building upon Bilbao et al. (2015)'s speculation about the relative importance of shallow and deep warming under different scenarios, we developed a bivariate pattern scaling, which uses a multiple linear regression with two predictors: GSAT and global-mean deep ocean temperature change. The two temperature changes can be generated by a two-layer energy-balance model (2LM) (Held et al., 2010;Winton et al., 2010), which has proved to be a useful tool for understanding the responses of climate system to climate forcing (Geoffroy et al, 2013a(Geoffroy et al, , 2013b)). Shallow and deep temperatures from a 2LM have previously been employed in an emulator to extend 21st century CMIP5 projections of GMTSLR to 2300 (M. D. Palmer et al., 2018), and M. D. Palmer et al. (2020) used GSAT from the two-layer model and univariate pattern scaling (based on GSAT) to emulate CMIP5 projections of DSL change.

In this study, we develop an emulator for DSL changes using both GSAT and deep-ocean temperature change projected by a 2LM. Here we drive the 2LM with radiative forcings from the Finite Amplitude Impulse Response model (FaIR), a SCM which includes a reduced-complexity carbon cycle and calculates atmospheric CO 2 concentrations, radiative forcing and temperature changes based on emissions (Millar et al., 2017;Smith et al., 2017Smith et al., , 2018)). FaIR was designed to more accurately reflect the temporal evolution of GSAT in response to a pulse emission, and it has been used in previous studies to produce observation-constrained future projections (Millar et al., 2017;Smith et al., 2017Smith et al., , 2018)). In this study, we develop an emulator for GMTSLR and DSL change using surface and deep-ocean temperature changes generated by the FaIR-2LM (Section 2.2). As the univariate pattern scaling fails to capture the delayed response of the deep ocean to warming, we employ FaIR-2LM and two-layer pattern scaling to project future DSL changes, taking into account uncertainty in climate sensitivity, and demonstrate their ability to interpolate between climate scenarios run by GCMs. Compared to M. D. Palmer et al. (2018Palmer et al. ( , 2020)), which also use a 2LM to emulate GMTSLR or DSL projections, our approach differs in: (1) employing radiative forcings calculated based on emissions; (2) applying a format of 2LM considering efficacy factor of deep ocean heat uptake; (3) using both surface and deep-ocean temperature for pattern scaling (more details are described in supporting information).

Section 2 describes data and methodology, including the details of FaIR-2LM, the calibration of the FaIR-2LM based on selected CMIP5 GCMs, the two-layer pattern scaling methodology, and the application of this system to emulate DSL projections. Section 3 evaluates the performance of the two-layer pattern scaling. Section 4 shows the resulting ensemble of DSL projections. Finally, Section 5 discusses and summarizes the results.

We use the zos variable from five CMIP5 general circulation models (GCMs) in RCP 2.6, 4.5, and 8.5 scenarios: MPI-ESM-LR, bcc-csm1-1, HadGEM2-ES, GISS-E2-R, IPSL-CM5A-LR. These five GCMs are used because they were used to calibrate the parameters of the 2LM by Geoffroy et al. (2013) and provide multi-century data (to 2300) for zos in all three scenarios. DSL is taken as zos with its global mean removed, consistent with the definition of Gregory et al. (2019). The drift is removed from DSL by subtracting a linear function of time fitted to the pre-industrial control simulation from each scenario experiments, at each grid point. In addition, we remove the climatology in a baseline period (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) from DSL. The global mean surface air temperature (GSAT) and GMTSLR from the five models in the three scenarios are also used to evaluate the performance of FaIR-2LM.

This study develops a hybrid SCM model by replacing the temperature module in FaIR 1.3 (Smith et al., 2018) with a 2LM. In FaIR 1.3, GSAT changes are the sum of two components, representing fast and slow responses to effective radiative forcing (ERF) (Equation 22in Smith et al., 2018). The fast and slow components of temperature changes in FaIR 1.3 mathematically depend on multiple coefficients (e.g., thermal response timescales) that are obtained from the ensemble mean of multiple CMIP5 models (Geoffroy, et al., 2013b). Since these components do not have an unambiguous physical meaning, it is challenging to link them to sea-level change. Therefore, we replace the temperature module in FaIR 1.3 by the 2LM to construct FaIR-2LM. In each step of FaIR-2LM, the 2LM is driven by radiative forcing from FaIR 1.3, and produces the GSAT anomaly, which feeds back to the FaIR carbon cycle (Figure S1).

We employ a 2LM that includes an efficacy term for deep ocean heat uptake (Geoffroy, et al., 2013a;Held et al., 2010;Winton et al., 2010):

(1)

where denotes the adjusted radiative forcing, C and C 0 are the heat capacity of the well-mixed upper layer and the deep ocean layer, respectively, and T and T 0 represent the global mean temperature anomalies of the upper and lower layer, respectively. Following Equation 22 in Geoffroy, et al. (2013b) and using C = 8.2 W yr m -2 K -1 and 0 C 109 W yr m -2 K -1 based on an average across multiple CMIP5 GCMs (Geoffroy, et al., 2013a), we estimate the average depths of the upper layer and lower layer are 86 m and 1141 m, respectively. T is equivalent to GSAT perturbation (Held et al., 2010). is the parameter for cli-mate feedback, is the coefficient of deep ocean heat uptake, and is the efficacy factor of deep ocean heat uptake, which represents the uneven spatial distribution of heat exchanges between the two layers.

To calibrate FaIR-2LM, we adjust parameter settings (listed in Table 1) based on previous studies (Forster et al., 2013;Geoffroy, Saint-Martin, et al., 2013a;Zelinka et al., 2014). The radiative forcing in FaIR-2LM is driven by the default emission trajectory for each scenario in FaIR 1.3, but scaled by two parameters determined for each GCM: (1) the radiative forcing of CO 2 doubling ( 2 2

) reported by Forster et al. (2013), and (2) the present-day aerosol forcing (af pd ) estimated in previous studies (Forster et al., 2013;Zelinka et al., 2014), or -0.9 W m -2 -the median of range estimated by the Fifth Assessment Report of Intergovernmental Panel on Climate Change (IPCC AR5) (Stocker et al., 2013)-for models not reported in previous studies. The five parameters in Equations 1 and 2 (i.e., , , , C, C 0 ) are the same as those in Geoffroy et al. (2013) for the corresponding GCMs.

GSAT produced by the calibrated FaIR-2LM is compared with that from the corresponding GCMs in the three scenarios (Fig. S2). For the five GCMs, the GSAT simulated by FaIR-2LM is close to the GSAT from the corresponding GCM, with the root mean square error (RMSE) determined over the entire simulation period in a range of 0.15-0.23 K for RCP2.6, 0.14-0.32 K for RCP4.5, and 0.20-0.43 K for RCP8.5.

GMTSLR is driven by the thermal expansion of sea water volume due to the increase in ocean heat uptake. To calibrate GMTSLR in FaIR-2LM to match a specific GCM, we first correct the drift in the GCM's GMT-SLR field by removing the linear trend in the pre-industrial control simulation, assuming the drift is not sensitive to the external forcing (Hobbs et al., 2016).

Then, we emulate GMTSLR based on the T and T 0 from FaIR-2LM following the approach described in Kuhlbrodt and Gregory (2012):

where is the expansion efficiency of heat in units of 10 -24 m J -1 . The value is calibrated by optimizing GMTSLR emulated from FaIR-2LM to match the GMTSLR simulated from the corresponding GCM.

Univariate pattern scaling is based on a linear relation between regional changes in a climate variable (DSL for this study) and global mean responses of climate change, such as GSAT (T ):

, , , , , , DSL t x y x y T t b x y t x y (4)

where x and y denote longitudes and latitudes, t represents the time, b is an intercept term, and is the residual term. Here, α captures the scaling relationship between DSL and GSAT (Figure 1). The five GCMs agree that the linear response of DSL to surface warming is positive over the Arctic and sub-polar Atlantic, and negative over the southeastern Pacific and the southern areas of Southern Ocean.

In the bivariate pattern scaling approach, we regress the DSL anomaly on both T (GSAT anomaly) and 0 T (deep-ocean temperature anomaly) from FaIR-2LM:

where i t denotes years in three scenarios, i = 1, 2, 3. For each GCM, we estimate the fields of α, β, b, and ε by regressing projections from all three emissions scenarios (RCPs 2.6, 4.5, and 8.5) on T and 0 T on a grid cell- GCMs.

Wm K , ,

Wyrm K C , and Geoffroy et al. (2013). The units for F 2 × co 2 and afpd are W m -2 same period (Figure 1). Consistent with the univariate scaling pattern, the five GCMs agree that the upper-layer response, represented by , is positively correlated with warming over the most areas of Arctic and northern edge of the Southern Ocean, and negatively correlated with warming over the southeastern Pacific YUAN AND KOPP 10.1029/2020MS002323 5 of 17 and the southern areas of Southern Ocean. The deep-layer response represented by is positively correlated with warming over the Indian and tropical and southern Pacific Oceans, and negatively correlated with warming over most areas of the Southern Ocean and Arctic. These reflect opposite behaviors between rapid and sustained changes in DSL over the Arctic, the Indian and tropical and southern Pacific Oceans, and a consistent DSL fall in both rapid and sustained changes over the Southern Ocean.

There is little agreement on either surface-or deep-layer slopes across the five GCMs over most parts of the Atlantic basin (Figure 1). This may reflect limited skill in simulating strong western boundary currents (e.g., the Atlantic Meridional Overturning Circulation (AMOC)) in the GCMs, which have a relatively coarse (∼1°) spatial resolution in ocean component (Small et al., 2014) and so poorly capture non-linear mesoscale processes in the ocean current (Hallberg, 2013). Near the eastern coast of North America, DSL is closely related to AMOC (Goddard et al., 2015), which is expected to weaken in a warming climate (Caesar et al., 2018). Low skill in capturing AMOC behavior can affect the DSL projections in the Atlantic basin as well as its coasts (van Westen et al., 2020). As the coefficients of pattern scaling depend on the simulations by the GCMs, they also do not explicitly resolve the non-linear mesoscale process of the ocean current. Therefore, we should interpret the DSL changes predicted by the two-layer emulator with cautions over the regions where non-linear mesoscale effects of ocean current are strong.

We use two steps to generate a probabilistic ensemble of DSL projections. First, we generate an ensemble of surface and deep-ocean temperature pairs using FaIR-2LM. The planetary energy balance at the top of the atmosphere (Zelinka et al., 2020) is:

where N is the radiative imbalance at the top of the atmosphere. The equilibrium climate sensitivity (ECS) is given by T when N = 0, and =

2 2 CO F . Therefore, is related to 2 2 CO F and ECS by 2 2 / CO F ECS (7) The uncertainty of 2 2 CO F is small relative to the spread of , while ECS largely determine the uncertainty of . Therefore, we adopt the best estimation in the Intergovernmental Panel on Climate Change Fifth Assessment Report (AR5) for 2 2 CO F = 3.71 2 W m (Collins et al., 2013). We produce initial distributions of ECS, , and based on the literature constraints (Fig. S4) outlined below: ECS: Based on multiple lines of evidence, the uncertainties of ECS estimated by AR5 are likely in the range 1.5°C-4.5°C with high confidence, extremely unlikely less than 1°C and very unlikely greater than 6°C (Collins et al., 2013). In the AR5 terminology, likely denotes a probability of at least 66%, very unlikely a probability of less than 10%, and extremely unlikely a probability of less than 5% (Mastrandrea et al., 2010). Therefore, we construct a log-normal distribution for ECS with parameters optimized to match a 5th percentile of 1°C, a 17th percentile of 1.5°C, an 83rd percentile of 4.5°C, and a 90th percentile of 6°C.

: We treat as normally distribution, with mean 0.67

W m K and standard deviation 0.15

W m K derived from the 16 GCMs in the CMIP5 archive (Geoffroy et al., 2013).

: As the efficacy factor of heat uptake is related to deep-ocean heat uptake (Held et al., 2010), we use instead of to maintain the covariance between and . We calculate the mean of 0.86

W m K and standard deviation 0.29

W m K of based on the products of and from 16 GCMs in CMIP5 archive (Geoffroy, et al., 2013a), The distribution of is constructed as a normal distribution with the multi-model mean and the multi-model standard deviation.

TCR: Under a zero-layer approximation which considers the 1%/yr increase in CO 2 until doubling scenario occurring on a timescale long enough that the upper ocean is in approximate equilibrium and short enough that the deep-ocean temperature has not yet responded substantially, Transient Climate Response (TCR) can be obtained by (Jiménez-de-la-Cuesta & Mauritsen, 2019): W m K based on the multi-model mean of GCMs from Coupled Model Intercomparison Project Phase 5 (CMIP5) archive (Geoffroy, et al., 2013a).

We then generate a 100,000-member ensemble of ECS, γ and γε based on these distributions via Monte Carlo sampling. As should be larger than 0, we discard parameter sets in which < 0 or > 2 0.86 to keep the mean of in parameter sets to be 0.86

W m K . Therefore, 99734 parameter sets are kept. An ensemble of is then computed by the best estimation of 2 2 XCO

and the ensemble of ECS based on Equation 6 (Fig. S4). The median (central 66% range) of is -1.39 (-2.4 to -0.8)

W m K . As the likely range of ECS estimated by AR5 is equivalent to the central 90% range of ECS estimated by CMIP5 GCMs, the uncertainty range of estimated by FaIR-2LM is larger than that estimated by ensemble of GCMs (Geoffroy et al., 2013). The spread of TCR is estimated as a diagnostic by substituting the ensemble of , , and best estimation of 2 2 XCO

into Equation 8. The uncertainty of TCR is in a central 66% range of 1.1-2.3 C, with a 95th percentile of 2.9 C. This is consistent with but slightly narrower than the TCR estimated by AR5, which is likely between 1 C and 2.5 C, and is extremely unlikely greater than 3℃.

We apply Latin hypercube sampling (LHS, Stein, 1987) to the parameter sets of , , by sampling 1,000 sets from the 99734 parameter sets. For each parameter, LHS divides the probability density function of the 99,734 samples into 1,000 portions that have equal area. A sample is taken from each portion randomly so that the 1,000 sample sets cover the multidimensional distribution of the three parameters. Finally, we applied 1,000 parameter sets together with the fixed parameters ( 2 2

, C, C 0 ) to the FAIR-2LM and generate a 1,000-member probabilistic ensemble of temperature pair time-series.

We compare the spread in GSAT projected by FaIR-2LM with the likely ranges estimated by AR5 for four different periods (Collins et al., 2013) (Table 3 and Figure S5). The mean of the probabilistic ensemble is slightly lower than the mean estimate of GSAT from AR5 in all four periods of RCP2.6 and RCP4.5, and in the 21st century for RCP8.5. Compared with AR5 likely ranges, the central 66% probability range of GSAT from FaIR-2LM is generally consistent: narrower in all four periods of RCP2.6, narrower in the first two periods but wider in the last two periods in RCP4.5, and wider in the first two periods but narrower in the last two periods in RCP8.5.

We project GMTSLR based on Equation 3 using the probabilistic ensemble of surface and deep-ocean temperature projections from FaIR-2LM. The C, C 0 and expansion efficiency of heat ( 24m 0.113 10

) used here are adopted from the multi-model ensemble mean of CMIP5 archive (Geoffroy, et al., 2013a;Kuhlbrodt & Gregory, 2012).

A projection of DSL is constructed as follows: 1) a pair of and is randomly picked with replacement from the pool of two-layer patterns produced in Section 2.3; 2) a temperature pair from the 1000 members is combined with the pair of and in an equation:

0 , , , , , DSL t x y x y T t x y T t b x y (9)

To evaluate the prediction skill of the two-layer pattern scaling, we compare the DSL changes simulated from a GCM with the DSL changes emulated by the two-layer pattern scaling ( DSL) based on FAIR-2LM using the key parameters (i.e. parameters in Table 1) in the same GCM. Two metrics are used: (1) absolute values of the residual differences between DSL changes and DSL changes during a period at each grid point, and ( 2) global average of the absolute values obtained from the metric 1 (Table S1). These two metrics are applied to both bivariate pattern scaling and univariate pattern scaling, to examine the improvement of bivariate approach comparing with the univariate approach.

In 2271-2290, for instance, the global-averaged climatology of DSL DSL (score obtained by the second metric) from the two-layer pattern scaling is less than that from the univariate pattern scaling (bottom row

Figure 2), with a reduction of 36%, 24%, and 34% in RCP2.6, RCP4.5, and RCP8.5, respectively. The spatial pattern of R = DSL DSL is derived from both approaches are various across GCMs (Figures S6-S10). The 5-model ensemble averaged climatology of R in both approaches is higher over high latitudes (e.g., Arctic, subpolar Northern Atlantic, Southern Ocean) than over middle to low latitudes, but is generally lower in two-layer pattern scaling than in univariate pattern scaling (first two rows Figure 2). As the pattern-scaling method cannot resolve DSL change due to unforced variability, the relatively large |R| over the high latitudes may be due to the relatively high unforced variability over these regions (Bilbao et al., 2015).

We further compare the time evolution of DSL predicted by the two-layer pattern scaling approaches with the evolution of DSL in corresponding GCMs through the period 1981-2290. As case studies, we pick two grid cells: one in the western Pacific near the Philippines (14.5°N, 127°E), and the other over the North Atlantic near the coast of New York City [NYC] (40°N, 73°W) (solid black dots in Figure 2). The grid point near the Philippines is selected because it is in the tropical Pacific, where DSL rise associated with the deep ocean temperature rise is strongest, while the grid point near the NYC is selected represent a coastal area that some projections find experiences significant DSL changes in response to changes in AMOC.

At the western Pacific grid cell, in RCP 2.6, the relationship between DSL and GSAT anomaly displays a hook-like shape, indicating continued rise in DSL as GSAT stabilizes and declines in response to negative emissions (Figure 3a). The delayed adjustment of DSL may be due to the continuous warming of deep layer (T 0 ) when GSAT is stabilized, because the ocean is not yet equilibrated with the elevated forcing. In response to changes in T 0 , the deep ocean density is still changing even without changing circulation in the deep ocean (England, 1995), so DSL continues to change. This hook-like shape is captured by the two-layer pattern scaling approach but not by the univariate pattern scaling. Compare to the DSL simulated by a GCM, the RMSE of the predicted DSL is smaller if using the two-layer pattern scaling approach than using the univariate pattern scaling approach. The average RMSE across the five GCMs is reduced 26% if we use the approach from the two-layer pattern scaling instead of the univariate pattern scaling (Table 3). Across the five GCMs, although the relationship between DSL and GSAT is diverse in RCP4.

5 and RCP8.5, DSL YUAN AND KOPP 10.1029/2020MS002323 9 of 17 projected by the two-layer technique is consistently closer than that predicted by the univariate technique to the DSL simulated by the GCMs. The sole exception is for bcc-csm1-1 in RCP8.5, for which the simulated DSL projection is quite linearly associated with GSAT. The average of RMSEs across the 5 GCMs decrease from the univariate pattern scaling approach to the two-layer pattern scaling approach by 35% in RCP4.5 and by 33% in RCP8.5 (Table 3).

At the North Atlantic grid cell, the relationship between DSL and GSAT also displays non-linear features for all the five models, especially in low-and moderate-emission scenarios (Figure 3b). These non-linear features, which may arise from the delayed response of deep branch of AMOC, cannot be captured by univariate pattern scaling but are captured to a large extent by the two-layer pattern scaling (lines in Figure 3). The value of the two-layer approach is highlighted by the clear non-linearity of the DSL response when viewed as a function of GSAT anomaly. Compare to the univariate approach, the RMSE between the DSL simulated by GCMs and the DSL predicted by the two-layer pattern scaling is smaller, with a reduction of 19%, 16%, and 13% for RCP2.6, RCP4.5, and RCP8.5, respectively. The method of two-layer pattern scaling generally has a better performance in emulating the DSL from the corresponding GCM than the univariate pattern scaling, as the two-layer pattern scaling includes one more predictor than univariate pattern scaling, allowing it to capture the delayed adjustment of DSL. The delayed adjustment of DSL due to the deep-ocean warming is important for the DSL projections, as different areas present different features that may reveal the regional variation in deep-ocean circulation (Held et al., 2010).

The procedure described in Section 2.4 allows us to produce a 1000-member probabilistic ensemble of DSL projections not only for the three CMIP5 scenarios: RCP2.6, RCP4.5, and RCP8.5, but also for any other scenarios with an emission pathway between these three scenarios. We demonstrate this capability using SSP3-7.0, a CMIP6 scenario that has forcing intermediate between RCP4.5 and RCP8.5 (O'Neill et al., 2016) and is closer than either to no-policy reference scenarios from most integrated assessment models (Riahi et al., 2017). The emission pathway of SSP3-7.0 used to drive the FaIR-2LM is taken from the Reduced Complexity Model Intercomparison Project (Nicholls et al., 2020).

The five projections using parameters calibrated to the five GCMs respectively are within the 66% range of the 1000-member ensemble for both surface and deep-ocean temperature in the three RCPs (Figure 4). By 2300, the median estimates (66% range) of the surface temperature relative the period of 1986-2005 are 0.5 C (0.2-1. 0 C) for RCP2.6, 2. 2 C (1.2-3.6 C) for RCP4.5, 7. 4 C (4.5-11. 7 C) for RCP8.5, and 5. 3 C (3.2-8.6 C) for SSP3-7.0.

Based on the projections of temperature pairs, we also produced projections of GMTSLR for the four scenarios (Figure 4). The spread of GMTSLR ensemble encapsulates the GMTSLR time series from the 5 GCMs (Figure 4). During the period of 2081-2100, the median estimates (66% range) of the GMTSLR relative the period of 1986-2005 are 0.12 (0.07-0.18) m for RCP2.6, 0.16 (0.10-0.24) m for RCP4.5, 0.19 (0.12-0.27) m for SSP3-7.0, and 0.24 (0.15-0.34) m for RCP8.5. This compares to Oppenheimer et al. ( 2019)'s projected median estimates (66% ranges) of 0.14 (0.10-0.18) m for RCP2.6, 0.19 (0.14-0.23) m for RCP4.5, and 0.27 (0.21-0.33) m for RCP8.5. By 2300, the median estimates (66% range) of GMTSLR relative to the period of 1986-2005 are 0.20 (0.12-0.33) m for RCP2.6, 0.43 (0.25-0.68) m for RCP4.5, 0.85 (0.50-1.33) m for SSP3-7.0, and 1.15 (0.69-1.76) m for RCP8.5. As climate warms, the projections of DSL changes increase along with the increase in GMTSLR. In the five GCMs, although the contributions of DSL changes to the local sterodynamic sea level (DSL + GMTSLR; Gregory et al., 2019) changes are small at some locations (i.e., regions marked by the light and dark gray shadings in Fig. S11), in others the ratio of DSL change with respect to the GMTSLR changes are fairly significant. For instance, DSL changes at some regions (e.g., Arctic, North Atlantic, and Southern Ocean) are greater than 50% of the GMTSLR during the period of 2271-2290 (identified by yellow contours in Fig. S11).

Compared with the GSAT and GMTSLR spread in 2300 estimated by M. D. Palmer et al. (2018), the FaIR-2LM projections have a slightly lower median for all the three RCPs. The 66% range of both surface temperature and GMTSLR estimated by FaIR-2LM is comparable to the 90% range of that estimated by M. Comparing the DSL projections between the period of 2081-2100 and the period of 2271-2290 (Figure 5), the median estimate is lower and the 66% range of uncertainty is narrower at the end of 21st century than that at the end of 23rd century in moderate-to high-emission scenarios (RCP4.5, SSP3-7.0 and YUAN AND KOPP RCP8.5). But in RCP2.6, the median estimate and 66% uncertainty range are comparable in magnitude between these two periods. In both periods, the median DSL anomaly projections across the four scenarios share many similar features (Figure 5). Over the Arctic region, a weak increase in DSL is observed over the Chukchi Sea and the Beaufort Sea in RCP2.6. In the higher emission scenarios, the increase in DSL extends to the whole Arctic basin with intensified amplitudes. The changes in DSL over the North Atlantic are dominated by a negative anomaly under RCP2.6, and display positive anomalies over much of the North Atlantic under RCP8.5 and SSP3-7.0. The ensemble spread of the 5th-95th range of DSL projections are relatively large over the Southern Ocean, Arctic and Subpolar Atlantic than other areas.

The large uncertainties over these areas, consistent with previous literatures (M. D. Palmer et al., 2020;Perrette et al., 2013;Yin, 2012), may be interpreted by the diverse characteristics simulated by GCMs that do not explicitly resolve non-linear mesoscale processes of the ocean current over these areas (van Westen et al., 2020).

At the illustrative grid point near Philippines over western Pacific (Figure 6a), the 66% range of the probabilistic ensemble encapsulates DSL projections from 2 of the 5 GCMs in the three RCPs, while the 90% range of the probabilistic ensemble contains DSL projections from all the 5 GCMs in the three RCPs, except for HadGEM2-ES in RCP2.6. At the grid point near NYC, the projected DSL changes estimated by the probabilistic ensemble exhibits a fat tail, with a median projection in RCP 8.5 of 0.13 m and a 95th percentile projection of 0.8 m by the end of 23rd century. By contrast, RCP 2.6 exhibits a much narrower range, with a median of 0 m and a 95th percentile of 0.08 m. The 66% range of the projected DSL uncertainties encapsulates 2 of 5 GCM projections. The 90% range of the probabilistic ensemble only encapsulates the DSL projections from three over five GCMs in RCP2.6 and RCP4.5, but encapsulates the DSL projections from all five GCMs in RCP8.5. The emulator fails to capture multidecadal variability in DSL, a limitation which would be expected because the emulator is constructed based on the pattern scaling approach.

YUAN AND KOPP 10.1029/2020MS002323 To compare with the DSL projections derived from two-layer pattern scaling, we also produced the DSL projections based on univariate pattern scaling following the same procedure (Figures S12 and S14). The median and 17th-83rd range of DSL projections derived from univariate pattern scaling are similar in patterns (Figure S12). However, differences of both median and spreads of the DSL projections between univariate and two-layer pattern scaling vary across regions, especially over high latitudes in high-emission scenarios (Figure S13). Specifically, compared to the univariate emulator, the median of DSL projected by the two-layer emulator is lower over the Pacific and Indian Ocean, and higher over Arctic, Atlantic and Southern Ocean in the period of 2081-2100, but the opposite in the period of 2271-2290. In addition, the 17th-83rd range of the DSL projected by the two-layer emulator is wider than that by the univariate emulator over the Arctic. The area-weighted increases in DSL spread over the Arctic are 0.02 m in RCP2.6, 0.03 m in RCP4.5, and 0.04 m in RCP8.5 during 2081-2100, and are 0.006 m in RCP2.6, 0.008 m in RCP4.5, 0.009 m in RCP8.5 during 2271-2290. For the grid cell near Philippines, despite the shift in median, the distributions of DSL projections derived from univariate pattern scaling exhibit a different shape from that derived from two-layer pattern scaling (Figure 6). In 2290, the 90% range of DSL projection from the univariate emulator is close to that from the two-layer emulator, slightly narrower by 0.03 m for RCP2.6, 0.02 m for RCP4.5, and 0.01 m for SSP3-7.0 and RCP8.5. The two-layer pattern scaling leads to not only a shift of the distribution but also a different shape of the distribution of the DSL projections, compared to that derived from univariate pattern scaling. There are more DSL projections simulated by GCMs encapsulated within the 90% range of the probabilistic ensemble of DSL projections derived from two-layer pattern scaling than that by the univariate pattern scaling (Figure 6). For the grid cell near the NYC, the shape of DSL distributions derived from univariate pattern scaling is similar to that derived from two-layer pattern scaling, except the spread of the DSL projections derived from univariate pattern scaling is slightly narrower than that derived from the two-layer pattern scaling, with the 90% range of DSL projection in 2290 narrower by 0.03 m for RCP2.6, 0.01 m for RCP4.5, and 0.03 m for SSP3-7.0 and RCP8.5 (Figure 6). The greatest difference is in RCP 2.6, where the difference in the 90% range is by far the largest compared to the overall range.

YUAN AND KOPP 10.1029/2020MS002323

13 of 17

We have developed a probabilistic ensemble of DSL projections through 2300 using a novel two-layer emulator. Replacing the climate module in the FaIR simple climate model with a two-layer energy-balance model, we developed FaIR-2LM, which produces projections of global average temperature in the wellmixed upper layer (T) for rapid responses to radiative forcing, and in the deep ocean layer (T 0 ) for delayed responses. Calibrated by the parameters for each GCMs, the GSAT (Figure S2) and GMTSLR (Figure S3) emulated by FaIR-2LM generally follow that from the corresponding GCM, with RMSE <0.43 K for GSAT and <0.05 m for GMTSLR. A two-layer pattern scaling based on surface and deep-ocean temperature is used to project DSL. During the period 2271-2290, for instance, the DSL predicted by the two-layer pattern scaling are closer to the DSL simulated by the corresponding GCM than that predicted by the univariate pattern scaling (Figure 2). At two selected grid cells (near the coast of Philippines and NYC), the time evolution of DSL projections predicted by the two-layer emulator more accurately reflects GCM behavior and captures non-linearities and non-stationarity in the relationship between DSL and global-mean warming, comparing with that predicted by the univariate technique (Figure 3).

By perturbating the key parameters, FaIR-2LM allows emulation of projected global-mean surface and deep-ocean temperature pairs and GMT-SLR for emissions scenarios (e.g., SSP3-7.0; Figures 4 and 5) beyond those run by the GCMs to which it is calibrated. Compared with the likely ranges assessed by AR5 in the RCP 2.6, 4.5 and 8.5, the FaIR-2LM performs well in emulating the GSAT spread (Table 2 and Figure S5). By 2300, the ensembles of GSAT and GMTSLR estimated by FaIR-2LM have a slightly lower median and a slightly wider 90% range than the estimations by M. D. Palmer et al. (2018), likely because we use the uncertainty of ECS from AR5, which has a larger range than that estimated by CMIP5 multi-model ensemble.

We produce probabilistic ensembles of DSL projections for four different emissions scenarios. Characteristics of median DSL projections during 2271-2290 include increases in DSL along most of the coast around the Pacific and Indian Oceans and a decrease in DSL over the Southern Ocean in all four scenarios, as well as increased DSL over the Arctic and along the North Atlantic Current in moderate to high emissions scenarios (Figure 5). The 66% range (17th-83rd percentile) of uncertainties are small over the middle and low latitudes, and are relatively large over the Southern Ocean, Arctic and North Atlantic, where the simulations of GCMs are diverse due to the challenges of capturing the complex physical processes, such as deep water formation in the subpolar Atlantic, the Antarctic circumpolar current, and ice-albedo feedback in polar regions

YUAN AND KOPP 10.1029/2020MS002323 14 of 17 Scenarios Western Pacific North Atlantic Univariate Two-layer Reduction Univariate Two-layer Reduction RCP2.6 0.01604 0.01192 25.7% 0.027 0.022 18.5% RCP4.5 0.01782 0.01152 35.4% 0.0246 0.0206 16.3% RCP8.5 0.02126 0.01422 33.1% 0.0305 0.0264 13.4%

Note. The averaged RMSEs and the reduction of RMSE from univariate pattern scaling approach to two-layer pattern scaling are calculated for the three RCP scenarios, respectively. Table 3 Comparison of the Distributions of GSAT Anomaly (Relative to 1986-2005) Projected by FaIR-2LM With the Distributions of Global-Mean Surface Temperature Assessed by AR5 (Collins et al., 2013) in RCP 2.6,RCP 4.5, Columns, units: C) (Flato et al., 2013;Landerer et al., 2014;Wang et al., 2014). The ensemble of DSL projections also allows us to examine the trajectories of the DSL projections and their uncertainties at specific locations (Figure 6). At selected locations in the North Atlantic and Western Pacific, the 90% range of DSL spread generally encapsulates the time series of DSL changes relative to the baseline period from the 5 GCMs.

The two-layer emulator provides a useful tool to explore the uncertainty of DSL projections over multiple centuries with computational resources that are much less than a GCM requires. It can be calibrated to match assessments of key values like the equilibrium climate sensitivity, and allows the flexibility of simulating forcing conditions intermediate between the RCPs as the patterns are common for different scenarios. However, we should note that the errors between the DSL predicted by two-layer emulator and DSL simulated by the corresponding GCMs are small in middle and low latitudes but relatively large in high latitudes (e.g., the Southern Ocean, Arctic, and subpolar Atlantic). In addition, the two-layer emulator cannot explicitly resolve the non-linear mesoscale effects of the ocean current due to the coarse resolutions of the CMIP5 GCMs that the two-layer pattern scaling relies on. Comparing with the predicted DSL derived from univariate approach, the improvement of using the two-layer approach on predicting DSL has a similar magnitude with the uncertainty of DSL projections over the middle and low latitudes in RCP2.6 and RCP4.5 scenarios during the period of 2271-2290. But the improvement is limited comparing to the uncertainty of DSL projections in RCP8.5. As the non-linear responses of DSL are more obvious in the RCP2.6 and RCP4.5 than in RCP8.5, the notable improvement of using the two-layer approach over the middle-and low-latitudes in the RCP2.6 and RCP4.5 highlight the advantage on improving the DSL projections in these two scenarios.

Loss of land ice (e.g., Greenland Ice Sheet and Antarctic Ice Sheet) is an important contributor not only to GMSL but also to RSL. Glacio-isostatic adjustment (GIA) caused by changes in ice sheet mass loading induces local vertical land motion and associated changes in local sea level. Mass loss of an ice sheet also reduces the gravitational attraction that pulls sea water toward it, causing water to migrate away with a distinct spatial pattern, or "fingerprint", of sea level change to the global ocean (Mitrovica et al., 2009). Freshwater flux from a melting ice sheet may drastically alter the salinity profile of the near-by ocean, bringing about complex feedbacks involving near-surface ocean stratification, sea ice formation, and corresponding changes in surface temperature, winds, and ocean currents (Bronselaer et al., 2018;Sadai et al., 2020)-each process could have its impact in RSL. Among these factors, only freshwater flux from the ice sheet can significantly affect both global climate and DSL (e.g., Golledge et al., 2019). Despite the importance of polar ice sheets, their contributions to DSL are not included in the current generation coupled climate models.

Our study, relying on outputs from climate models participating the CMIP5 project, thus cannot take into account the effects of evolving ice sheets on DSL. In more comprehensive analyses, the effect of land-ice loss should be considered.