Monte Carlo drift correction – quantifying the drift uncertainty of global climate models

Abstract

Global climate models are susceptible to drift, causing spurious trends in output variables. Drift is often corrected using data from a control simulation. However, internal climate variability within the control simulation introduces uncertainty to the drift correction process. To quantify this drift uncertainty, we develop a probabilistic technique: Monte Carlo drift correction (MCDC). MCDC samples the standard error associated with drift in the control time series. We apply MCDC to an ensemble of global climate models from the Coupled Model Intercomparison Project Phase 6 (CMIP6). We find that drift correction partially addresses a problem related to drift: energy leakage. Nevertheless, the energy balance of several models remains suspect. We quantify the drift uncertainty of global quantities associated with the Earth's energy balance and thermal expansion of the ocean. When correcting drift in a cumulatively integrated energy flux, we find that it is preferable to integrate the flux before correcting the drift: an alternative method would be to correct the bias before integrating the flux, but this alternative method amplifies the drift uncertainty. Assuming that drift is linear likely leads to an underestimation of drift uncertainty. Time series with weak trends may be especially susceptible to drift uncertainty: for historical thermosteric sea level rise since the 1850s, the drift uncertainty can range from 3 to 24 $\mathrm{mm}$, which is of comparable magnitude to the impact of omitting volcanic forcing in control simulations. Derived coefficients – such as the ocean's expansion efficiency of heat – can also be susceptible to drift uncertainty. When evaluating and analysing global climate model data that are susceptible to drift, researchers should consider drift uncertainty.

1 Introduction

Global climate models are susceptible to drift , causing spurious trends in output variables such as thermosteric sea level rise. Drift may also influence derived climate indices . Drift is not caused by transient external forcing . Instead, drift in long-term centennial-scale simulations is caused by two primary factors: insufficient spin-up and model errors .

First, drift may occur due to insufficient spin-up of equilibrium simulations . Following a perturbation to equilibrium forcing, the deep ocean responds slowly, taking millennia to approach a quasi-equilibrium state . Therefore, variables influenced by the deep ocean are especially susceptible to drift . The slow response of the deep ocean is consistent with the physics of the real climate system. This source of drift can be reduced by increasing the spin-up length of equilibrium simulations.

Second, drift may arise due to errors – including biases – in the climate model . Even when climate models are spun up for long periods under equilibrium forcing, modelled variables may be inconsistent with a true quasi-equilibrium state – spurious trends remain . Of particular relevance here is the problem of energy leakage : unphysical sinks and sources of energy contribute to inconsistent biases in the top-of-atmosphere radiative flux and the sea surface heat flux Sect. 4.1;. Such biases indicate defects in the tuning procedure and the parameterisations included in the global climate model .

Fortunately, drift can often be corrected when analysing global climate model data . Correcting drift may also address the problem of energy leakage . Most drift correction methods rely on control simulations . These control simulations are forced by equilibrium forcing consistent with pre-industrial conditions using 1850 as the reference year . Such control simulations provide initial conditions for historical simulations.

A commonly used drift correction method involves fitting a linear trend to the entire control time series of a state variable . This linear trend can then be subtracted from transient time series – such as time series produced by historical or projection simulations – to correct the drift. However, internal climate variability in the control simulation introduces uncertainty to the drift correction process .

In this paper, we quantify the drift uncertainty associated with the drift correction of global climate model data. To do this, we propose a Monte Carlo drift correction (MCDC) framework (Sect. 3.2). We explore three alternative MCDC methods. Using MCDC, we produce drift-corrected time series of excess system energy ($\mathrm{\Delta }E$), excess ocean heat ($\mathrm{\Delta }H$), and thermosteric sea level rise ($\mathrm{\Delta }Z$). Using the drift-corrected time series, we quantify the drift uncertainty associated with $\mathrm{\Delta }E$, $\mathrm{\Delta }H$, and $\mathrm{\Delta }Z$ during the historical period. We also quantify the drift uncertainty associated with the fraction of energy absorbed by the ocean ($\mathit{\eta }$) and the expansion efficiency of heat ($\mathit{ϵ}$). Our overarching aim is to quantify drift uncertainty. To support this aim, we ask two primary research questions.

• Do the three alternative MCDC methods produce similar estimates of drift uncertainty? If not, which method is preferable?

• What is the contribution of drift uncertainty to $\mathrm{\Delta }E$, $\mathrm{\Delta }H$, $\mathrm{\Delta }Z$, $\mathit{\eta }$, and $\mathit{ϵ}$?

2 Data

We use global climate model data produced for the Coupled Model Intercomparison Project Phase 6 CMIP6; and the CMIP6-endorsed Scenario Model Intercomparison Project ScenarioMIP;. CMIP6 includes contributions from many climate modelling groups globally.

We analyse annual mean global variables that shed light on energy balance and thermosteric sea level rise (Appendix A). We derive the variables from annual mean CMIP6 diagnostic variables as follows.

• Net downward top-of-atmosphere radiative flux ($E^{\prime }$) is calculated as $E^{\prime }=\mathrm{rsdt}-\mathrm{rsut}-\mathrm{rlut}$, where $\mathrm{rsdt}$ is the incoming shortwave radiation, $\mathrm{rsut}$ is the outgoing shortwave radiation, and $\mathrm{rlut}$ is the outgoing longwave radiation. Each top-of-atmosphere flux ($\mathrm{rsdt}$, $\mathrm{rsut}$, and $\mathrm{rlut}$) is multiplied by the atmosphere grid cell area, (areacella) then summed globally.

• Net downward sea surface heat flux ($H^{\prime }$) corresponds to the variable $\mathrm{hfds}$, assuming there is no flux correction . One model – MRI-ESM2-0 – has flux correction ($\mathrm{hfcorr}$) data, so we calculate $H^{\prime }$ as $H^{\prime }=\mathrm{hfds}+\mathrm{hfcorr}$ for this model. Each sea surface flux ($\mathrm{hfds}$ and $\mathrm{hfcorr}$) is multiplied by the ocean grid cell area (areacello), then summed globally.

• Thermosteric sea level rise ($\mathrm{\Delta }Z$) corresponds to the variable zostoga, a global mean diagnostic.

Excess system energy ($\mathrm{\Delta }E$) is calculated by integrating $E^{\prime }$ cumulatively (Eq. A1). Excess ocean heat ($\mathrm{\Delta }H$) is calculated by integrating $H^{\prime }$ cumulatively (Eq. A2). We use the 1850s (1850–1859) as the reference period for $\mathrm{\Delta }E=\mathrm{0}$, $\mathrm{\Delta }H=\mathrm{0}$, and $\mathrm{\Delta }Z=\mathrm{0}$. The fraction of energy absorbed by the ocean ($\mathit{\eta }$) is calculated as the linear regression coefficient of $\mathrm{\Delta }H$ versus $\mathrm{\Delta }E$ (Eq. A3). The expansion efficiency of heat ($\mathit{ϵ}$) is calculated as the linear regression coefficient of $\mathrm{\Delta }Z$ versus $\mathrm{\Delta }H$ (Eq. A4). The coefficients $\mathit{\eta }$ and $\mathit{ϵ}$ are both estimated using ordinary least squares (Appendix B).

Our MCDC framework uses the control simulations . These control simulations follow equilibrium forcing representative of the year 1850. The time at which a historical simulation branches from a control simulation is called the branch time. We use the branch-time metadata to assign corresponding dates to each control simulation so that the year 1850 of the control simulation corresponds to the start of the historical simulation.

We apply MCDC to the historical simulations . The historical simulations follow transient forcing for the period 1850–2014. The 1850–2014 transient forcing includes both natural forcing (including volcanic eruptions) and anthropogenic forcing (including greenhouse gas concentrations and aerosol emissions). We also apply MCDC to ScenarioMIP's four Tier 1 scenarios for the period 2015–2100: SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 . These projection scenarios are based on the narratives of the Shared Socioeconomic Pathways SSPs;. The projection simulations are essentially continuations of the historical simulation, beginning where the historical simulation ends. Therefore, when processing the projection time series and applying MCDC, we prepend the historical period (1850–2014) before the projection period (2015–2100) so that the combined time series covers the period 1850–2100. For example, the cumulative integration that produces the $\mathrm{\Delta }E$ SSP5-8.5 time series begins in 1850 (not 2015), and $\mathrm{\Delta }E$ is subsequently referenced to the 1850s decadal mean. We subsequently remove the historical period from the projection time series before we calculate $\mathit{\eta }$ and $\mathit{ϵ}$ to distinguish the projection simulations from one another and from the historical simulation.

Within the CMIP6 ensemble, we search for model variants that have monthly data for the required diagnostic variables ($\mathrm{rsdt}$, $\mathrm{rsut}$, $\mathrm{rlut}$, areacella, $\mathrm{hfds}$, areacello, and zostoga) and scenario (control, historical, SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5). It is common practice to select only one ensemble member per model e.g.. For each model, we use only the first available “r1i1” variant – in practice, this means we use either the “r1i1p1f1” or “r1i1p1f2” variant, depending on the model (Table S1 in the Supplement). These constraints provide an ensemble of 16 models (Table S1).

We focus on one model as an illustrative example: the UK Earth System Model UKESM1;. With a high equilibrium climate sensitivity of 5.4 $\mathrm{K}$ , UKESM1 is a “hot model” that may overestimate future warming . Nevertheless, as a state-of-the-art global climate model, UKESM1 is well-suited for our present purpose of exploring drift uncertainty. Furthermore, UKESM1 has the longest control time series length (1100 years) of any model within the ensemble (Table S1). Three of the figures presented in this paper show results for UKESM1 (Figs. 1, 2, and 4). Nevertheless, our analysis also includes the larger ensemble (Figs. 3, 5, and S1–S3 in the Supplement; Tables 1 and S2–S6).

3 Drift correction methods

3.1 Conventional drift correction using a single estimate of drift

considered whether to use all or only part of the control time series. They recommended using the entire control time series to estimate the drift “to minimize the contamination of the drift estimate by internal variability” p. 8597. To derive a single best estimate of drift, we should use the entire control time series.

also considered different drift correction methods, including linear, quadratic, and cubic drift correction approaches. Quadratic and cubic drift correction both risk overfitting any curvature in the control time series . On the other hand, such curvature may contain valuable information about non-linearity in the drift. Among recent studies focusing on the Earth's energy budget or sea level change, a few consider the possibility that results may be sensitive to alternative linear, quadratic, and/or cubic models of drift e.g.. However, most studies use only a single statistical model of drift: either linear , quadratic e.g., or cubic e.g.. We observe that researchers generally select a statistical model a priori before analysing any data. They do not generally compare the a posteriori performance of alternative statistical models.

3.2 Monte Carlo drift correction (MCDC)

To quantify drift uncertainty, we propose a Monte Carlo drift correction (MCDC) framework. Following , MCDC estimates the drift using the entire control time series. MCDC also samples the standard error associated with the estimated drift. MCDC proceeds as follows.

• Select a statistical model of drift, such as a linear trend or quadratic polynomial. Alternative models of drift correspond to alternative MCDC methods: for example, a linear trend corresponds to linear-method MCDC (see below).

• Using ordinary least squares, fit the statistical model to the entire control time series of a variable of interest (e.g. $\mathrm{\Delta }E$). This provides the best estimate of the drift, corresponding to a conventional drift correction approach.

• Estimate the standard error associated with each parameter of the fitted model. To account for autocorrelation in the residuals, we use a heteroskedasticity and autocorrelation consistent covariance matrix .

• For each parameter, randomly draw samples from a Gaussian distribution with mean equal to the best estimate (step 2) and standard deviation equal to the standard error (step 3). Use these samples to construct drift samples. Here, we draw 1500 drift samples, which is sufficient to quantify drift uncertainty using the 2nd–98th inter-percentile range.

• For each scenario of interest (e.g. historical), subtract each drift sample from the time series to produce drift-corrected time series. Each drift sample will produce a different drift-corrected time series, enabling quantification of drift uncertainty.

Different statistical models of drift (step 1) correspond to different MCDC methods. In this paper, we apply three MCDC methods: integrated-bias-method MCDC, linear-method MCDC, and agnostic-method MCDC.

Integrated-bias-method MCDC can be applied to $\mathrm{\Delta }E$ and $\mathrm{\Delta }H$, both of which are derived by cumulatively integrating a flux ($E^{\prime }$ or $H^{\prime }$). MCDC is first applied to the flux by modelling the “drift” as a constant: in other words, we assume that the flux contains a constant bias. This application of MCDC to the flux could be referred to as Monte Carlo bias correction. Each bias-corrected $E^{\prime }$ (or $H^{\prime }$) time series is subsequently integrated cumulatively to produce an integrated bias-corrected $\mathrm{\Delta }E$ (or $\mathrm{\Delta }H$) time series.

Linear-method MCDC models drift in a state variable (e.g. $\mathrm{\Delta }E$, $\mathrm{\Delta }H$, or $\mathrm{\Delta }Z$) using a linear trend. In other words, we assume that the underlying drift is linear. For $\mathrm{\Delta }E$ and $\mathrm{\Delta }H$, the assumption of linearity is consistent conceptually with the assumption that the underlying flux contains a constant bias. The assumption of linearity will likely contribute to an underestimation of drift uncertainty (Sect. 5.2).

Agnostic-method MCDC relaxes the assumption of linearity. In addition to sampling the uncertainty associated with the parameters of a given statistical model, agnostic-method MCDC also samples the uncertainty associated with the choice between alternative statistical models: we assume that linear, quadratic, and cubic models of drift are equally valid. This corresponds to the common practice of selecting one of these alternative statistical models a priori (Sect. 3.1). For each of the three models of drift, we draw 500 drift samples. We then combine these samples, producing 1500 drift samples in total.

When applying a quadratic or cubic model of drift, the branch time must be known: which part of the control time series parallels the historical time series? For earlier generations of global climate model ensembles, this branch-time metadata may have been either unavailable or unreliable . To address this, proposed an alternative method to identify the branch time. Here, we assume that the branch-time metadata that accompany the CMIP6 data are reliable. When defining and fitting each statistical model of drift, we use the year 1850 as the reference year that corresponds to the origin.

Figure 1

Uncorrected time series and bivariate relationships for the UK Earth System Model (UKESM1). These uncorrected time series have not yet been corrected for drift. Using data from the control and historical simulations, time series are shown for the following global variables: (a) top-of-atmosphere radiative flux ($E^{\prime }$), (b) sea surface heat flux ($H^{\prime }$), (c) excess system energy ($\mathrm{\Delta }E$), (d) excess ocean heat ($\mathrm{\Delta }H$), and (e) thermosteric sea level rise ($\mathrm{\Delta }Z$). Global total $E^{\prime }$ and $H^{\prime }$ are expressed in global mean units of watts per square metre ($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$). In (a) and (b), the legend shows the mean of the entire time series and the standard error of the mean. In (a–e), only the segment of the control time series that corresponds to the historical period (1850–2014) is shown. (The full length of the entire control time series is 1100 years.) For the historical simulation, bivariate relationships are also shown: (f) $\mathrm{\Delta }H$ versus $\mathrm{\Delta }E$ and (g) $\mathrm{\Delta }Z$ versus $\mathrm{\Delta }H$. Interpretation is offered in Sect. 4.1.

[Figure omitted. See PDF]

4 Results

4.1 Uncorrected time series

An example of drift is illustrated in Fig. 1. For UKESM1's control simulation, the top-of-atmosphere radiative flux ($E^{\prime }$) is close to zero (Fig. 1a), consistent with quasi-equilibrium radiative balance at the top of the atmosphere. The UKESM1 team have achieved their stated goal of keeping $E^{\prime }$ close to zero when tuning the model . Therefore, the drift in excess system energy ($\mathrm{\Delta }E$) is small (Fig. 1c).

However, the sea surface heat flux ($H^{\prime }$) has a negative bias (Fig. 1b), showing that the ocean loses heat spuriously. When we integrate the flux cumulatively, the negative bias in $H^{\prime }$ integrates into a negative trend in excess ocean heat ($\mathrm{\Delta }H$; Fig. 1d). This illustrates the connection between bias and drift: a constant bias in a flux will drive a linear trend in the cumulatively integrated flux . In this example, the drift in $\mathrm{\Delta }H$ also dominates the historical time series, obscuring the anthropogenic warming signal during the 20th century (Fig. 1d). Anthropogenic forcing only becomes strong enough to offset the negative bias in $H^{\prime }$ towards the end of the 20th century (Fig. 1b).

If energy is conserved within the modelled climate system, the $\mathrm{\Delta }E$–$\mathrm{\Delta }H$ relationship should be approximately linear (Eq. A3). In contrast, for UKESM1's historical simulation, the drift in $\mathrm{\Delta }H$ drives a strange $\mathrm{\Delta }E$–$\mathrm{\Delta }H$ relationship (Fig. 1f): for most of the historical period, $\mathrm{\Delta }E$ remains close to zero, but $\mathrm{\Delta }H$ decreases. This reveals energy leakage. Furthermore, the inconsistency between $\mathrm{\Delta }E$ and $\mathrm{\Delta }H$ reveals that the energy leakage occurs somewhere between the top of the atmosphere and the sea surface i.e. outside the ocean;.

In this example, thermosteric sea level rise ($\mathrm{\Delta }Z$) also exhibits drift (Fig. 1e). The drift in $\mathrm{\Delta }Z$ obscures the anthropogenic sea level rise signal during the 20th century. The drift will also contaminate future projections. Furthermore, the $\mathrm{\Delta }H$–$\mathrm{\Delta }Z$ relationship (Fig. 1g) is inconsistent with Eq. (A4).

Figure 2

Different drift correction methods applied to excess system energy ($\mathrm{\Delta }E$) using the UKESM1 control and historical simulations. The first row (a–c) shows integrated-bias-method MCDC results: (a) drift samples derived using the integrated bias method, (b) integrated-bias-method drift-corrected control time series, and (c) integrated-bias-method drift-corrected historical time series, plotted alongside the uncorrected time series. The second row (d–f) shows corresponding linear-method MCDC results. The third row (g–i) shows corresponding agnostic-method MCDC results. These drift correction methods are described in Sect. 3.2.

[Figure omitted. See PDF]

4.2

Excess system energy ($\mathrm{\Delta }E$)

Figure 2 illustrates the application of MCDC to UKESM1's $\mathrm{\Delta }E$ time series. For each MCDC method, we produce 1500 drift samples using the uncorrected control time series (Fig. 2a, d, and g). For UKESM1, these drift samples all indicate positive drift during the historical period. We subtract these drift samples from the uncorrected control time series to produce drift-corrected control time series (Fig. 2b, e, and h). In comparison to the uncorrected time series, the trends of these drift-corrected control time series are much closer to zero, as expected. We also subtract the drift samples from the uncorrected historical time series to produce drift-corrected historical time series (Fig. 2c, f, and i).

The spread among the drift samples and the drift-corrected time series indicates drift uncertainty. For UKESM1, the drift uncertainty is largest when integrated-bias-method MCDC is used (Fig. 2a–c). The linear-method drift uncertainty is much smaller (Fig. 2d–f). This demonstrates that integrated-bias-method MCDC and linear-method MCDC produce different estimates of drift uncertainty, even though these two methods have similar underlying assumptions – namely, a constant bias in $E^{\prime }$ drives a linear trend in $\mathrm{\Delta }E$. We ascribe the difference between the two methods to the size of the standard error: even if autocorrelation is accounted for, the standard error of the trend of $\mathrm{\Delta }E$ is smaller than the standard error of the mean of $E^{\prime }$ because integrating cumulatively effectively averages over the substantial inter-annual variability in $E^{\prime }$. Therefore, integrated-bias-method MCDC provides a much weaker constraint on the drift correction. To minimise drift uncertainty, it is preferable to use linear-method MCDC.

However, the underlying drift may be non-linear. If so, a linear model of drift may be inappropriate. Therefore, linear-method MCDC will likely underestimate drift uncertainty. By partially accounting for the possibility of non-linearity, agnostic-method MCDC should provide a more accurate estimate of drift uncertainty. For UKESM1, agnostic-method drift uncertainty is larger than linear-method drift uncertainty yet smaller than integrated-bias drift uncertainty (Fig. 2; Table S2).

Figure 3

Drift-corrected excess system energy ($\mathrm{\Delta }E$) during the historical period and the corresponding drift uncertainty. We calculate $\mathrm{\Delta }E$ as the decadal mean for the 2000s relative to the 1850s. Each panel shows results for one CMIP6 model (Table S1). Within each panel, each box plot shows results for one MCDC method. The central line shows the median, the box shows the interquartile range, the whiskers show the 2nd–98th inter-percentile range, and the dots show outliers beyond the range of the whiskers. Corresponding results for $\mathrm{\Delta }H$ and $\mathrm{\Delta }Z$ are shown in Figs. S1–S2.

[Figure omitted. See PDF]

Table 1

CMIP6 ensemble median and range (minimum–maximum) for different sources of uncertainty. For each drift correction method, drift uncertainty corresponds to the 2nd–98th inter-percentile range of the drift-corrected data. Model uncertainty corresponds to the inter-model range. Scenario uncertainty corresponds to the inter-scenario range. The statistics for $\mathit{\eta }$ and $\mathit{ϵ}$ are based on the 21st century projection simulations (2015–2100). The ensemble statistics shown here correspond to the summary statistics shown in Tables S2–S6. For further details, see Tables S2–S6.

		$\mathrm{\Delta }E$ (2000s; $\mathrm{YJ}$)	$\mathrm{\Delta }H$ (2000s; $\mathrm{YJ}$)	$\mathrm{\Delta }Z$ (2000s; $\mathrm{mm}$)	$\mathit{\eta }$ (unitless)	$\mathit{ϵ}$ ($\mathrm{mm}{\mathrm{YJ}}^{-\mathrm{1}}$)
Drift uncertainty	Int.-bias	0.07 (0.03–0.15)	0.07 (0.03–0.16)		0.03 (0.01–0.07)
	Linear	0.02 (0.00–0.06)	0.02 (0.00–0.06)	2 (0–8)	0.01 (0.00–0.03)	1 (0–3)
	Agnostic	0.09 (0.03–0.17)	0.08 (0.03–0.20)	10 (3–24)	0.06 (0.01–0.14)	7 (1–22)
Other uncertainty	Model	0.56	0.61	64	0.17 (0.16–0.18)	12 (11–13)
	Scenario				0.01 (0.01–0.08)	7 (4–10)

We extend our analysis to other global climate models within the CMIP6 ensemble by quantifying the drift uncertainty associated with $\mathrm{\Delta }E$ during the historical period (2000s relative to 1850s; Fig. 3; Table S2). For all models, the linear-method drift uncertainty is substantially smaller than the integrated-bias-method drift uncertainty. The agnostic-method drift uncertainty is always larger than the linear-method drift uncertainty and is often larger than the integrated-bias-method drift uncertainty. When averaged across the ensemble, agnostic-method MCDC produces the largest estimate of drift uncertainty for $\mathrm{\Delta }E$, with a median of 0.09 $\mathrm{YJ}$ (Table 1). The ensemble maximum drift uncertainty is 0.17 $\mathrm{YJ}$, approximately twice as large as the ensemble median (Table 1; Fig. 3h).

Differences between global climate models also drive differences in $\mathrm{\Delta }E$ (Fig. 3). We refer to the inter-model range as model uncertainty . For $\mathrm{\Delta }E$, the model uncertainty of 0.56 $\mathrm{YJ}$ is more than 3 times as large as the ensemble maximum drift uncertainty (Table 1).

4.3

Excess ocean heat ($\mathrm{\Delta }H$)

For each CMIP6 model, the drift uncertainty associated with $\mathrm{\Delta }H$ (Fig. S1; Table S3) is similar to the drift uncertainty associated with $\mathrm{\Delta }E$. The ensemble statistics are very similar to those of $\mathrm{\Delta }E$ (Table 1). Therefore, we can draw very similar conclusions for the $\mathrm{\Delta }H$ results as we have for the $\mathrm{\Delta }E$ results. This is not surprising because $\mathrm{\Delta }H$ is closely related to $\mathrm{\Delta }E$: both consist of cumulatively integrated fluxes, and $\mathrm{\Delta }H$ should track $\mathrm{\Delta }E$ closely if the fraction of excess energy absorbed by the ocean ($\mathit{\eta }$) is close to 1.0 (Sect. 4.5; Appendix A).

4.4

Thermosteric sea level rise ($\mathrm{\Delta }Z$)

For $\mathrm{\Delta }Z$ over the historical period (2000s relative to 1850s), the agnostic-method drift uncertainty is always larger than the linear-method drift uncertainty, as expected (Fig. S2; Table S4). (Integrated-bias-method MCDC is not applicable to $\mathrm{\Delta }Z$, which does not consist of a cumulatively integrated flux.) The agnostic-method drift uncertainty ranges from 3 to 24 $\mathrm{mm}$, with an ensemble median of 10 $\mathrm{mm}$ (Table 1). This is smaller than the model uncertainty of 64 $\mathrm{mm}$ (Table 1). On the other hand, the drift uncertainty is of comparable magnitude to the impact of omitting volcanic forcing in control simulations: found that neglecting pre-industrial volcanic forcing leads to an underestimate of 5–30 $\mathrm{mm}$ over the period 1850–2000. Furthermore, for models with relatively large drift uncertainty, the drift uncertainty may influence the extent to which the historical time series agrees with estimates of $\mathrm{\Delta }Z$ derived from observations, such as the reconstructions of and . For two CMIP6 models, the agnostic-method drift uncertainty is large enough to obscure the positive $\mathrm{\Delta }Z$ signal during the historical period: the 2nd–98th percentile range includes negative values of $\mathrm{\Delta }Z$ (Fig. S2a and n).

Figure 4

Drift-corrected bivariate relationships and regression coefficients ($\mathit{\eta }$ and $\mathit{ϵ}$) for the UKESM1 simulations. (a–b) Agnostic-method drift-corrected excess ocean heat ($\mathrm{\Delta }H$) versus excess system energy ($\mathrm{\Delta }E$) for (a) the historical simulation and (b) the projection simulations. (c) Estimates of the fraction of excess energy absorbed by the ocean ($\mathit{\eta }$) and the corresponding drift uncertainty. The coefficient $\mathit{\eta }$ is calculated as the linear regression coefficient of drift-corrected $\mathrm{\Delta }H$ versus $\mathrm{\Delta }E$. For each combination of MCDC method and projection scenario, the central line shows the median, the box shows the interquartile range, the whiskers show the 2nd–98th inter-percentile range, and the dots show outliers beyond the range of the whiskers. (d–e) Agnostic-method thermosteric sea level rise ($\mathrm{\Delta }Z$) versus $\mathrm{\Delta }H$ for (d) the historical simulation and (e) the projection simulations. (f) Estimates of the expansion efficiency of heat ($\mathit{ϵ}$) and the corresponding drift uncertainty. The coefficient $\mathit{ϵ}$ is calculated as the linear regression coefficient of drift-corrected $\mathrm{\Delta }Z$ versus $\mathrm{\Delta }H$. Corresponding estimates of $\mathit{\eta }$ and $\mathit{ϵ}$ for other members of the CMIP6 ensemble are shown in Figs. 5 and S3.

[Figure omitted. See PDF]

Figure 5

Drift-corrected estimates of the fraction of excess system energy absorbed by the ocean ($\mathit{\eta }$) and the corresponding drift uncertainty. Each panel shows results for one CMIP6 model. Within each panel, each box plot shows results for one combination of MCDC method and projection scenario. The central line shows the median, the box shows the interquartile range, the whiskers show the 2nd–98th inter-percentile range, and the dots show outliers beyond the range of the whiskers. The horizontal line at $\mathit{\eta }=\mathrm{1.0}$ shows the theoretical maximum value of $\mathit{\eta }$.

[Figure omitted. See PDF]

4.5

Fraction of excess energy absorbed by the ocean ($\mathit{\eta }$)

If energy is conserved and if the fraction of excess energy absorbed by the ocean ($\mathit{\eta }$) is constant, then excess ocean heat ($\mathrm{\Delta }H$) should be a linear function of excess system energy ($\mathrm{\Delta }E$): $\mathrm{\Delta }H=\mathit{\eta }\mathrm{\Delta }E$ (Eq. A3). However, as noted in Sect. 4.1 above, the uncorrected $\mathrm{\Delta }E$–$\mathrm{\Delta }H$ relationship for UKESM1's historical time series is non-linear and reveals energy leakage (Fig. 1f). Can drift correction address this problem? After we apply agnostic-method MCDC, the trend-corrected $\mathrm{\Delta }E$–$\mathrm{\Delta }H$ relationships are linear (Fig. 4a–b). Estimates of the coefficient $\mathit{\eta }$ are positive and less than 1.0 (Fig. 4c). Therefore, for UKESM1, the drift-corrected $\mathrm{\Delta }E$–$\mathrm{\Delta }H$ relationships are generally consistent with Eq. (A3).

When applying MCDC, we produce 1500 drift-corrected time series for each variable, scenario, model, and MCDC method. Therefore, we also calculate 1500 estimates of $\mathit{\eta }$, so we can quantify drift uncertainty. For all three MCDC methods, the drift uncertainty in $\mathit{\eta }$ is relatively large during the historical period (Figs. 4c and 5): the relatively weak global warming signal in the $\mathrm{\Delta }E$ and $\mathrm{\Delta }H$ historical time series leaves them susceptible to drift uncertainty, which in turn influences the $\mathit{\eta }$ estimates. The drift uncertainty in $\mathit{\eta }$ is much smaller during the 21st century projection period (Figs. 4c and 5): global warming drives clear trends in the $\mathrm{\Delta }E$ and $\mathrm{\Delta }H$ projection time series, constraining the $\mathit{\eta }$ estimates. We average across the projection simulations to summarise the drift uncertainty in $\mathit{\eta }$ (Table S5; Table 1). For UKESM1, the agnostic-method drift uncertainty is approximately 0.01 (Table S5). Across the CMIP6 ensemble, the agnostic-method drift uncertainty ranges from 0.01 to 0.14, with a median of 0.06 (Table 1).

The $\mathit{\eta }$ estimates differ slightly between the different projection simulations (Figs. 4c and 5). The differences arise from the response of the climate system to different forcing scenarios. We refer to the inter-scenario range as the scenario uncertainty. For the UKESM1 projection simulations, the scenario uncertainty is approximately 0.01, similar to the drift uncertainty (Table S5). For most models in the CMIP6 ensemble, the scenario uncertainty is similarly small (Fig. 5; Table S5), with an ensemble median of 0.01 (Table 1): in general, $\mathit{\eta }$ does not depend strongly on the scenario. However, for two models, the scenario uncertainty is 0.08 and 0.07, respectively (Table S5; Fig. 5c–d): for these models, $\mathit{\eta }$ is sensitive to the choice of scenario.

Estimates of $\mathit{\eta }$ vary between different models within the ensemble (Fig. 5). The model uncertainty of 0.17 is the largest source of uncertainty (Table 1). For most of the CMIP6 models, the $\mathit{\eta }$ estimates lie in the range 0.9–1.0, consistent with our understanding of the climate system Appendix A;: even if the ocean were to absorb all excess energy entering the climate system, this would lead to a theoretical maximum of 1.0. However, some of the $\mathit{\eta }$ estimates are larger than 1.0. For four models (25 % of the ensemble), the median $\mathit{\eta }$ estimates are consistently greater than 1.0 for all projection scenarios (Fig. 5h, k, l, and m): even after drift correction, the energy balance of these models remains problematic.

4.6

Expansion efficiency of heat ($\mathit{ϵ}$)

As the ocean warms, it expands, leading to thermosteric sea level rise. We expect the relationship between $\mathrm{\Delta }H$ and $\mathrm{\Delta }Z$ to be approximately linear: $\mathrm{\Delta }Z=\mathit{ϵ}\mathrm{\Delta }H$ (Eq. A4). As noted in Sect. 4.1 above, the uncorrected data from UKESM1's historical simulation are inconsistent with Eq. (A4) (Fig. 1g). However, after we apply agnostic-method MCDC, the drift-corrected $\mathrm{\Delta }H$–$\mathrm{\Delta }Z$ relationships are approximately linear (Fig. 4d and e). Drift correction successfully establishes meaningful $\mathrm{\Delta }H$–$\mathrm{\Delta }Z$ relationships that can be interpreted using Eq. (A4).

However, in disagreement with Eq. (A4), the drift-corrected $\mathrm{\Delta }H$–$\mathrm{\Delta }Z$ relationships during the historical period exhibit hysteresis-like behaviour (Fig. 4d). $\mathrm{\Delta }H$ and $\mathrm{\Delta }Z$ are referenced to the 1850s decadal mean, so the $\mathrm{\Delta }H$–$\mathrm{\Delta }Z$ relationships begin near the origin. $\mathrm{\Delta }H$ and $\mathrm{\Delta }Z$ then become negative for much of the historical period, driving the relationships towards the lower left of Fig. 4d. $\mathrm{\Delta }H$ and $\mathrm{\Delta }Z$ become positive later in the historical period, driving the relationships towards the upper right. However, $\mathrm{\Delta }H$ and $\mathrm{\Delta }Z$ do not necessarily both become positive at the same time, so the relationships do not necessarily pass through the origin, as demonstrated by the “max intercept” and “min intercept” lines. It is possible that the $\mathrm{\Delta }E$–$\mathrm{\Delta }H$–$\mathrm{\Delta }Z$ relationships may depend on the sign of the forcing see. However, the hysteresis-like behaviour in Fig. 4d – and also in Fig. 4a – does not reveal systemic hysteresis: for different MCDC samples, the intercept can be either negative or positive. Regardless of the explanations for this hysteresis-like behaviour, we use ordinary least squares with an intercept – not regression through the origin – to estimate the linear regression coefficient (Appendix B).

The linear regression coefficient is $\mathit{ϵ}$, the expansion efficiency of heat. As was the case for $\mathit{\eta }$, the drift uncertainty in $\mathit{ϵ}$ is much larger during the historical period than the 21st century (Figs. 4f and S3) due to the comparatively weak forcing signal during the historical period. When drift uncertainty is considered, the historical period provides only a weak constraint on the value of $\mathit{ϵ}$.

For the UKESM1 projection simulations, the estimates of $\mathit{ϵ}$ depend on the scenario, ranging from approximately 117 $\mathrm{mm}{\mathrm{YJ}}^{-\mathrm{1}}$ (SSP1-2.6) to 126 $\mathrm{mm}{\mathrm{YJ}}^{-\mathrm{1}}$ (SSP5-8.5; Fig. 4f). This scenario uncertainty of 9 $\mathrm{mm}{\mathrm{YJ}}^{-\mathrm{1}}$ is much larger than the agnostic-method drift uncertainty of 1 $\mathrm{mm}{\mathrm{YJ}}^{-\mathrm{1}}$ (Table S6), leading to distinct non-overlapping box plots for the projection scenarios (Fig. 4f). For all models in the CMIP6 ensemble, the $\mathit{ϵ}$ estimates are smallest for SSP1-2.6 and largest for SSP5-8.5 (Fig. S3). This suggests that the thermal expansion of the ocean is sensitive to the forcing scenario and is not merely a linear function of $\mathrm{\Delta }H$ see also. However, the strength of the dependence on scenario varies: across the ensemble, the scenario uncertainty varies from 4 to 10 $\mathrm{mm}{\mathrm{YJ}}^{-\mathrm{1}}$, with a median of 7 $\mathrm{mm}{\mathrm{YJ}}^{-\mathrm{1}}$ (Table 1).

The agnostic-method drift uncertainty varies over a wider range, from 1 to 22 $\mathrm{mm}{\mathrm{YJ}}^{-\mathrm{1}}$, with a median of 7 $\mathrm{mm}{\mathrm{YJ}}^{-\mathrm{1}}$ (Table 1). The linear-method drift uncertainty is much smaller, suggesting that linear-method MCDC underestimates the drift uncertainty. The model uncertainty is 12 $\mathrm{mm}{\mathrm{YJ}}^{-\mathrm{1}}$. When analysing $\mathit{ϵ}$ across the CMIP6 ensemble, drift uncertainty, scenario uncertainty, and model uncertainty all constitute substantial sources of uncertainty.

5 Discussion

5.1 Energy balance

Following and , we have considered the energy budgets of global climate models. Inconsistent relationships between excess system energy ($\mathrm{\Delta }E$) and excess ocean heat ($\mathrm{\Delta }H$; Fig. 1f) are symptomatic of energy leakage outside the ocean domain . In agreement with and , we have shown that drift correction can partially address this problem of energy leakage. Drift-corrected $\mathrm{\Delta }E$–$\mathrm{\Delta }H$ relationships are approximately linear (Fig. 4a and b). For all 16 models analysed here, the coefficient $\mathit{\eta }$ – the fraction of excess energy absorbed by the ocean – is positive (Fig. 5). For most of these models, $\mathit{\eta }$ is also less than 1.0, consistent with energy conservation. However, even after drift correction, the energy balance of several models remains suspect: energy leakage is revealed by $\mathit{\eta }$ estimates greater than 1.0. For four models, the median $\mathit{\eta }$ estimates are greater than 1.0 for all projection scenarios (Fig. 5).

When evaluating the performance of global climate models, the coefficient $\mathit{\eta }$ is a useful diagnostic. This diagnostic could be used alongside other diagnostics relating to energy balance, such as the diagnostics considered by and .

5.2 Drift uncertainty

previously explored the contribution of drift to uncertainty in global climate sensitivity indices. To do this, added drift (via spurious forcing) and inter-annual variability (via noise) to a simple two-timescale response model. In contrast, we have quantified drift uncertainty directly from global climate model data using Monte Carlo drift correction (MCDC).

We have compared alternative methods of correcting drift in cumulatively integrated fluxes (such as $\mathrm{\Delta }E$ and $\mathrm{\Delta }H$). When applying linear-method MCDC, the flux is first integrated cumulatively, and then the spurious trend is quantified and subtracted. When applying integrated-bias-method MCDC, the bias in the flux is first quantified and subtracted, and then the bias-corrected flux is integrated cumulatively. Conceptually, these two alternative approaches are similar: a constant bias in flux (e.g. $E^{\prime }$) will drive a spurious trend in cumulatively integrated flux (e.g. $\mathrm{\Delta }E$) . In practice, however, we find that these two alternative methods produce differing results: the linear-method drift uncertainty is smaller than the integrated-bias-method drift uncertainty because integrating cumulatively effectively averages over the substantial inter-annual variability, resulting in a smaller standard error. In other words, compared with integrated-bias-method MCDC, linear-method MCDC provides a much stronger constraint on the appropriate drift correction to apply. When correcting drift in an integrated variable, it is generally preferable to integrate the variable first and then apply the drift correction afterward. Nevertheless, we recognise that there may be cases for which it makes sense to correct the bias before integrating: for example, when performing an analysis that involves both $E^{\prime }$ and $\mathrm{\Delta }E$, we may prefer to use the integrated bias method for consistency. Researchers should clarify whether they correct drift before or after integration. Although we primarily refer to cumulative integration across time here, such clarity should also extend to other calculations: for example, clarify that they calculate spatial integrals before correcting drift.

When estimating drift uncertainty, our preferred method is agnostic-method MCDC. If the drift is non-linear, then linear-method MCDC will likely underestimate the drift uncertainty. Therefore, agnostic-method MCDC – which assigns equal prior probabilities to linear, quadratic, and cubic models of drift – should provide a more accurate estimate of drift uncertainty. For the variables analysed in this paper, we have found that agnostic-method drift uncertainty is larger than linear-method drift uncertainty, as expected. We emphasise that agnostic-method corresponds to an assumption that linear, quadratic, and cubic models of drift are all equally plausible. This assumption corresponds to the common practice of selecting a statistical model of drift a priori before fitting the model to data (Sect. 3.1).

A possible further step would be to fit and compare alternative statistical models of drift a posteriori using measures such as the Bayesian information criterion. We could then select the best statistical model(s) for a specific time series. When applying this “best-fit” approach, we could also consider additional statistical models of drift, including signal processing filters e.g.. This best-fit approach would be a sensible way to correct drift. Application of the best-fit approach should lead to a reduction in drift uncertainty.

It may be important to quantify drift uncertainty when analysing time series with a relatively weak trend. In particular, historical time series may be susceptible to drift uncertainty. For the historical period, the agnostic-method drift uncertainty in $\mathrm{\Delta }Z$ can be as large as 24 $\mathrm{mm}$, which is of comparable magnitude to the impact of omitting volcanic forcing in control simulations . Therefore, drift uncertainty deserves similar attention to that given to volcanic forcing.

We hypothesise that drift uncertainty may influence the extent to which a historical simulation agrees with an observation-based estimate. Drift uncertainty should be considered when comparing historical simulations with observation-based estimates of thermosteric sea level change such as the reconstruction of. Drift uncertainty should also be considered when using ocean-related variables – such as ocean heat content  – as an emergent constraint for global climate models.

We further hypothesise that drift uncertainty may be large for drift-susceptible time series that exhibit large internal variability, such as regional dynamic sea level . Large variability will lead to a large standard error when fitting the statistical model of drift. Therefore, drift uncertainty should be considered when analysing time series with large internal variability, provided that the time series is known to be susceptible to drift.

On the other hand, drift is often weak for time series unrelated to the deep ocean . Therefore, many analyses – such as multi-model mean hemispheric partitioning of upper-ocean heat – may be insensitive to drift . In such cases, drift correction may be unnecessary and drift uncertainty may be ignored.

In this study, we have focused on global variables. When correcting drift in multivariate regional contexts – such as might be done in preparation for dynamical downscaling – drift correction approaches may be informed by additional considerations to ensure consistency . Furthermore, we have focused on correcting drift in long-term climate simulations, which are dominated by external forcing. In contrast, shorter-term decadal simulations are further complicated by observation-based initialisation; hence, decadal simulations require different drift correction methods . The quantification of drift uncertainty in such contexts remains an open question.

6 Conclusions

We have developed a probabilistic technique: Monte Carlo drift correction (MCDC). MCDC has enabled us to quantify the drift uncertainty of global climate models.

We have also considered a problem related to drift: energy leakage. Although energy leakage is partially addressed by drift correction, the energy balance of some CMIP6 models remains suspect. Following drift correction, a useful diagnostic for model evaluation is provided by the coefficient $\mathit{\eta }$ (the fraction of excess energy absorbed by the ocean): is $\mathit{\eta }\le \mathrm{1}$?

We draw four conclusions about drift correction.

• When correcting drift in cumulatively integrated fluxes (e.g. $\mathrm{\Delta }E$, the excess system energy), it is generally preferable to integrate the fluxes before correcting the drift rather than correcting the bias before integrating the fluxes.

• Linear-method MCDC may underestimate drift uncertainty due to non-linearity in the underlying drift. If we assume that linear, quadratic, and cubic statistical models of drift are equally plausible, then agnostic-method MCDC provides a more accurate estimate of drift uncertainty.

• When analysing $\mathrm{\Delta }Z$ (thermosteric sea level rise) during the historical period, drift uncertainty is relatively large. For the period 1850s–2000s, the agnostic-method drift uncertainty ranges from 3 to 24 $\mathrm{mm}$. This is comparable to the impact of omitting volcanic forcing in control simulations .

• When using 21st century projection scenarios to estimate $\mathit{ϵ}$ (the expansion efficiency of heat), agnostic-method drift uncertainty of 1 to 22 $\mathrm{mm}{\mathrm{YJ}}^{-\mathrm{1}}$ constitutes an important source of uncertainty alongside scenario uncertainty and model uncertainty.

We propose two hypotheses – each with an accompanying question – to be tested in future work. First, drift uncertainty will influence the extent to which a historical simulation agrees with observation-based estimates. In light of this, can ocean-related variables still be used as emergent constraints for the evaluation of global climate models? Second, drift uncertainty will be large for regional variables with large internal variability, such as dynamic sea level change. What are the implications for climate projections?

Finally, we offer one recommendation: when evaluating and analysing data that are prone to drift, researchers should consider the potential influence of drift uncertainty.

Appendix A Equations connecting energy fluxes, excess energy, and thermosteric sea level rise

If the Earth's climate system were in equilibrium, the incoming and outgoing energy fluxes would balance . In reality, the energy fluxes vary due to internal variability, leading to a quasi-equilibrium state that approximates equilibrium when averaged over long timescales. When the Earth's climate system experiences an external forcing – such as that caused by rising concentrations of greenhouse gases – the energy fluxes no longer balance.

A net downward top-of-atmosphere radiative flux ($E^{\prime }$) heats the Earth's climate system , increasing the excess system energy ($\mathrm{\Delta }E$). Most of this excess system energy is transferred to the ocean via a net downward sea surface heat flux ($H^{\prime }$), increasing the ocean heat content  – here, we refer to the change in ocean heat content as the excess ocean heat ($\mathrm{\Delta }H$). As the ocean warms, the water expands due to thermal expansion, causing thermosteric sea level rise $\mathrm{\Delta }Z$;. To a first approximation, the relationships between these variables can be described using relatively simple equations. Following , we describe the equations below.

We begin by noting that the total excess energy that enters the Earth's climate system can be calculated by integrating the global total top-of-atmosphere radiative flux ($E^{\prime }$) cumulatively over time: A1$$\text{Excess system energy}\equiv \mathrm{\Delta }E\equiv \underset{t_{\mathrm{0}}}{\overset{t}{\int }}{E}^{\prime }\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau },$$ where $t_{\mathrm{0}}$ is the reference time. Expressed in global mean units, $E^{\prime }$ is estimated to be 0.47 $\pm$ 0.1 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ for the period 1971–2018 and is increasing .

On annual to decadal timescales and longer, most excess system energy is absorbed by the ocean . The change in ocean heat content can be calculated by integrating the global total sea surface heat flux ($H^{\prime }$) cumulatively over time: A2$$\text{Excess ocean heat}\equiv \mathrm{\Delta }H\equiv \underset{t_{\mathrm{0}}}{\overset{t}{\int }}{H}^{\prime }\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }.$$

A constant geothermal heat flux of approximately 0.1 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ also contributes to ocean heat content. We ignore this geothermal heat flux – which is much smaller than $E^{\prime }$ – to focus on the relationship between $E^{\prime }$ and $H^{\prime }$. The drift-corrected $\mathrm{\Delta }H$ time series should be insensitive to the inclusion or exclusion of the geothermal heat flux: a constant geothermal heat flux would essentially modify the bias in $H^{\prime }$ and the linear drift in $\mathrm{\Delta }H$, both of which are removed by drift correction.

If we assume that the fraction of excess energy absorbed by the ocean ($\mathit{\eta }$) is constant, then the excess ocean heat can be written as a linear function of the excess system energy: A3$$\mathrm{\Delta }H=\mathit{\eta }\mathrm{\Delta }E.$$

For the period 1971–2018, $\mathit{\eta }$ is estimated to be approximately 0.89–0.91 .

As the ocean warms, the water expands, causing thermosteric sea level rise . The thermal expansion of water varies between locations and depths, increasing with temperature, salinity, and pressure . Remarkably, the thermal expansion coefficient is an order of magnitude larger in the warm low-latitude ocean than it is in the cold high-latitude ocean . Nevertheless, in practice, we can use a globally representative coefficient: $\mathit{ϵ}$, the “expansion efficiency of heat” . Global mean thermosteric sea level rise ($\mathrm{\Delta }Z$) can then be written as a linear function of excess ocean heat : A4$$\mathrm{\Delta }Z=\mathit{ϵ}\mathrm{\Delta }H.$$

For the period 1995–2014, $\mathit{ϵ}$ is estimated to be 121 $\pm$ 1 $\mathrm{mm}{\mathrm{YJ}}^{-\mathrm{1}}$ .

Appendix B

Estimating $\mathit{\eta }$ and $\mathit{ϵ}$ using ordinary least squares

Informed by Eqs. (A3) and (A4), we estimate $\mathit{\eta }$ and $\mathit{ϵ}$ using linear regression. Should we use ordinary least squares with an intercept or regression through the origin ? Regression through the origin is controversial yet sometimes justifiable . In this particular case, Eqs. (A3) and (A4) provide strong theoretical support for the application of regression through the origin. We would expect Eqs. (A3) and (A4) to hold across the full range of the climate model time series, including the origin (because $\mathrm{\Delta }E$, $\mathrm{\Delta }H$, and $\mathrm{\Delta }Z$ are all referenced to the 1850s decadal mean).

However, the drift-corrected $\mathrm{\Delta }E$–$\mathrm{\Delta }H$ and $\mathrm{\Delta }H$–$\mathrm{\Delta }Z$ relationships often exhibit hysteresis-like behaviour during the historical period (Figs. 4a and d, especially the “max intercept” and “min intercept” lines): when the time series transition from negative values to positive values, the $\mathrm{\Delta }E$–$\mathrm{\Delta }H$ and $\mathrm{\Delta }H$–$\mathrm{\Delta }Z$ relationships do not necessarily pass through the origin. Such hysteresis-like behaviour during the historical period influences the starting point of the projection time series. This undermines the appropriateness of regression through the origin. Furthermore, we are interested primarily in linear relationships within the range of data for any given scenario. For example, when estimating $\mathit{\eta }$ and $\mathit{ϵ}$ for a projection scenario, we are interested in the relationship over the period 2015–2100. In light of these considerations, we use ordinary least squares with an intercept.

Code and data availability

The CMIP6 global climate model data can be downloaded from the Earth System Grid Federation (ESGF). The analysis code used to produce the figures and tables can be downloaded from 10.5281/zenodo.8219778 .

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-16-6593-2023-supplement.

Author contributions

BSG: conceptualisation, data curation, formal analysis, investigation, methodology, software, visualisation, writing (original draft preparation, review and editing). ZYK: validation, writing (review and editing). DS: validation, writing (review and editing). BPH: funding acquisition, supervision (supporting), writing (review and editing). JD: funding acquisition, supervision (supporting), writing (review and editing). LYC: funding acquisition, project administration, resources, supervision (lead), writing (review and editing).

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP6. We thank the climate modelling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the multiple funding agencies who support CMIP6 and ESGF. We thank Damien Irving, the anonymous referee, and the editor for their helpful comments. This work comprises EOS contribution number 547.

Financial support

This research project is supported by the National Research Foundation of Singapore and the National Environment Agency of Singapore under the National Sea Level Programme Funding Initiative (award no. USS-IF-2020-3).

Review statement

This paper was edited by Sergey Gromov and reviewed by Damien Irving and one anonymous referee.