3.2.1 Default configuration: local temperature anomalies as a function of global temperature anomalies

The MESMER emulator was introduced in Beusch et al. (2020a) to emulate multi-ESM initial-condition ensembles for a specific emission pathway in a spatially resolved manner. This first MESMER configuration has thus far been successfully applied to emulate temperature anomalies for the highest-emission scenarios of both the CMIP5 (Beusch et al., 2020a) and the CMIP6 (Beusch et al., 2020b) ensembles.

MESMER is an ESM-specific emulator that needs to be calibrated for each emulated ESM separately to capture the unique characteristics of that ESM, both in terms of local forced warming and local variability around the forced warming (Beusch et al., 2020a). Within MESMER, local temperature anomalies ΔT for a specific climate model “m” at every point in space “s” and time “t” are emulated as follows:

The local temperature anomaly is a direct function of global mean temperature change ΔT^glob and a spatiotemporally correlated noise term η. At every point in space, a multiple linear regression is employed to relate ΔT^glob information to local temperature anomalies. Forced global temperature change ΔT^glob,fr and internal global temperature variability ΔT^glob,iv serve as predictors. The associated regression coefficients are β^fr and β^iv. β^intercept constitutes the intercept term. ΔT^glob,fr is obtained by first applying locally weighted scatterplot smoothing (LOWESS) to the ΔT^glob time series and subsequently adding volcanic spikes to the smooth forced temperature time series via a linear regression of the global temperature anomaly residuals to stratospheric aerosol optical depth. The remaining internal global variability ΔT^glob,iv and the residual local variability η are both modeled as autoregressive processes. For η, spatially correlated innovation terms are drawn to account for local to regional cross-correlations between grid points. Beusch et al. (2020a) describe the full algorithm and additionally carry out a thorough evaluation of the resulting emulations. In short, their evaluation highlights that MESMER reliably emulates grid-point-level forced warming and internal climate variability. However, they also show that MESMER's emulations are increasingly under-dispersive for larger regional averages since MESMER's local residual variability module reduces the magnitude of the empirical covariance between grid cells as a function of their geographical distance. At the global land level, MESMER emulations become reliable again, due to the global-scale variability captured through the ${\mathit{\beta }}_{\mathrm{m},\mathrm{s}}^{\mathrm{iv}}\cdot \mathrm{\Delta }{T}_{\mathrm{m},\mathrm{t}}^{\mathrm{glob},\mathrm{iv}}$ term.

In this work, MESMER is for the first time trained on and applied to the full range of emission scenarios covered within the CMIP6 ScenarioMIP projections instead of a single emission pathway. The main additional assumption when extending MESMER from emulating a single emission pathway to emulating a range of scenarios is that the estimated parameters are scenario-independent. This assumption of universal scenario applicability of the calibrated parameters is also regularly made in the well-established pattern scaling literature (Tebaldi and Arblaster, 2014), although care is nonetheless required, particularly for strong mitigation scenarios (Goodwin et al., 2020).

To obtain robust MESMER parameter estimates for each ESM, MESMER is trained on all available ensemble members of each available scenario and equal weight is given to each scenario. The scenarios employed for training are – subject to availability – SSP5-8.5, SSP3-7.0, SSP4-6.0, SSP2-4.5, SSP5-3.4-over, SSP4-3.4, SSP1-2.6, SSP1-1.9, and the historical scenario. Note that we consider the historical time period as its own scenario during the MESMER training. Since projections for different SSPs usually branch from the same historical members, the historical years would receive more weight if the historical and the future time period were concatenated during training.

Many ESMs have a different number of ensemble members (n_e) for each scenario, and the scenarios (historical and projections) have a different number of time steps (n_t), resulting in a different number of samples ($n_{\mathrm{s}}=n_{\mathrm{e}}\cdot n_{\mathrm{t}}$) per scenario. Therefore, when estimating the coefficients of the multiple linear regression, all simulations from all scenarios are pooled and each sample (e.g., the temperature anomaly for a specific location, year, and scenario) is weighted by $\mathrm{1}/n_{\mathrm{s}}$, i.e., by the inverse number of samples available for the respective scenario. The same approach is applied when estimating the empirical spatial covariance matrix of the residual local variability.

For the autoregressive processes, the coefficients are determined for each ensemble member individually and are subsequently averaged for each scenario before averaging across all scenarios. For the autoregressive process describing the global variability, the order of the process must additionally be selected before the parameters are fit. The order is first chosen for each ensemble member individually based on the Bayesian information criterion and then the median order is identified for each scenario. Lastly, the median order over all scenarios is selected to describe the global variability of the ESM at hand.

The ESM-specific forced global warming time series for each scenario is obtained by using data solely from that scenario, since the forced global warming module relies on a simple statistical smoothing and volcanic spikes approach (see Beusch et al., 2020a, for details). To facilitate the task of the forced global warming module, it operates directly on the average global warming across all ensemble members for a given ESM and scenario rather than on individual ensemble members.

3.2.2 Additional predictors configuration: local temperature anomalies as a function of global temperature anomalies and global ocean heat uptake

To account for potential nonlinearities in the climate system (Mitchell, 2003) and for changing forced response warming patterns when moving from a transient to an equilibrium climate (King et al., 2020), we introduce two additional predictors into the grid-point-level multiple linear regressions of MESMER (Eq. 1): (i) squared forced global temperature change (ΔT^glob,fr)² is used to represent nonlinear feedbacks, and (ii) the forced trend in global ocean heat uptake ΔOHU^glob,fr is used as a proxy for how close the climate system is to equilibrium. ΔOHU^glob,fr is obtained with the same LOWESS approach used to derive ΔT^glob,fr, but no volcanic spikes are added as they are already accounted for in the global temperature response. Hence, the grid-point-level temperature anomalies are emulated as follows:

where β^fr2 and β^fr,ohu are the newly introduced regression coefficients.

3.2.3 Evaluating both MESMER configurations

Before continuing with the methods, we evaluate the two proposed MESMER configurations. Specifically, we assess whether the additional predictors configuration offers considerable benefits over the default configuration. When analyzing climate information, it is often helpful to distinguish between forced response and natural variability. Here, we evaluate both MESMER configurations in terms of their performance in emulating single ESMs with respect to these two quantities. The two can be separated in a straightforward fashion within the additive MESMER framework (Eqs. 1 and 2).

The forced local temperature anomaly ΔT^fr of the default MESMER configuration consists of the following terms:

In the additional predictors configuration of MESMER, the contribution of the additional predictors enters the forced local temperature anomaly too, resulting in the following equation:

Figure 3Latitudinally averaged grid-point-level mean of emulated local forced warming and performance with respect to ESM simulations for both MESMER configurations over the last 30 years of each scenario. The error shown here represents the absolute deviation between the mean of the emulated forced warming and the mean temperature anomaly value across all available ESM initial-condition members for the scenario at hand and is normalized by the absolute value of the mean of the emulated forced warming.

Download

Figure 3 visualizes the latitudinally averaged mean of the local forced warming over the last 30 years for each individual ESM and scenario for both MESMER configurations. For a given scenario, the ESMs' latitudinal patterns exhibit considerable differences. For example, the magnitude of the Arctic amplification (Serreze and Barry, 2011) clearly differs between the ESMs. These regional differences are additionally visualized in maps for select high- and low-emission scenarios for each ESM in the Supplement (Figs. S1–S5). Overall, the patterns of forced warming are very similar in the default and additional predictors configurations (Fig. 3) and inter-ESM differences are considerably larger than MESMER inter-configuration differences.

Nevertheless, with the additional predictors configuration, consistent but small improvements are achieved for most ESMs (Fig. 3). The magnitude of the improvements depends on the ESM, the scenario, and the latitude. The improvements mostly occur in the highest-emission scenario (SSP5-8.5) and in the strong mitigation scenarios (SSP5-3.4-over, SSP1-2.6, and SSP1-1.9). This is an expected consequence of the additional predictors targeting nonlinearities and the transition to an equilibrium climate.

Overall, the performance is excellent for both MESMER configurations (Fig. 3): the latitudinally averaged absolute error with respect to the local warming signal rarely reaches 10 % in the future projections. The SSP1 scenarios have a tendency for higher percentage errors because of their lower overall forced warming signal. This is amplified in the historical period, since its low regional warming leads to small absolute deviations, translating into large percentage errors. In the Supplement (Figs. S1–S5), error maps can be found that further highlight regional differences in the errors for select high- and low-emission scenarios.

Sensitivity experiments are additionally carried out to quantify the impact on emulation performance when using a reduced number of scenarios during training. For MESMER's default configuration, a generally very similar, although slightly reduced, performance is achieved by only training on a single high-emission (SSP5-8.5) and low-emission (SSP1-2.6) scenario (Fig. S6). If only a single future scenario and the historical time period are used for training, a high-emission future scenario results in a better overall performance. This is because training on a low-emission scenario requires MESMER to extrapolate to warm climates when emulating higher-emission scenarios. Overall, the SSP5-8.5-trained default MESMER configuration often performs comparably to MESMER trained on all scenarios but experiences minor performance reductions in strong mitigation scenarios whose climate moves towards equilibrium conditions. The SSP1-2.6-trained emulator, on the other hand, struggles to emulate high-emission scenarios. Hence, for MESMER's default configuration even a single (high-emission) scenario largely suffices to emulate a wide range of different emission pathways as long as their emissions do not exceed the ones used for training. The best emulation performance across all scenarios is, however, achieved if both high- and low-emission scenarios are included in the training.

If the additional predictors configuration of MESMER is applied instead, it is more important that at least one high-emission scenario and one low-emission scenario are available during training because training on a single emission pathway leads to overfitting on the scenario type and thus poor emulation skill in other scenarios (Fig. S7). This mainly occurs because the cross-correlations between the three predictors ($\mathrm{\Delta }{T}_{\mathrm{m},\mathrm{t}}^{\mathrm{glob},\mathrm{fr}}$,$(\mathrm{\Delta }{T}_{\mathrm{m},\mathrm{t}}^{\mathrm{glob},\mathrm{fr}}{)}^{\mathrm{2}}$, and $\mathrm{\Delta }{\mathrm{OHU}}_{\mathrm{m},\mathrm{t}}^{\mathrm{glob},\mathrm{fr}}$) are fundamentally different in the high-emission transient climate change scenarios compared to the low-emission equilibrium-approaching scenarios.

After considering the forced response in detail, we now turn our attention to the local temperature variability for MESMER trained on all available scenarios. The ESM-specific emulated local temperature variability around the forced warming consists of the combination between the local response to the global variability and the residual local variability for both MESMER configurations:

Figure 4Latitudinally averaged grid-point-level standard deviation of emulated local variability and performance with respect to ESM simulations for both MESMER configurations over the full scenario period. The standard deviation of the emulated local variability is computed based on 600 variability emulations for each ESM. To obtain local variability from ESM simulations, the emulated local forced warming is subtracted from every ESM simulation. Subsequently, the standard deviation of these estimates of ESM simulations' local variability is computed while employing all ESM initial-condition ensemble members available for the scenario at hand. The error shown here represents the absolute deviation between the standard deviation of the emulated variability and the standard deviation of the ESMs' variability normalized by the standard deviation of the emulated variability.

Download

Figure 4 shows the latitudinally averaged standard deviation of the emulated local variability analogously to Fig. 3 for each ESM and scenario available for that ESM. However, the standard deviations are computed over the full scenario length instead of only over the last 30 years to obtain more robust estimates. Additionally, the standard deviation of the emulated variability is identical for every scenario since no scenario dependence is integrated in the variability emulation. Differences in the error between different scenarios for individual ESMs hence solely occur because the ESMs' variability differs from simulation to simulation.

Generally, the emulated variability is smallest in the tropics and largest in the northern high latitudes, but considerable inter-ESM differences in the magnitude and spatial distribution exist (Fig. 4). In the Supplement Figs. S8–S12, spatially resolved maps of the standard deviation of the emulated local variability are shown that further highlight region-specific characteristics. MESMER inter-configuration differences are minimal in terms of improvements with respect to capturing the characteristics of the local variability of the actual ESM simulations (Fig. 4). However, the spatially resolved maps can help pinpoint differences in the errors of the two configurations in some regions (Figs. S8–S12).

The percentage errors in the emulated local variability (Fig. 4) are generally a bit larger than the forced local warming ones (Fig. 3). In part, the larger errors are caused by the fact that percentage differences are considered and that standard deviation of local variability is mostly a small quantity. This impression is reinforced by the fact that the largest percentage differences occur most often in the low-variability tropics. A possible physical explanation of the remaining errors is revealed by consulting maps of the error in the standard deviation (Figs. S8–S12). For certain ESMs and regions, the standard deviations are overestimated in the high-end scenario but underestimated in the low-end scenario (or the other way round), indicative of changes in local variability across climate states (Olonscheck and Notz, 2017), which are not accounted for in either MESMER configuration. In line with this finding, the highest deviations are generally observed for SSP5-8.5, which experiences the most extreme climatic change and thus also the largest changes in variability (Fig. 4). Given that SSP5-8.5 is designed to represent an unlikely high-risk future (Hausfather and Peters, 2020), one could justify excluding it from the local variability training to further improve the local variability emulations for the other scenarios. However, since the absolute differences in standard deviations are nevertheless rather small, we continue to use all scenarios for training.

This detailed MESMER evaluation reveals that the additional predictors bring systematic but small improvements. The main forced warming signal can already be successfully extracted based on forced global temperature change alone and local variability is generally emulated similarly well in both MESMER configurations. In the following, we thus use the default MESMER configuration instead of the additional predictors configuration. This minimizes the risk of poorly calibrated parameters when only a limited number of scenarios are available for training in which different predictors are strongly correlated.

3.3.1 ESM-specific emulations

When emulating a specific ESM with MAGICC, a single forced global temperature anomaly ΔT^glob,fr time series is obtained for every emission scenario. Here, we employ ESM-specific MAGICC output published as part of the Reduced Complexity Model Intercomparison Project (RCMIP) Phase 1 (Nicholls et al., 2020) for two ESMs, CanESM5 and CNRM-CM6-1. Note that this output was generated with MAGICC v7.1.0-beta but that very similar results would be obtained with MAGICC v7.5.1, which is employed throughout all other parts of our study.

To obtain full global realizations ΔT^glob for these ESMs, the stochastically generated, ESM-specific global variability of MESMER ΔT^glob,iv is added to MAGICC's ESM-specific ΔT^glob,fr time series. For this study, 600 ΔT^glob,iv emulations are created with MESMER for each ESM. Hence, the MAGICC–MESMER ΔT^glob ensemble contains 600 realizations per ESM.

The spatially resolved forced warming fields ΔT^fr are obtained by combining MAGICC's ΔT^glob,fr time series with the associated ESM's local forced response parameters provided by MESMER (Eq. 3). The resulting ΔT^fr field time series is combined with 600 local variability emulations ΔT^iv for that ESM provided by MESMER (Eq. 5), leading to a 600-member ESM-specific MAGICC–MESMER ensemble of ΔT (Eq. 1).

The ESM-specific MAGICC–MESMER ensembles can be regarded as a direct approximation of very large ESM initial-condition ensembles (Deser et al., 2012, 2020), which can be provided at a negligible computational cost for any emission scenario of interest.

3.3.2 Globally constrained probabilistic emulations

To derive trustworthy probabilistic climate projections that thoroughly sample the climate response and natural variability uncertainty space for any emission scenario, observational constraints are often required. At each point of the MAGICC–MESMER emulation chain, an observational constraint could theoretically be introduced.

Several studies employing fundamentally different approaches have all demonstrated that the CMIP6 ESMs which exhibit the strongest future forced global warming are not consistent with observationally constrained forced global warming estimates (Forster et al., 2020; Beusch et al., 2020b; Brunner et al., 2020b; Tokarska et al., 2020; Nicholls et al., 2021c; Ribes et al., 2021). However, in terms of regional forced warming response to the global forced warming, i.e., in terms of regionally averaged ${\mathit{\beta }}_{\mathrm{m},\mathrm{s}}^{\mathrm{fr}}$ (see Eq. 1), most CMIP6 ESMs perform in an observationally consistent manner in most regions (Beusch et al., 2020b). Therefore, the MAGICC–MESMER constrained probabilistic projections are constrained solely at the global level and span the full regional ESM response uncertainty in this study.

The probabilistic MAGICC output used in this study follows the HadCRUT5 (Morice et al., 2021) calibration of MAGICC presented in Nicholls et al. (2021c) and consists of 600 forced global temperature change ΔT^glob,fr time series per scenario. It is assumed that global- and regional-scale performance are sufficiently decoupled to allow for combining MAGICC's probabilistic output with each of MESMER's ESM-specific local parameter sets. This assumption is supported by Beusch et al. (2020b), who demonstrate that there is no direct relation between an ESM's performance skill for global-scale forced warming response to emissions and regional-scale forced warming response to forced global warming. Hence, similar to the ESM-specific emulation approach, a full initial-condition ensemble is created for each combination of the 600 ΔT^glob,fr time series with the 25 ESM-specific parameter sets. For the global temperature realizations ΔT^glob, this means that each of the 600 MAGICC ΔT^glob,fr time series are combined with 600 global variability ΔT^glob,iv realizations of 25 ESMs, resulting in a total number of 9 million ΔT^glob realizations. For the local temperature anomaly field realizations ΔT, each of MAGICC's ΔT^glob,fr time series is combined with the ESM-specific local parameter sets of MESMER (Eq. 1). This results in 600 local forced temperature change ΔT^fr field time series for each ESM-specific parameter set (Eq. 3). Each of these local forced temperature field time series is subsequently combined with 600 local variability ΔT^iv realizations from MESMER (Eq. 5), yielding 600 different initial-condition ensembles with 600 members each (Eq. 1). Thus, 360 000 probabilistic ΔT realizations are contributed from each set of ESM-specific local parameters, and the full probabilistic ensemble contains 9 million emulations.

The complete sampling approach we employ here ensures a thorough sampling of the full uncertainty space but may quickly lead to computer memory issues. Alternatively, a broad – albeit sparser – sampling of the full climate response and natural variability uncertainty space could also be achieved by randomly combining single MAGICC time series with single ESM-specific parameter sets.