A new energy-balance approach to linear filtering for estimating effective radiative forcing from temperature time series

Reliable estimates of historical effective radiative forcing (ERF) are important for understanding the causes of past climate change and for constraining predictions of future warming. This study proposes a new linear-filtering method for estimating historical radiative forcing from time series of global mean surface temperature (GMST), using energy-balance models (EBMs) fitted to GMST from CO2quadrupling general circulation model (GCM) experiments. We show that the response of any k-box EBM can be represented as an ARMA(k, k− 1) (autoregressive moving-average) filter. We show how, by inverting an EBM’s ARMA filter representation, time series of surface temperature may be converted into radiative forcing. The method is illustrated using three-box EBM fits to two recent Earth system models from CMIP5 and CMIP6 (Coupled Model Intercomparison Project). A comparison with published results obtained using the established ERF_trans method, a purely GCM-based approach, shows that our new method gives an ERF time series that closely matches the GCM-based series (correlation of 0.83). Time series of estimated historical ERF are obtained by applying the method to a dataset of historical temperature observations. The results show that there is clear evidence of a significant increase over the historical period with an estimated forcing in 2018 of 1.45±0.504 Wm−2 when derived using the two Earth system models. This method could be used in the future to attribute past climate changes to anthropogenic and natural factors and to help constrain estimates of climate sensitivity.


Introduction
The estimation of historical radiative forcing, a measure of the net change in the energy balance of the climate system in response to an external perturbation, is a matter of strong scientific interest, as evidenced by the dedication of a whole chapter to this topic in the most recent assessment report from the Intergovernmental Panel on Climate Change (IPCC) (Myhre et al., 2013). Historical radiative forcing; equilibrium climate sensitivity (ECS), the long-term increase in global mean surface temperature (GMST) caused by a doubling of atmospheric CO 2 concentration; and transient climate response (TCR), the increase in GMST after 70 years of CO 2 increasing at 1 % per year, are the three main quantities that inform long-range forecasts of global warming. Estimation of historical forcing is therefore of particular relevance to climate policy decision-makers. In this paper we develop a new method for estimating radiative forcing by inverting simple climate models, a method that can also be of benefit to the detection and attribution of climate change.
The term radiative forcing refers to a change in Earth's energy balance relative to some predefined baseline value, usually chosen to represent pre-industrial conditions. In this study we use the effective radiative forcing (ERF), whose definition is given in Myhre et al. (2013): "the change in net [top-of-the-atmosphere] downward radiative flux after allowing for atmospheric temperatures, water vapour and clouds to adjust, but with surface temperature or a portion of surface conditions unchanged." The adjustments of atmospheric temperatures and other variables, included in the definition Published by Copernicus Publications. 92 D. P. Cummins et al.: A new energy-balance approach to linear filtering of ERF, are commonly referred to as "rapid adjustments" and take place over much shorter timescales than resultant changes in surface temperature (Vial et al., 2013;Chung and Soden, 2015;Smith et al., 2018).
It is infeasible to calculate historical radiative forcing from observational data using the raw definition of ERF, as relevant climate variables, in particular top-of-the-atmosphere (TOA) net downward radiative flux, were unobserved for most of the historical period. Techniques have therefore been developed for diagnosing radiative forcing from general circulation model (GCM) experiments. Forster et al. (2016) describe methods, such as ERF_trans, which use GCMs to simulate the historical period with and without emissions of various forcing agents and hence calculate associated changes in Earth's energy balance. Resulting series of estimated forcings are strictly conditional on the climate model used. This approach is computationally expensive, and there is also some unavoidable noise contamination of results due to internal variability in model output. More recently, Andrews and Forster (2020) have combined GCM simulations with instrumental observations of historical GMST and global heat uptake to constrain historical ERF.
Alternative approaches, based on simple climate models, have been used to estimate historical radiative forcing from the observational record (e.g. Tanaka et al., 2009;Urban and Keller, 2010;Padilla et al., 2011;Aldrin et al., 2012;Urban et al., 2014;Johansson et al., 2015;Ljungqvist, 2015). These studies used Bayesian inference (or non-linear Kalman filtering in the case of Padilla et al., 2011) to jointly estimate some or all of the parameters of a simple climate model together with corresponding series of historical forcing. A common aim of these studies is to constrain ECS, with uncertainty in ECS being found dependent on corresponding uncertainty in historical forcing (Tanaka et al., 2009). As well as being computationally cheaper, methods based on simple climate models have the appeal of directly incorporating observational data, as opposed to GCMs whose dependence on the historical record (through parameter tuning) can be more subtle.
The present study is motivated by the potential application of simple climate models and forcing estimation techniques to the detection and attribution problem. Simple climate models have previously been used for detection and attribution of changes in GMST (e.g. Otto et al., 2015;Rypdal, 2015;Haustein et al., 2017). These studies regressed time series of historical temperature observations on predicted temperature series from linear impulse-response models, thus applying a simplified version of the widely used optimal fingerprinting methodology (Hasselmann, 1997;Allen and Tett, 1999;Allen and Stott, 2003), which is more typically applied to high-dimensional gridded datasets of observations and GCM output.
In the context of detection and attribution, surface temperature is an observable proxy for radiative forcing, which is itself not directly observable. The suitability of this proxy for regression analysis is reduced by two artefacts of the cli-mate system's thermal inertia: (i) a delayed temperature response to changes in radiative forcing and (ii) strong temporal autocorrelation in natural temperature variability. Within the framework of linear impulse-response models, instantaneous surface temperature is simply a convolution of previous changes in radiative forcing (Good et al., 2011). We therefore propose that detection and attribution of changes in GMST using simple climate models could be improved by first deconvolving temperature observations to obtain series of estimated radiative forcing and then by performing the traditional regression step on radiative-forcing time series. In this way we might reasonably expect to eliminate (or at least reduce) both the time delay in the regressed series and also the temporal autocorrelation in the regression residuals.
This paper is a proof-of-concept study in which we develop a method for performing the proposed deconvolution of temperature time series using k-box energy-balance models (EBMs). To do this, we exploit a known correspondence between linear ordinary differential equations (ODEs) and discrete-time autoregressive moving-average (ARMA) time series models. Specifically, we use the fact that a system of k first-order linear ODEs has a discrete-time representation as an ARMA(k, k − 1) filter (Spolia and Chander, 1974;Chang et al., 1982). This correspondence has been used before in the context of three-box EBMs by Grieser and Schönwiese (2001) as a convenient means of integrating the continuoustime model over a time series of discrete forcing inputs and by Stern (2005) as a way to enforce energy-balance constraints on estimated parameters of an ARMA model.
Here we propose a novel application of the three-box EBM filter in its inverted form. Li and Jarvis (2009) showed that the three-box model is simply a physically motivated parameterization of a causal, linear input-output system which maps a time series of radiative forcings onto a corresponding time series of surface temperatures. In this paper we show how the k-box model's equivalent ARMA filter representation may be derived, and we describe how, by inverting the ARMA filter, radiative forcing may be obtained from a temperature time series.

The k-box energy-balance model
An energy-balance model (EBM) is a simple climate model where the time evolution of global temperatures is explained by changes in Earth's radiative imbalance. The simplest class of EBM, known as the k-box (or k-layer) model, consists of two components.
Firstly, a linear relation between GMST anomaly T 1 and TOA net downward radiative flux change N determines the equilibrium response to radiative forcing F (Gregory et al., 2004). In its most basic form, where κ 1 denotes the "climate feedback" parameter, which is conventionally denoted by λ. The notation used here is consistent with Fredriksen and Rypdal (2017) and Cummins et al. (2020). Secondly, a system of k vertically stacked boxes recreates the thermal inertia of the ocean mixed layer and deep ocean, determining the characteristic timescales over which the response unfolds (e.g. Gregory, 2000;Held et al., 2010;Geoffroy et al., 2013). For k = 3: where T i and C i denote the temperature and heat capacity of each box and κ i controls the rate of heat transfer between boxes.
To account for temporal variation in the relationship between T 1 and N , which is present in some climate model experiments, Eqs. (1) and (3) may be modified to include a so-called efficacy factor ε (Held et al., 2010). Equation (1) then becomes while Eq. (3) becomes Internal variability in surface temperature, i.e. fluctuations in T 1 not attributable to changes in F, may be modelled by adding a white-noise process ξ (t) ∼ N(0, σ 2 ξ ) to the righthand side of Eq. (2) (Hasselmann, 1976). For a schematic diagram of the EBM described in this section, see Fig. 1.

ARMA filter representation
For a time-invariant linear system (such as the k-box model), the surface temperature response to a general forcing F(t) can be written as where R(v) is the discrete-time impulse-response function.
Since an impulse is the derivative of a step change, the discrete-time impulse-response function is equal to the oncedifferenced response to a step change in forcing (e.g. a CO 2quadrupling experiment). Defining the "backshift" operator Let It can be shown (see Appendix A) that for the k-box model with φ and θ being polynomials in B of degree k and k − 1 respectively. Thus we have i.e. the linear map from F(t) onto T 1 (t) is an ARMA(k, k−1) filter. The representation in Eq. (10) as a ratio of two polynomials gives us the corresponding inverse filter which is ARMA(k − 1, k). The AR (autoregressive) and MA (moving-average) coefficients can be efficiently calculated from the discrete-time impulse response using software for computing Padé approximants such as the "Pade" package in R (Adler, 2015). Note that the above result for k-box EBMs is a special case of a long-standing and more general result: a system of k first-order linear ODEs has a discrete-time representation as an ARMA(k, k −1) filter (Spolia and Chander, 1974). EBM-derived ARMA filters have previously been used in the climate literature in the "forward" direction, i.e. to convert time series of annually averaged forcings into temperatures (Grieser and Schönwiese, 2001) and for the purpose of EBM parameter estimation (Stern, 2005). Here we apply the ARMA filter in the reverse direction, to find the forcing input required by an EBM to yield a given time series of temperatures.
4 Three-box climate model fits Cummins et al. (2020) developed a maximum-likelihood method for estimating parameters of k-box EBMs from abrupt CO 2 -quadrupling GCM experiments. Using their method, three-box models have been fitted to the two most recent Earth system models (ESMs) from the UK Met Office: HadGEM2-ES (Hadley Centre Global Environmental Model) from CMIP5 (Coupled Model Intercomparison Project) and HadGEM3-GC3.1-LL from CMIP6 (see Fig. 2). Table 1 contains maximum-likelihood estimates of model parameters, as well as estimates of equilibrium climate sensitivity (ECS); transient climate response (TCR); characteristic timescales τ 1 , τ 2 , and τ 3 ; and the coefficients of the ARMA polynomials φ(B) and θ (B). The ARMA coefficients were calculated from the discrete-time impulse responses of the fitted models.
Note that ARMA models can alternatively be estimated using Bayesian inference (e.g. Monahan, 1983).

ARMA filter validation
Numerical estimates of the ARMA coefficients enable, in theory, conversion of time series of surface temperatures into corresponding series of radiative forcings. The properties of this proposed temperature-forcing conversion were investigated using CMIP6 historical runs from the HadGEM3-GC3.1-LL climate model.
Model surface temperatures from 1850 to 2014 were averaged annually (January-December), globally and over four ensemble members. Anomalies were calculated by subtracting from the whole series the mean absolute temperature in the first 50 years (1850-1899). The final temperature series was fed into the fitted HadGEM3-GC3.1-LL ARMA filter using the digital filter implementation in the "signal" package in R (signal developers, 2014).
The resulting filtered forcing series was compared (see  a time-mean historical volcanic forcing to the HadGEM3-GC3.1-LL piControl (pre-industrial control), meaning that they report an ERF of around 0.2 W m −2 in 1850. To avoid an offset due to our choice of temperature baseline, we add 0.2 W m −2 to the forcing series computed from the ARMA filter.
Results from the two methods agree strongly, and the filtered forcing series explains a majority of the variance in the ERF_trans series: the coefficient of determination is R 2 = 0.69 (correlation of 0.83). Reassuringly the ARMA filter correctly infers the rate of forcing increase in recent decades. However, there appears to be a systematic disagreement in years with strong negative forcing. Negative forcing spikes associated with volcanic eruptions are smaller when calculated using the ARMA filter. This is especially evident in the case of Krakatoa (1883).
The observed discrepancy in inferred volcanic forcing is not entirely unexpected. The ARMA filter is derived from a three-box EBM which is known to be unable to resolve temperature responses on timescales significantly shorter than 1 year. Because the data used here are annual averages, there is also scope for error due to discretization, as a volcanic eruption might occur earlier or later in a given year. Finally, it may also be the case that the GCM's temperature response to volcanic forcing deviates from the linearity assumption of the ARMA filter.

Filtering the historical temperature record
By applying the box-model ARMA filters to time series of historical surface temperatures, we can obtain series of esti-mated historical forcings. The Cowtan and Way 2.0 (CW2.0) historical temperature series is an updated version of the dataset described in Cowtan and Way (2014), which is a modification of the HadCRUT4 (Hadley Centre Climatic Research Unit Temperature) series (Morice et al., 2012), corrected for coverage bias. Because CW2.0 is a blend of surface air temperatures (SATs) and sea surface temperatures (SSTs), it is not directly comparable with the pure-SAT 4 × CO 2 datasets used to calibrate the EBMs. The CW2.0 temperatures series have therefore been scaled by 1.09 according to Richardson et al. (2016) for use in this study. The HadGEM2-ES and HadGEM3-GC3.1-LL ARMA filters have been applied to the scaled CW2.0 temperatures 2 , yielding estimated historical forcing series for the period 1850-2018 (see Fig. 4).
It can be seen from panel (d) of Fig. 4 that the forcing series generated by the two ARMA filters are very similar. The HadGEM3-GC3.1-LL ARMA filter reports slightly higher forcing in the latter years of the historical period and slightly more extreme negative forcing in years surrounding major volcanic eruptions.
The reconstructed forcing series are very noisy, since the natural variability which contaminates historical surface temperatures is amplified by the ARMA filter. However, unlike noise in the temperature series, the filtered noise is essentially uncorrelated in time. This follows from the fact that the three-box EBM successfully accounts for the thermal inertia (memory) of the system. The white-noise-like properties of natural variability in the filtered forcing series mean that the long-term trend can be extracted using regression techniques. A generalized additive model (GAM) was fitted to the estimated forcing series using the "mgcv" package in R (Wood, 2011). Approximate 95 % confidence intervals (±1.96 standard errors) for the estimated trend indicate a dramatic acceleration in radiative forcing in the second half of the twentieth century, with an estimated forcing increase of 1.45 ± 0.504 W m −2 between 1850 and 2018.
The filtered forcing series and its subsequent decomposition into signal and noise are subject to multiple sources of uncertainty which must be taken into consideration. As well as the (substantial) internal climate variability, there is observational uncertainty in the historical temperature series and parameter uncertainty in the fitted impulse responses. Given a standard assumption of zero-mean errors, the aggregate noise from observational error and internal climate variability should be well accounted for by GAM regression smoothing. Uncertainty in the fitted impulse responses is also of limited concern: since the EBMs were fitted to abrupt CO 2quadrupling experiments, which have a good signal-to-noise ratio, uncertainty in the fitted impulse responses is in fact quite small (see Fig. 2). It should however be noted that, under forcing scenarios less drastic than 4 × CO 2 , fitting of exponential response functions can suffer from ill-conditioning (De Groen and De Moor, 1987;Kaufmann, 2003). Of arguably greater concern is uncertainty due to inter-model variation between GCMs, which can be seen as a measure of confidence in a particular GCM's ability to represent the true climate. Inter-model variation in GCM output can arise through the use of different parameter values and/or model structures in GCMs (Flato et al., 2013). A full multi-model ensemble, perhaps across all the GCMs in CMIP6, would be required to quantify this satisfactorily, as HadGEM2-ES and HadGEM3-GC3.1-LL are clearly not independent samples. Our numerical results should therefore be seen as conditional on those two climate models.
While we argue that the use of a post hoc GAM regression is reasonable for the reasons given above, a more integrated approach to uncertainty quantification in future analyses might be achieved using Bayesian methods. A Bayesian alternative to GAM smoothing of the filtered forcing series is the "latent-force model" approach (Álvarez et al., 2009;Särkkä et al., 2018). In a latent-force model, uncertainty in the unobserved historical forcing input to a system of ODEs is represented in continuous time using a Gaussian process. Another alternative is the use of sequential methods based on the Kalman filter (Kalman, 1960). ARMA models have natu- ral representations as linear Kalman filters (De Jong and Penzer, 2004), and sequential filtering methods can account for uncertainty in both model parameters and unknown inputs simultaneously, by augmenting the system state vector with estimates of the uncertain parameters (Lourens et al., 2012;Yu and Chakravorty, 2015). Similar techniques have previously been used in the climate literature: for example, by An-nan et al. (2005) and Padilla et al. (2011), who used Kalman filtering to sequentially tune parameters of an Earth model of intermediate complexity (EMIC) and a two-box EBM respectively, and by Cohen and Wang (2014), who estimated historical time series of global black carbon emissions.

Summary
A framework has been developed for estimating radiative forcing from time series of surface temperatures using GCMcalibrated k-box EBMs. There is a known correspondence between k-box EBMs and ARMA models. We have shown how, by inverting a k-box EBM's equivalent ARMA filter representation, a convenient mapping from temperatures onto forcings may be obtained. The method has been validated using historical simulations from the HadGEM3-GC3.1-LL climate model and found to generally perform well, with the notable exception being negative forcing due to volcanic eruptions, which was underestimated. Three-box EBMs fitted to the HadGEM2-ES and HadGEM3-GC3.1-LL climate models have been combined with the Cowtan and Way 2.0 temperature dataset to produce estimates of radiative forcing for the historical period (1850-2018). Forcing estimates calculated using the two models' ARMA filters are very similar. A significant increase in radiative forcing over the historical period has been detected at the 95 % level.
This study has developed a method which has the potential to improve the detection and attribution of past temperature changes by comparing patterns of forcings rather than temperatures. Future work will examine how performance of the method generalizes to the other GCMs in CMIP6 and will address in more detail the question of uncertainty quantification. By combining ARMA filter estimates of historical radiative forcing with known observational constraints, it may be possible to further constrain climate sensitivity metrics, such as ECS and TCR, and hence constrain projections of future warming.

Appendix A: ARMA filter derivation
The step response of the k-box EBM is a linear combination of saturating exponentials (Geoffroy et al., 2013;Tsutsui, 2020;Cummins et al., 2020). The discrete-time impulse response R(v) is a linear function of the step response and is a sum of decaying exponentials: where a i , τ i > 0 and r i = e −1/τ i . Then It follows that the operator (B) is a ratio of polynomials θ (B) and φ(B) of degree k − 1 and k respectively. Cases k = 1, 2, and 3 are shown below.
Case k = 1 so φ(B) and θ (B) have respective degrees one and zero; i.e. the system is an ARMA(1, 0) or simply AR(1) filter.
so φ(B) and θ (B) have respective degrees two and one; i.e. the system is an ARMA(2, 1) filter.
Author contributions. All authors contributed to the development of the statistical methodology and the interpretation of the results. DPC performed the data analysis and was responsible for writing the paper. DBS and PAS proofread and edited the paper.
Competing interests. The authors declare that they have no conflict of interest.

Acknowledgements.
We are grateful to Chris E. Forest and the anonymous referees for their useful comments. We thank Timothy Andrews of the Met Office Hadley Centre for providing forcing data for the HadGEM3-GC3.1-LL model, originally calculated in Andrews et al. (2019). We acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modelling groups for producing and making available their model output. For CMIP the US Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. David B. Stephenson wishes to thank Richard E. Chandler for a helpful discussion about the AR representation of an exponential MA response at the 14th International Meeting on Statistical Climatology.
Review statement. This paper was edited by Chris Forest and reviewed by four anonymous referees.