The Fractional Energy Balance Equation for Climate projections through 2100

We produce climate projections through the 21 century using the fractional energy balance equation (FEBE) which is a generalization of the standard EBE. The FEBE can be derived either from Budyko-Sellers models or phenomenologically by applying the scaling symmetry to energy storage processes. It is easily implemented by changing the integer order of the storage (derivative) term in the EBE to a fractional value near 1/2. 5 The FEBE has two shape parameters: a scaling exponentH and relaxation time τ ; its amplitude parameter is the equilibrium climate sensitivity (ECS). Two additional parameters were needed for the forcing: an aerosol re-calibration factor α to account for the large aerosol uncertainty, and a volcanic intermittency correction exponent ν. A Bayesian framework based on historical temperatures and natural and anthropogenic forcing series was used for parameter estimation. Significantly, the error model was not ad hoc, but was predicted by the model itself: the internal variability response to white noise internal forcing. 10 The 90% Confidence Interval (CI) of the shape parameters wereH = [0.33,0.44] (median=0.38), τ = [2.4,7.0] (median=4.7) years compared to the usual EBEH = 1, and literature values τ typically in the range 2-8 years. We found that aerosols were too strong by an average factor α= [0.2,1.0] (median=0.6) and the volcanic intermittency correction exponent was ν = [0.15,0.41] (median=0.28) compared to standard values α= ν = 1. The overpowered aerosols support a revision of the global modern (2005) aerosol forcing 90% CI to a narrower range [-1.0,-0.2]Wm−2.The key parameter ECS = [1.6, 2.4]K (median = 2.0K) 15 compared with the IPCC AR5 range [1.5, 4.5]K (median = 3.2K). Similarly, we found the transient climate sensitivity (TCR) = [1.2, 1.8]K (median = 1.5K) compared to the AR5 range TCR = [1.0, 2.5]K (median =1.8K). As commonly seen in other observational-based studies, the FEBE values are therefore somewhat lower but still consistent with those in IPCC AR5. Using these parameters we made projections to 2100 using both the Representative Carbon Pathways (RCP) and Shared Socioeconomic Pathways (SSP) scenarios and shown alongside the CMIP5/6 MME. The FEBE hindprojections (1880-2019) 20 closely follow observations (notably during the “hiatus”, 1998-2015). Overall the FEBE were 10-15% lower but due to their smaller uncertainties, their 90% CIs lie completely within the GCM 90% CIs. The FEBE thus complements and supports the GCMs. 1 https://doi.org/10.5194/esd-2020-48 Preprint. Discussion started: 4 August 2020 c © Author(s) 2020. CC BY 4.0 License.

observational-based studies, the FEBE values are therefore somewhat lower but still consistent with those in IPCC AR5.
Using these parameters we made projections to 2100 using both the Representative Carbon Pathways (RCP) and Shared Socioeconomic Pathways (SSP) scenarios and shown alongside the CMIP5/6 MME. The FEBE hindprojections (1880-2019 20 closely follow observations (notably during the "hiatus ", 1998-2015). Overall the FEBE were 10-15% lower but due to their smaller uncertainties, their 90% CIs lie completely within the GCM 90% CIs. The FEBE thus complements and supports the GCMs.

25
The Earth is a complex, heterogenous system with turbulent atmospheric and oceanic processes operating over scales ranging from millimeters up to planetary scales. When considered by time scale, there are three main regimes: the weather, macroweather and climate ( (Lovejoy and Schertzer, 2013), (Lovejoy, 2013)). From dissipation times up until scales of about ten days -the lifetime of planetary structures -fluctuations in the temperature and other atmospheric quantities increase with time scale: this is the weather regime. Beyond this -in macroweather -fluctuations generally decrease with scale: longer and 30 longer averages tend to converge. Eventually, in the industrial epoch at a scale of ≈ 20 years, this is reversed and fluctuations again tend to increase. This marks the beginning of the climate regime, in the pre-industrial epoch the transition is at centuries or millennia and the regime continues up to Milankovitch scales ( (Lovejoy, 2015), (Lovejoy, 2019b)).
A major challenge is determining the Earth's decadal and centennial response to anthropogenic and natural perturbations.
At the moment projection uncertainties -famously exemplified in the range 1.5-4.5K for a CO 2 doubling -are quite large 35 so that for many purposes, including the development of mitigation policies, the development of complementary approaches is desirable. In considering alternatives, perhaps the most important point is that although perturbations to the Earth system can be quite varied, when compared to the mean solar radiation, over the past and future decades, those of interest are of the order of only a few percent. This allows diverse forcings to be conveniently approximated by their equivalent radiative forcings. It also explains why -in spite of their highly nonlinear weather dynamics -that to a good approximation, GCM macroweather 40 and climate responses to external perturbations are typically linear (as quantified for CMIP5 models in (Hébert and Lovejoy, 2018) but with stochastic internal variability.
In order to construct macroweather and climate models, beyond linearity and stochasticity, we require additional model constraints, the classical one being energy balance. Starting with the first Energy Balance Models (EBMs) proposed by (Budyko, 1969) and (Sellers, 1969), EBMs have been extensively used for understanding the climate ( (North, 1975), (North et al., 1981, 45 Review)). In this paper, we will only consider EBMs for the globally averaged temperature that can be obtained by latitudinally averaging Budyko-Sellers models. The resulting "zero dimensional" energy balance equation (EBE) is a first order linear differential equation, it can alternatively be obtained by considering the Earth to be a uniform slab of material ("box") radiatively exchanging heat with outer space. Such box models usually involve at least two boxes and they assume Newton's law of cooling as well as ad hoc assumptions relating surface temperature gradients to the rate of heat exchange. 50 Energy conservation is an important symmetry principle, yet when implemented in box type models, it violates another symmetry: scale invariance. This is because box models are integer ordered differential equations whose response functions (Green's functions) are exponentials. In order to respect the scaling, these "climate response functions" (CRFs) have therefore been postulated to be scaling (power laws). However the use of pure power law CRFs (e.g. (Rypdal, 2012), (Myrvoll-Nilsen et al., 2020) leads to divergences: the "runaway Green's function effect" (Hébert and Lovejoy, 2015) so that if the Earth 55 is perturbed by even an infinitesimal step function forcing, its temperature monotonically increases without ever attaining thermodynamic equilibrium: its equilibrium climate sensitivity (ECS) is infinite. Whereas the classical EBMs conserve energy but violate scaling, the pure scaling CRF models are scaling but violate energy conservation. Such models can only make projections by using forcings that start and then return to zero. (Hébert et al., 2020) proposed taming the divergences by cutting off the power law CRFs at small scales. The resulting model 60 was scaling at long times and when forced by step functions, it reaches thermodynamic equilibrium. With this truncated power law CRF and using Bayesian techniques (Hébert et al., 2020) were able to make climate projections through 2100 with the IPCC RCP scenario forcings that were coherent with the MME 90% confidence interval. Furthermore, using the historical part of each GCM simulation, they were able to accurately reproduce the corresponding GCM climate projection: for the Earth's globally averaged temperature, both models were thus effectively equivalent. The main drawback was that the CRF model 65 truncation was somewhat ad hoc. In addition, due to the truncation, it was only useful at decadal or longer scales.
To make more realistic models, the key issue is energy storage. Storage is a consequence of imbalances in incoming short wave and outgoing long wave radiation and it must be accounted for in applications of the energy balance principle. As pointed out in ( (Lovejoy, 2019b), (Lovejoy, 2019a)) and developed in  it is sufficient that the scaling principle not be applied to the CRF, but rather to the storage term in the EBE. Rather than energy being stored by uniformly heating a box, 70 energy is instead stored in a hierarchy of structures from small to large, each with time constants that are power laws of their sizes. This conceptual shift can be implemented simply by changing the integer order of the storage (derivative) term in the EBE to a fractional value: the Fractional Energy Balance Equation (FEBE). While  derived the FEBE in a phenomenological manner, ( (Lovejoy, 2020b), (Lovejoy, 2020c), (Lovejoy, 2020a)) showed how it could instead be derived from the continuum mechanics heat equation used in the Budyko-Sellers models. Indeed, by extending Budyko-Sellers models 75 from 2D to 3D (i.e. to include the vertical) and imposing the (correct) conductive -radiative surface boundary conditions, one immediately obtains fractional order equations for the surface temperature. In other words, nonclassical fractional equations turn out to be necessary consequences of the standard Budyko-Sellers approach.
To understand the FEBE's key new features, recall that linear differential equations can be solved with Green's functions; in the classical integer ordered case, these are based on exponentials. However in the general case where one or more terms are of 80 fractional order, they are instead based on "generalized exponentials", themselves based on power laws. In the FEBE, there are two distinct power law regimes with a transition at the relaxation time (estimated to be of the order of a few years, see below).
While the low frequency Green's function can be very close to Hébert's truncated power law, the high frequency regime is able to produce internal variability coherent with the observed scaling and fractional Gaussian noise used by , (Del Rio Amador and Lovejoy, 2019), (Del Rio Amador and Lovejoy, 2020)) for skillfully forecasting the stochastic 85 (internal) variability at monthly, seasonal, interannual (macroweather) scales. In short, there are theoretical arguments as well as empirical evidence that the FEBE accurately models the Earth's temperature response to both internal and external forcing over macroweather and climate time scales. This paper has 3 sections: methods and materials, results, and conclusions. In methods and materials section, we introduce the standard EBE and generalize it to the FEBE, describe the radiative forcing, temperature and GCM simulations that will be 90 used, and finally we introduce Bayesian inference for determining the model and forcing parameters. In the results section, we present the probability distribution functions for the parameters, estimate Equilibrium Climate Sensitivity (ECS) and Transient Climate Response (TCR), and finally we produce global projections to 2100 using the Representative Carbon Pathways (RCPs) and Shared Socioeconomic Pathways (SSPs) which are a little cooler than the CMIP5 GCMs Multi-Model Ensemble (MME) and the currently available CMIP6 GCMs with which they are compared. The zero-dimensional FEBE may be written: ( (Lovejoy, 2019a)), ) Where T (t) is the Earth temperature anomaly with respect to a reference temperature 100 (T (−∞) = 0), τ is the relaxation time, λ the climate sensitivity, F (t) the anomalous external radiative forcing, and H the order of the ("Weyl") fractional derivative: (Γ is the Gamma function). If this derivative is integrated by parts and the limit H → 1 is taken, using (T (−∞) = 0), −∞ D H t T = dT dt so that we recover the standard box EBE. The FEBE thus generalizes the EBE and the power law FEBE 105 relaxation time τ generalizes the exponential EBE relaxation time. Physically T /λ corresponds to the linearized energy flux (rate per area) of long wave black body surface emission and since F is the flux of shortwave forcing, we see that the imbalance F − T /λ represents the flux of energy conducted into storage, it is proportional to the derivative term. At any instant, the storage is thus proportional to the power law weighted integral of past imbalances. The exception is the box model with H = 1 where the relationship is instantaneous.

110
If we solve the FEBE using Green's functions, we obtain: Where G δ,H is the impulse (Dirac) response Green's function, for the FEBE it is given by: 125 At high frequencies (t τ ), important for modelling and predicting the internal variability we have: If we consider the response to Gaussian white noise forcing, then G δ,H ∝ t H−1 implies that T (t) is approximately a fractional Gaussian noise (fGn) with statistical scaling exponent H I = H − 1/2 (when forced by a Gaussian white noise, the FEBE response is exactly a fractional Relaxation noise, see (Lovejoy, 2020a) To see if this is compatible with the value estimated from the low frequency response to external forcings consider the low frequency behaviour (t τ ), important for modelling and projecting the multidecadal responses to external forcing: (note Γ(−H) < 0 for 0 < H < 1). In the box, H = 1, case we have exactly G Θ,1 (t) = 1 − e −t/τ so that when H < 1, the 135 exponential approach to thermodynamic equilibrium is replaced by a power law. (Hébert et al., 2020)  We consider natural and anthropogenic sources of external forcing: solar and volcanic, and anthropogenic forcing. We use the approximate carbon dioxide concentration to forcing relationship (Myhre et al., 1998): Where F CO2 is the forcing due to carbon dioxide, ρ is the concentration of carbon dioxide and ρ 0 is the preindustrial concentration of carbon dioxide which we take to be 277ppm.

Greenhouse Gas Forcing
The global climate is warming and most of the observed changes are due to increases in the concentration of anthropogenic greenhouse gases (GHGs) (IPCC, 2013). Future anthropogenic forcing is prescribed in the Representative Concentration Path-150 ways (RCPs), established by the IPCC for CMIP5 simulations: we considered RCP 2.6, RCP 4.5, and RCP 8.5 (Meinshausen et al., 2011b). RCP 6.0 was omitted in this study since fewer CMIP5 modelling groups performed the associated run. In the CMIP6 simulations the anthropogenic forcings are prescribed in the Shared Socioeconomic Pathways (SSPs) we investigate SSP 126 (strong mitigation) and SSP 585 (strong emission) scenarios, the counterparts to the most distinct previous scenarios RCP 2.6 and 8.5.

155
The RCP scenarios are derived from estimates of emissions computed by a set of Integrated Assessment Models (IAM), these emissions are converted to concentrations using the Model for the Assessment of Greenhouse-gas Induced Climate Change (MAGICC, (Meinshausen et al., 2011a)), while for the SSP scenarios the emissions are converted to forcings using the Finite Amplitude Impulse Response (FAIR) model ( (Smith et al., 2018a)). These scenarios will allow us to compare our results from the FEBE with CMIP5/6 simulations.

160
The wide spread of the scenarios allows for the investigation of the consequences of various future policies, from strong mitigation (RCP 2.6, SSP 126) to business as usual (RCP 8.5, SSP 585) shown in fig. 1  585 is a continuously rising radiative forcing pathway, "business as usual", in which the radiative forcing levels by the end of the 21st century at approximately 8.5W m −2 , and most closely follows current emissions.
In this paper we use the forcing due to carbon dioxide equivalent, F CO2 Eq , as the measure of our anthropogenic forcing,

Aerosol Forcing
Aerosols are a strong component of radiative forcing associated with anthropogenic emissions, resulting from a combination of direct and indirect aerosol effects. There exists high uncertainty of the aerosol forcing, arising from a poor understanding 175 of how clouds respond to aerosol perturbations, compared to the fairly well constrained GHG forcing, thus following ( (Padilla et al., 2011), (Hébert et al., 2020)) we introduce the aerosol linear scaling factor α to account for our poor knowledge of aerosol forcing.
We obtained the CMIP5 aerosol forcing by subtracting from the total CO 2 EQ forcing the combined effective radiative forcing from the gases controlled by the Kyoto protocol, F Kyt , and from those controlled under the Montreal protocol, F M tl . F M tl is 180 given in CFC-12 equivalent concentration and we use the relation from (Ramaswamy et al., 2001) to convert this to W m −2 .
The total amount of aerosol forcing in 2005 given at the 90% confidence interval (CI) in the IPCC AR5 is [−1.9, −0.1]W m −2 , but since then attempts have been made to better constrain this value; (Stevens, 2015) demonstrates that an aerosol forcing more negative than −1W m −2 is implausible. Using results from (Murphy et al., 2009)

185
The prescribed CMIP6 SSP aerosol forcing, F Aer SSP , contain contributions from aerosol-radiation interactions and from cloud interactions: F ari and F aci (Smith et al., 2018a). F ari includes the direct radiative effect of aerosols, in addition to rapid adjustments due to changes in the atmospheric temperature, humidity and cloud profile (formerly the semi-direct effect), and is calculated using multi-model results from Aerocom (Myhre et al., 2013). F aci describes how aerosols affect clouds in the radiation budget and is calculated from the aerosol model of (Stevens, 2015), which includes a logarithmic dependence of F aci 190 on sulphates, black carbon and organic carbon emissions. Year

Solar Forcing
The other external forcings considered are solar and volcanic; there exist other natural forcings such as mineral dust and sea salt, but they are small and will be implicitly included with the internal variability. We use the CMIP5 recommendation for solar forcing, F Sol , a reconstruction obtained by regressing sunspot and faculae time series with total solar irradiance (TSI) 195 (Wang et al., 2005), shown in fig. 2. Following (Meinshausen et al., 2011b), solar forcing anomaly is calculated as the change in solar constant over the average value of the two year 11-year solar cycles from 1882 to 1904 divided by 4 (average insolation) and multiplied by 0.7 (representing planetary co-albedo). To extend solar forcing to the future we reproduce solar cycle 23 (the last one prior to 2008) as a proxy for future solar forcing.

200
The volcanic forcings series, F V ol , used in this study was generated from the volcanic optical depths, τ V . Over the 1850 to 2012 period we use the approximate relation: F V ol ≈ −27W m −2 τ V , obtained from the Goddard Institute for Space Science (GISS) website (Sato, 2012). We follow (Hébert et al., 2020), extending the series to 1765 using the optical depth reconstruction (Crowley et al., 2008), and setting volcanic forcing to zero for the future.
It has been shown by (Lewis and Curry, 2015) that volcanic forcing must be scaled down by 40-50% in order to produce a 205 comparable effect on surface temperature, and thus most EBMs linearly scale volcanic forcing. However the amplitude of the volcanic forcing is not the only issue, volcanic forcings are highly intermittent (spiky). The intermittency can be quantified in a multi-fractal framework ( (Lovejoy and Schertzer, 2013), and (Lovejoy and Varotsos, 2016)). Since linear response models do not alter the intermittency, the volcanic series must first be non-linearly transformed before being introduced into a linear response framework. With the effective volcanic forcing F V olν , the volcanic intermittency correction exponent ν and the mean 210 of the whole volcanic series F V ol , we follow (Hébert et al., 2020) using a non-linear relation to change the intermittency so that the transformed signal can be linearly related to the temperature: The normalization is such that the mean is unchanged: The volcanic intermittency correction exponent, ν, required to reduce the intermittency parameter of the volcanic forcing, C 1,F V , to equal the corresponding parameter of the 215 temperature response, C 1,T V , can be calculated theoretically using: where α M F is the multifractality index of the volcanic forcing (Lovejoy and Schertzer, 2013), C 1 is the codimension of the mean. For the volcanic forcing intermittency, C 1,F V ≈ 0.16, the temperature response intermittency, C 1,T V ≈ 0.03, and α M F ≈ 1.5, we find an approximate but plausible theoretical estimate of the volcanic intermittency correction exponent ν ≈ 0.3 220 ((Lovejoy and Schertzer, 2013) table 11.8, (Lovejoy and Varotsos, 2016)).
Working in a linear framework we write the forcing series as the sum of anthropogenic and natural forcings:

Surface Air Temperature Data and CMIP5/6 Simulations
We used five historical records of surface air temperature for our analysis each spanning the period 1880-2019, with median Year -12 is shown alongside two non-linearly transformed versions. Linearly damped by a constant 0.5 coefficient (black), and non-linearly transformed using equation with ν = 0.3 (red). The solar forcing F Sol (orange) has been shifted down by -9 and amplified by a factor of 10 for clarity. Adapted from (Hébert et al., 2020) the Cowtan and Way dataset uses HadCRUT4 as raw data, but address missing data that would lead to bias especially at high latitudes by infilling missing data using an optimal interpolation algorithm (kriging); we use the dataset with land air temperature anomalies interpolated over sea-ice. The GISTEMP dataset combines the Global Historical Climate Network version 3 (GHCNv3) land surface air temperature records with the Extended Reconstructed Sea Surface Temperature version 235 4 (ERSST) along with the temperature dataset from the Scientific Community on Antarctic Research (SCAR) and is compiled by the Goddard Institute for Space Studies; the NOAA National Climate Data Center uses GHCNv3 and ERSST but applies different quality controls and bias adjustments. The final data, BEST, makes use of its own land surface air temperature product along with a modified version of HadSST.
The CMIP5 models selected have monthly historical simulation outputs available over the 1860 to 2005 period along with 240 outputs of scenario runs from 2005 to 2100 for RCP 2.6, RCP 4.5, and RCP 8.5, summarized in table A1. The available CMIP6 model outputs for the historical and SSP scenarios are provided in (Forster et al., 2020), and climate sensitivity of models is taken from (Flynn and Mauritsen, 2020), summarized in table A2.

Bayesian Parameter Estimation
In this section we establish a procedure to estimate the probability distribution associated with the climate sensitivity: λ, model 245 parameters: τ , H and the forcing parameters: α, ν. To estimate them, we relate the forcing to surface air temperature data using the FEBE in a multi-parameter Bayesian technique. To apply Bayesian inference we require temperature observations, a statistical model that relates forcing data to temperature, and prior information about the model parameters (priors). Bayesian inference is chosen due to its ability to better constrain model parameters by using information from different sources including data and models.

250
In this framework, each parameter combination (H, τ for G δ,H and α, ν for F as well as λ) produces a time-dependent forced response which is associated with a likelihood that depends on how well the corresponding model output matches the observational temperature records over the historic period. To see how this works, recall that the FEBE describes the temperature response to the total external deterministic forcing F (t) and internal stochastic forcing γ(t): 255 Where T ext , T int are the responses. Any given set of parameters (H, τ for G δ,H and α, ν for F as well as λ) implies a series T ext (t); and hence when compared to observation temperature series, it implies a series of residuals: The residuals are thus equal to the internal temperature variability, i.e. the response to the internal forcing γ(t). Here we make the usual assumption that γ(t) is a Gaussian white noise so that T res (t) = T int (t) is a fractional Relaxation noise process (fRn, 260 (Lovejoy, 2019a)). However, for scales shorter than the relaxation time τ (of the order of years), the fRn process is very close to a fractional Gaussian noise (fGn) process (due to the approximation G δ,H ≈ G δ,high,H (t), eq. 8). The fGn approximation takes into account the strong power law correlations induced by the fractional derivative term in the FEBE and it is valid except at the low frequencies that only weakly influence the likelihood function. It is already a much more accurate residual model than the standard Auto Regressive order one (AR (1)) assumption which models temporal correlation as a (short range) exponential 265 function. Rather than making an ad hoc assumption about the statistics of the residuals, in our approach the statistics are given by the model itself. It should be noted that if the residuals are fGn, then the uncertainties are larger than for AR(1) (or other exponential decorrelation models) in spite of the fact that our residuals are more realistic.
We therefore estimate the likelihood of any parameter set from the maximum likelihood function that the residuals are an fGn process: The likelihood function (L) corresponds to the probability ("P r") of observing the series T (t) conditioned on the parameters: λ, H, τ, α, ν (right hand side), assuming the residuals are a fGn process with parameter H, and zero mean.

Results
In this section, using Bayes' theorem as described above, we derive probability density functions (PDFs) for the model and forcing parameters of the FEBE from the mean likelihood functions of the five observational datasets. We treat the different observational datasets as dependent due to the use of overlapping raw data, with the differences between series coming partly 295 from the different processing of the raw data by different teams. This corresponds to putting the datasets into a Bayesian framework where each has equal a priori probability: HadCRUTv4, C&W, GISTEMP, NOAAGlobalTemp and BEST.
P r(λ, H, τ, α, ν|T (t)) = 1 Following IPCC methodologies, we report the "very likely" confidence interval at the 90% confidence level throughout this work along with median estimates for the all ensemble spreads. The complete suite of model and forcing parameters and  The second model parameter is the relaxation time scale τ that characterizes the approach to equilibrium. τ is a difficult parameter to determine since from the point of view of parameter estimation, it is inversely correlated with λ: a large τ can be somewhat compensated by a smaller λ and vice versa. Since the one box (EBE) model is the special H = 1 case of the FEBE, we can use box model estimates of τ for our prior distribution. For this purpose, we used the (Geoffroy et al., 2013) two-box 315 relaxation times found by fitting a dozen GCM outputs: 4.1 ± 1.1 years (the fast box τ ). This prior is close to other box model fast relaxation times: τ = 8.5 ± 2.5 years (Schwartz, 2008), τ ≈ 4 years (Held et al., 2010), τ = 4.3 ± 0.6 years (Rypdal and Rypdal, 2014). So as not to be very constraining, we doubled the (Geoffroy et al., 2013) parameter spreads, choosing an N[4,2] prior distribution of τ . With this prior, we obtained the a posteriori median value of 4.7 years and 90% CI of [2.4,7.0] years when using F Aer RCP , and nearly identical results using F Aer SSP . we find both cases that α < 1. The result from two independent aerosol forcing series again shows that the forcing associated with aerosols is still widely uncertain and overpowered, supporting post-AR5 studies that found aerosol forcings simulated by 330 GCMs were unrealistic ( (Zhou and Penner, 2017), (Sato et al., 2018), (Bellouin et al., 2020)), and that aerosol forcing was weaker when climate feedbacks were allowed (Nazarenko et al., 2017).

Volcanic Intermittency Correction Exponent ν
The volcanic intermittency correction exponent ν, given an uninformative uniform prior, ν ∈ [0, 1] (where ν = 0 implies a constant mean forcing and the original series if ν = 1), was found to have a posterior median value of 0.28 with 90% CI 335 of [0.15,0.41] when using F Aer RCP and similar median value 0.28 with 90% CI of [0.16,0.40] when using F Aer SSP . Both contain the theoretically calculated ν within their 90% CI (ν = 0.32). This result confirms that volcanic forcing is generally overpowered since ν = 1 has nearly null probability as seen in figure 4. Thus, the original volcanic series described without the intermittency correction does not reproduce well, within the FEBE model presented, the cooling events observed in instrumental records following eruptions: the cooling would be overestimated.

340
In fig. 5 (left) we compare the total forcing series, F T ot (t) (black), IPCC AR5, eq. (13) where α = ν = 1, with the adjusted forcing series, F T ot (α, ν; t) (blue). During the historical period, the intermittency and strength of the strong volcanic events is greatly reduced and in the recent past the median adjusted forcing series is higher than the unadjusted forcing due to the reduced aerosol forcing strength. This adjusted forcing series consequently contributes to a lower climate sensitivity, presented in the following section, due to the historic negative forcings of volcanoes and aerosols being adjusted to closer match historical 345 observations, eliminating the need for a high climate sensitivity to compensate.

Equilibrium Climate Sensitivity
Climate sensitivity is a key measure used in climate modelling. Two standard types of climate sensitivity are used for intermodel comparisons: Equilibrium Climate Sensitivity (ECS) and Transient Climate Response (TCR). Climate sensitivity de-360 scribes the amount of warming in the atmosphere associated with increases in atmospheric carbon dioxide (CO 2 ).
If atmospheric CO 2 were increased to double pre-industrial concentrations and then held there, the planet would only slowly reach a new thermodynamic equilibrium. This delay is largely because the world's oceans take a long time to heat up in response to the enhanced greenhouse effect. The Equilibrium Climate Sensitivity (ECS) is the amount of warming achieved when the entire climate system reaches 'equilibrium' or the stable temperature response to a doubling of CO 2 , it is equal to λ 365 when measured in units of K/(CO 2 doubling).
The PDF for ECS shown in figure 8 (left), for both aerosols series was found to have a 90% CI of [1.6, 2.4]K and a median value of 2.0K when using F Aer RCP , and median of 1.8K and 90% CI [1.5,2.2] using F Aer SSP . These results are lower than those found in the CMIP5 MME which had a best value of 3.2K, but our 90% CI bounds are more narrow, laying within the CMIP5 MME range of [1.9, 4.5]K. Although when we consider the expanded ECS 90% CI of [1.5, 4.5]K considered in 370 (IPCC, 2013), which takes into account both the CMIP5 MME and historical estimates, we see that the FEBE estimates are wholly within this range and much less uncertain. For the CMIP6 MME which has a 90% CI of [2.0, 5.5]K and mean estimate 3.7K, our best estimate using the corresponding F Aer SSP is slightly below the lower confidence due to the upward shift of ECS estimates seen in CMIP6 models (Zelinka et al., 2020), but again has a more narrow CI.
The climate sensitivity parameter, λ, refers to the equilibrium change in the annual GMST following a unit change in 375 radiative forcing. By the definition of the temperature response to external forcings in eq. 7, the climate sensitivity parameter is the equilibrium climate sensitivity. The two are equivalent to within a constant factor: the number of W m −2 per CO 2 doubling, the standard value being 3.71W m −2 /(CO 2 doubling) IPCC (2013). Its inverse is the climate feedback parameter, the increase in radiation to space per unit of global warming. We find λ to have a median value of 0.56 K(W m −2 ) −1 with 90% CI [0.45,0.67] K(W m −2 ) −1 using F Aer RCP (fig. 7). When using F Aer SSP we find median 0.52K(W m −2 ) −1 and 90% CI 380 [0.43,0.61] K(W m −2 ) −1 . Both on the lower end of the CMIP5 MME climate sensitivity parameter of 1 ± 0.5K(W m −2 ) −1 but within the 90% CI although both estimates are below the CMIP6 MME 90% CI [0.63, 1.50]K(W m −2 ) −1 which has been criticized as being too high ( (Zelinka et al., 2020), (Tokarska et al., 2020), (Flynn and Mauritsen, 2020)).

Transient Climate Response
Conventionally, TCR quantifies the changes that would occur if CO 2 levels increase by 1% (compounded) per year until they 385 double (≈ 70 years). Since the CO 2 forcing is logarithmically dependent on CO 2 concentration, the TCR is then simply the global temperature increase that has occurred at the point in time that a linearly increasing forcing reaches double pre-industrial levels.
The derived PDF for TCR is shown in fig. 8 (middle). Our TCR was found to have a 90% CI of [1.2, 1.8]K with a median of 1.5K when using F Aer RCP , while when using F Aer SSP we find a median 1.4K and 90% CI of [1.1,1.6]K. Both estimates 390 are lower and more constrained, but within the 90% CI given by the CMIP5 MME: a 90% CI of [1.2, 2.4]K and a best value of 1.8K, and by the CMIP6 MME: 90% CI of [1.2, 2.8]K with best value of 2.0K.
The TCR-to-ECS ratio is a non-dimensional measure of the fraction of committed warming already realised after a steady increase in radiative forcing, in this case a doubling of CO 2 , this quantity is generally referred to as realised warming fraction (RWF) ( (Stouffer, 2004), , (Millar et al., 2015)). A model with a low RWF will indicate that global 395 warming may continue for centuries after emissions have stopped. We present the TCR-to-ECS ratio in fig. 8 (right), having a 90% CI [0.70, 0.78] and median 0.73 using F Aer RCP . Similar results are found using F Aer SSP , a median of 0.72 and 90% CI [0.71, 0.79]. The TCR-to-ECS ratio is higher than both generations of MME 90% CI, a consequence of lower ECS and TCR values, and similar uncertainty.
In the next section we show that with a lower and more constrained climate sensitivity parameter (figs. 7 and 8), the adjusted 400 forcing and long memory process of the FEBE ( fig. 5) produces future projections that tend to be cooler than the CMIP5/6 projections, yet within their 90% CI. The associated 90% CI (bars under the axis), the CMIP5 MME 90% CI (dark gray shading), and the CMIP6 MME 90% CI (light gray shading). With the full suite of model and forcing parameters, we now use the FEBE to reconstruct the forced temperature response over 405 the historical period and make projections for the coming century using forcings prescribed by the RCP and SSP scenarios. In this section we present comparisons of the FEBE observation-based projections with those from the GCMs in the CMIP5 MME, and all the currently available CMIP6 GCMs analyzed by (Forster et al., 2020). The CI provided for the MME corresponds to the spread between the different GCMs.
From the multidimensional parameter space we draw from the multinormal distribution (eq. 18) using a Monte Carlo method 410 to generate realizations of the forced temperature response (eq. 3) using a numerical convolution. Once we have our ensemble of projections, we remove the pre-industrial baseline  and calculate the desired confidence intervals of the forced response, the historical reconstructed period (hindprojections) presented in fig. 9 and projections to 2100 in fig. 10.
Throughout 1880-2005 (the overlapping observational and CMIP5 MME period), the reconstructed forced temperature response from the FEBE, using F Aer RCP , and the mean of the CMIP5 MME are close, seen in fig. 9). The spread between 415 the forced temperature response and the five observational datasets decreases as we move closer to the present. Between 1915Between -1960 period the CMIP5 MME is consistently warmer than the historical reconstruction and historical temperature records, although generally by less than 0.05K. The slowdown in global warming during the first decade of the 21st century, termed as 'the pause' or 'the hiatus' ((Kaufmann et al., 2011), (Meehl et al., 2011), (Medhaug et al., 2017)), is tracked closely by the FEBE reconstruction while the CMIP5 MME overshoots (by 0.1K to 0.2K), a well studied divergence between GCMs and 420 observations ( fig. 9 lower plot). This supports (Lovejoy, 2015) who found that the hiatus could be well predicted by a stochastic fGn model (comparable with the present hindprojection) and concluded that the problem was GCM overprojection. . The reconstruction of the forced temperature response over the historical period (1880-2019) using the FEBE, with parameters calibrated using FAer RCP (blue), compared with the CMIP5 MME projection (black) and mean of 5 observational temperature series (red); 90% CI are indicated (shaded). The whole historical period shown (top), along with the 'hiatus/pause' where the CMIP5 MME median is consistently above the temperature records while the FEBE median tracks the observations.

Projections through to 2100
The future pathways only begin to diverge into their respective scenarios roughly two decades after their beginning (RCPs begin in 2000, SSPs begin in 2010), so at first the temperature increase in each case is nearly identical. Further into the future, 425 the warming rate begins to depend more on the specified scenario, the highest being in RCP 8.5/SSP 585 ( fig. 10 bottom) while significantly lower in RCP 2.6/SSP 126 ( fig. 10 top), particularly after about 2050 when the global surface temperature response stabilizes (and declines thereafter). Of particular interest are the low emissions scenarios, RCP 2.6/SSP 126, demonstrating the potential of strong mitigation policies and speculative negative emission technologies where anthropogenic forcing starts decreasing around the mid-2040s. In the CMIP5 MME, the temperature stays below 2K throughout the 21st century, whereas 430 the corresponding median FEBE temperature projection never exceeds 1.5K. In comparison, at 2100, the FEBE projection reaches While the forcing of the (perhaps most realistic) middle scenario, RCP 4.5, stabilizes in the mid 2060s, the temperature projections continue rising throughout the 21st century for both FEBE and the CMIP5 MME ( fig. 10 middle). At 2100 the FEBE and CMIP5 MME projects a temperature reaching 1.9K [1.6,2.2]K, and 2.6K [1.8, 3.2]K respectively. A key point to note is that the FEBE RCP 4.5 projection remains below 2K of warming by 2100, while the CMIP5 MME is well beyond this with SSP scenarios to be warmer than their RCP counterparts.
The discrepancy in aerosol forcing strength ( fig. 1) at 2100 between F Aer RCP and F Aer SSP in the strong mitigation scenario (RCP 2.6/SSP 126) is 0.6W m −2 . With the total GHG forcings being nearly identical, the the SSP 126 scenario is closer to a RCP 3.3 scenario than the RCP 2.6 scenario. Therefore for RCP 2.6, we can roughly attribute 30% of the increased future 455 warming shown in CMIP6 in comparison with CMIP5 (figs. 10, 11) to the strong reduction of future aerosol emissions. For the MME, in RCP 2.6, this corresponds to about 0.5K in 2100. The corresponding values for RCP 8.5 is 7% and 0.6K in 2100.
Although the FEBE projections are consistently about 15% cooler than the CMIP5 MME, the FEBE 90% CI lies entirely within the corresponding CMIP5 CI. Being in complete agreement, the two projections methods mutually support each other and are thus complementary. When compared to CMIP6 projections although most of 90% CIs overlap, the median CMIP6 460 temperatures are by nearly 65% warmer than the corresponding FEBE median. The main reason being their overpowered aerosols ( (Zelinka et al., 2020), (Flynn and Mauritsen, 2020)) and previously mentioned removal of future aerosol forcing.  Figure 11. The forced temperature response projected using the FEBE, with parameters calibrated using FAer SSP (blue), compared with the CMIP6 projection (black); 90% CI are indicated (shaded). The projections until 2100, for SSP 126 (top), and SSP 585 (bottom), are shown.
( (Schurer et al., 2017), (Smith et al., 2018b), (Iseri et al., 2018)) have given evidence that there are important tipping points which could lead to irreversible changes in major ecosystems and the planetary climate if certain threshold in warming are exceeded. Therefore we calculate the probability of temperature exceeding 1.5K and 2.0K for both the FEBE and the CMIP5/6 465 MME shown in figs. 12 and 13.
According to the FEBE for the low emission scenario RCP 2.6, it is unlikely to exceed the 1.5K threshold in 2100 (< 10%) while it is much more likely to exceed this threshold according to CMIP5 MME (probability of 67%). The FEBE has a negligible probability of exceeding 2K while the CMIP5 MME has a 26% chance. While in the SSP 126 scenario, the FEBE peaks at below 50% of exceeding 1.5K and has a negligible probability of exceeding 2K as before; it is nearly inevitable to 470 cross 1.5K threshold according to the CMIP6 MME, while the 2K threshold being exceeded hovers around 60% even under strong mitigation.
In the RCP 4.5 scenario, the probability of the FEBE exceeding the 1.5K threshold is extremely likely (> 95%) and it occurs much later than that projected by the CMIP5 MME: occurring 20 years after the GCMs: 2070. Similarly, the 2K overshoot, as projected by the FEBE will be avoided with a probability of < 40% but will most likely not be avoided according to the 475 CMIP5 MME with an 89% probability of exceeding 2K.
For the final high emission, business as usual, RCP 8.5 and SSP 585 scenarios; both the FEBE and CMIP5/6 MME project that exceeding the 1.5K threshold is virtually inevitable by 2100. Although in the FEBE projection it is extremely likely that this threshold is exceeded nearly 15 years after the CMIP5/6 MME projections of 2040. The same is found for the 2K threshold, with both the FEBE and CMIP5/6 MME exceeding the threshold about 15 years after the 1.5K threshold.  Figure 13. The probability for the global mean surface temperature of exceeding a 1.5K threshold (top), and a 2K (bottom) are given as a function of years for the FEBE, using FAer SSP (blue), and for the CMIP6 MME (black). The two SSP scenarios are considered for each case: SSP 126 (solid), and SSP 585 (circles).

Conclusions
Multidecadal projections have large uncertainties with ECS confidence limits of 1.5-4.5K essentially unchanged for decades.
For policy makers, the most deleterious consequence is that projections emanating from quite diverse future scenarios have significant overlap. For example, up until 2050, the RCP 2.6 and 8.5 scenarios can both claim to respect the 2K threshold -albeit with rather different probabilities. Large overlaps imply a disconnection between policies (scenarios) and outcomes 485 (temperatures). Now that governments have committed themselves to keeping industrial epoch temperature increases to below 2K (and aim at 1.5K), we face an uncertainty crisis.
One way of reducing this uncertainty is by developing complementary types of models. In this paper we directly constructed a model in the macroweather regime (roughly one month and up) based on the physically principles of energy conservation discovered ( (Lovejoy, 2020b), (Lovejoy, 2020c)) that the FEBE could be derived as a consequence of classical (Budyko, 1969) and (Sellers, 1969), Energy Balance Models (EBMs) that have been regularly used to determine the Earth's latitudinal temperature variations, its stability to perturbations and to study past and future climate states. The key was to introduce a vertical coordinate that allows for the application of the correct conductive-radiative surface boundary condition, needed for correctly determining the energy storage. A surprising consequence is that even the classical (integer ordered) continuum 495 mechanics heat equation used by Budyko and Sellers implies that the surface temperature obeys a fractional ordered energy balance equation. The FEBE's fractional storage terms imply that the system has a long memory so that when calibrated by observational data, its responses to past forcings are constrained to respect the historical climate.
The FEBE is a parsimonious model with only two shape parameters: an exponent H, and relaxation time scale τ ; the classical EBE is the H = 1 special case. In order to make FEBE projections, we use a Bayesian parameter estimation approach. (Hébert 500 et al., 2020) used a Climate Response Function that at long times (> τ ) was close to the corresponding FEBE Green's function, but that was a different power law at shorter time scales. While the (Hébert et al., 2020) CRF was justified on the basis of scale invariance, the FEBE has a stronger physical basis since it respects both scaling as well as energy conservation, and the long memory is explicitly situated in energy storage mechanisms. From the practical (projection) point of view, the main advantage is that the FEBE directly handles the short time scales (down to a month or less). This allows the FEBE to directly 505 take into account the internal variability: a stochastic white noise forcing and the FEBE response. The ability to model the forced response to both external and internal forcing improves FEBE parameter estimates and contributes to lowering the corresponding projection uncertainties.
Bayesian inference allows for a robust probabilistic parameter characterization. The basic external forcings were those prescribed for the historical part of the CMIP5/6 GCMs and these were constrained by five monthly, global resolution empirical 510 temperature series (since 1880). The internal forcing was assumed to be a Gaussian white noise and, since to a good approximation, the FEBE white noise response is a fractional Gaussian noise (fGn), the latter was taken as the Bayesian inference error model.
In order to estimate the parameters, the forcing series required two adjustments. The most important was the aerosol recalibration parameter α which linearly scales the aerosol forcing to take into account the increasing evidence that the CMIP5 515 and CMIP6 aerosol cooling was too strong ( (Padilla et al., 2011), (Hébert et al., 2020), (Zelinka et al., 2020), (Tokarska et al., 2020), (Flynn and Mauritsen, 2020)). The former aerosol series (F Aer RCP ) was based both on uncertain data but also on uncertain modelling assumptions, especially about the direct and indirect effects of aerosols. Whereas the latter (F Aer SSP ) is based on global sulphate production and derived from an alternative model than that in CMIP5.
From our analysis we find a more constrained aerosol forcing. For the F Aer RCP we found a median recalibration factor α of 520 0.6 with 90% CI [0.2,1.0]. Like (Stevens, 2015) this supports a revision of the global modern (2005) (Hébert et al., 2020) and the phenomenological HEBE ( (Lovejoy, 2020b), (Lovejoy, 2020c)). The relaxation time scale τ , characterizing the approach to equilibrium was found to have median value 530 of 4.7 years and 90% CI of [2.4,7.0] years when using F Aer RCP , and nearly identical results using F Aer SSP . The estimated relaxation time is comparable to other box model fast relaxation times ( (Geoffroy et al., 2013), (Schwartz, 2008), (Held et al., 2010), (Rypdal and Rypdal, 2014)) as the one box (EBE) model is a special case of the FEBE where H = 1.
In comparison to IPCC AR5 (and to the CMIP6 MME), we find lower likely ranges for ECS and TCR when using the FEBE with F Aer RCP (or F Aer SSP ). For projections, perhaps the most important parameter is λ, the equilibrium climate sensitivity 535 (ECS) which determines the temperature following an increase in forcing. We see a lower median for the ECS in comparison to the IPCC AR5 (and CMIP6 MME) estimates for their corresponding forcings, the 90% CI range is reduced from [1.5, and the median estimate decreases from 1.8K (2.0K) to 1.5K (1.4K). Several recent observation-based studies ( (Otto et al., 2013), (Skeie et al., 2014), (Johansson et al., 2015), (Lewis and Curry, 2015), (Lewis and Curry, 2018)) have also reported lower ECS upper bounds.
The forcings and parameters combined with the RCP and SSP scenarios allow us to make projections through to 2100, we did this for RCP 2.6 (SSP 126), 4.5, and 8.5 (SSP 585). Overall, the observational based FEBE projections had uncertainties 545 that are smaller by more than a factor of two in comparison to the CMIP5/6 MME uncertainties. However, the two modelling approaches have quite different sources of uncertainty. Whereas the CMIP5/6 uncertainty is purely due differences in the climates of the GCMs ("structural uncertainties"), the FEBE uncertainty depends largely on the uncertainty of the historical forcings and temperatures, in particular, those associated with aerosols. In fact, a byproduct of the model and Bayesian framework is that we are able to more tightly constrain aerosol forcing, supporting recent literature findings of weaker historical 550 aerosol cooling. As a consequence, the FEBE projections are consistently a little cooler than those of the CMIP5/6 MME but still within the uncertainty bounds of the latter, effectively complementing the GCMs.
During the Paris Conference in 2015, nations of the world strengthened the United Nations Framework Convention on Climate Change by agreeing to holding "the increase in the global average temperature to well below 2 • C above pre-industrial levels and pursuing efforts to limit the temperature increase to 1.5 • C". According to our projections, crossing either of these 555 thresholds is delayed with respect to the CMIP5/6 MME projections but will eventually be passed if strong mitigation is not implemented. To avert a 1.5K warming, drastic cuts would have to be made to global greenhouse emissions, similar to that in RCP 2.6 (and SSP 126), for which we found <10% (<50%) probability of exceeding 1.5K in comparison to the CMIP5 (CMIP6) MME which projects a 67% (>80%). Both the FEBE and CMIP5/6 projections have temperatures surpassing 1.5K in scenarios with weak or no mitigation: RCP 4.5, RCP 8.5 and SSP 585, albeit the FEBE projects this occurring nearly two 560 decades later than the GCMs. The 2K threshold is projected to be avoided by both the FEBE and CMIP5/6 MME if we follow low emission scenarios of RCP 2.6 and SSP 126. The opposite is true for any other emission scenarios; the exceeding of the 2K threshold in this coming century would surely occur. Thus our model reinforces that only strong mitigation scenarios, such as RCP 2.6 and SSP 126, will avoid excursions over these 1.5K and 2K thresholds, although it remains to be seen whether negative emission technologies are feasible and whether the appropriate policies are implemented.

565
In this paper we have presented the FEBE model and projections to 2100 which is an observational-based model physically based on energy and scale symmetries, complementary to GCMs. Future work will explore the full (regional, 2D) FEBE model (Lovejoy, 2020c), which hopefully will constrain and improve future projections. Once the full suite of CMIP6 projections are released, the FEBE could be also used to help understand the generational differences between CMIP models. Extensions to precipitation may also be possible at global and regional scales since the FEBE model is consistent with space-time scaling 570 processes in historical precipitation data (de Lima and Lovejoy, 2015).