Attribution in the presence of a long-memory climate response

Multiple, linear regression is employed to attribute variability in the global surface temperature to various forcing components and prominent internal climatic modes. The purpose of the study is to asses how sensitive attribution is to long-range memory (LRM) in the model for the temperature response. The model response to a given forcing component is its fingerprint and is different for a zero response time (ZRT) model and one with an LRM response. The fingerprints are used as predictors in the regression scheme to express the response as a linear combination of footprints. For the instrumental period 1880–2010 CE (Common Era) the LRM response model explains 89 % of the total variance and is also favoured by information-theoretic model selection criteria. The anthropogenic footprint is relatively insensitive to LRM scaling in the response and explains almost all global warming after 1970 CE. The solar footprint is weakly enhanced by the LRM response, while the volcanic footprint is reduced by a factor of 2. The natural climate variability on multidecadal timescales has no systematic trend and is dominated by the footprint of the Atlantic Multidecadal Oscillation. The 2000–2010 CE hiatus is explained as a natural variation. A corresponding analysis for the last millennium is performed, using a Northern Hemisphere temperature reconstruction. The Little Ice Age (LIA) is explained as mainly due to volcanic cooling or as a long-memory response to a strong radiative disequilibrium during the Medieval Warm Anomaly, and it is not attributed to the low solar activity during the Maunder Minimum.


Introduction
There will always be variability in the Earth's climate, even in the absence of external forcing like variation in solar irradiance, volcanic eruptions, or human-induced changes.The nature of internal climate variability is analogous to the change in weather, just extrapolated to longer spatial and temporal scales.This "song of nature" is comprised of a cacophony of frequencies corresponding to the natural modes of the climate system and forms a background spectrum with a pink-noise character.This means that the power spectral density (PSD) of global temperature to a crude approximation has the form S(f ) ∼ 1/f for frequencies f corresponding to periods from months to millennia.The shape of this spectrum implies that internal variability on low frequencies (long timescales) is strong, and this constitutes a problem when we want to detect climate signals and trends with external causes.Another complication is that there are also inter-nal modes that stand out from this noise, and the separation of these modes from the noise background is not unique and depends on how the noise is modelled.
Signal detection means establishing the statistical significance of a trend, an oscillation, or a spatiotemporal pattern.This is successfully done if we can establish that it is very unlikely that the pattern, or fingerprint, has arisen by chance from the internal background noise.Once fingerprints have been successfully detected, the next issue is to assess their relative weights, or footprints, in the total climate signal.This process is what we call attribution.A particular footprint can in some cases be perceived as the result of a particular cause, such as a well-identified radiative forcing.In that case the footprint can be thought of as the global temperature response to this particular forcing.However, attribution does not have to be causal, which is the case if the footprint is the global temperature manifestation of an internal climate Published by Copernicus Publications on behalf of the European Geosciences Union.
mode.For such a mode, a particular climatic variable or index can serve as a particularly sensitive gauge for this specific mode, and its contribution to the variance of the global temperature signal is the mode's footprint.
A standard method in attribution studies is that of multiple linear regression.The idea is to separate the climate signal into a number of components assumed to represent the climate response to individual forcings in addition to a few prominent internal modes.Each of these components has a certain characteristic fingerprint.In order to determine these fingerprints we need models of some sort.Full-scale AOGCMs (atmosphere-ocean general circulation models) can be used, but often also simpler, conceptual models are useful.The rationale for attribution studies is that even the most advanced climate models may estimate wrongly the magnitude of individual responses, even though they have got the fingerprints right.Hence, we may write the total climate signal T (t) as a linear combination of the fingerprints.The validity of the linear approximation for global climate variables has been documented in AOGCM studies by Meehl et al. (2004).Consider, for instance, the global temperature T (t) and the fingerprints of various forcings and internal modes.Then we may, for instance, select the following model for the explained global surface temperature (this is also called the response variable or the predictand): where S(t), V (t), H (t) are the fingerprints of solar, volcanic, and human-induced (anthropogenic) forcing, and A(t) and E(t) are the fingerprints of the Atlantic Multidecadal Oscillation (AMO) and the El Niño-Southern Oscillation (ENSO), respectively.In regression theory the fingerprints are also called predictors.The fitting parameters (or regressors) f sun , f volc , . . .represent the weight of each fingerprint in the total response and can be estimated by minimizing the leastsquare error with respect to the observed data.These weights take into account that we may not have modelled the magnitude of the individual forcings correctly, or that we have overlooked, or modelled incorrectly, climate feedbacks that operate differently for each forcing.A third possible cause of changed weights is incorrect modelling of the temporal response to the forcing.This will give rise to distorted fingerprints.A measure of how successfully the method attributes variability to the various forcing components is to compute how much of the observed variance is explained by the model.One common problem with this approach is that if there are many causal factors to consider, and hence many parameters to fit, there is a risk of overfitting.This means that a good fit can be obtained even when the result is unphysical.Another problem is that the fingerprints of forcing in general are distorted and delayed by inertia in the climate response caused by slow heat exchange between the ocean sur-face layer and the deep ocean, sea ice, and ice sheets.This inertia may, for instance, lead to a small response to the relatively fast solar cycle forcing, while the response to slow trends in solar irradiance may be stronger but considerably delayed.
Delay effects are generally not accounted for in the regression model Eq. ( 1) if the model defining the fingerprint does not involve a dynamic response to forcing.Some authors include delays by introducing a fixed time shift which is different for each fingerprint (Lean andRind, 2008, 2009;Foster and Rahmstorf, 2011).In these papers delays are introduced for the sole purpose of improving the fit, and they increase the number free parameters in the regression model.Under any circumstance, the delay introduced for volcanic, solar, and ENSO fingerprints is a few months and hence is not detectable in the present analysis, which deals with annual data.However, a decadal delay of the anthropogenic fingerprint found by Lean andRind (2008, 2009) was explicitly represented as a result of ocean heat uptake by Canty et al. (2013).In the present paper this delay due to heat exchange with the deep ocean is represented by the longmemory response.The response function to all forcing components is assumed to have the same shape and involves distortion, not just shifts, of the forcing signals.A conceptual stochastic-dynamic model of such a long-memory dynamic response is described in Rypdal and Rypdal (2014), where it is shown that for the global temperature, this model provides results that are essentially indistinguishable from those obtained from the Coupled Model Intercomparison Project Phase 5 (CMIP5) ensemble of general circulation models for the industrial period with historical forcing.
In its most simple form the stochastic-dynamic model is a zero-dimensional energy-balance model (EBM) of the form where T (t) is a perturbation of the surface temperature from an equilibrium state, F det (t) is the total deterministic forcing, σ w(t) is a white-noise stochastic forcing, and −(1/τ )T (t) the radiation imbalance at the top of the atmosphere.The solution if T (0) = 0 is where the response function G(t) = c exp(−t/τ ) represents the impulse response to a delta-function forcing, and hence τ is the characteristic damping time (time constant).It depends on the effective heat capacity C eff of the combined land and ocean surface layer and the climate sensitivity S as τ = C eff S. The first term T det (t) on the right-hand side is the temperature response to the known (deterministic) forcing.
The second term T stoch (t) is the Ornstein-Uhlenbeck (OU) stochastic process, which in discrete time reduces to the firstorder autoregressive (AR(1)) process.This process is stationary and has an autocorrelation function (ACF) of the form C(t) ∼ exp(−t/τ ).The PSD of this process has the shape of a Lorentzian distribution; it is flat (S(f ) ∼ f 0 ) for f τ −1 and decays as S(f ) ∼ f −2 for f τ −1 .If Eq. ( 3) were a good model for the global surface temperature, the residual T obs (t) − T det (t) should correspond to T stoch and hence be successfully modelled as an OU process.In Rypdal and Rypdal (2014), however, it was shown that this residual does not have a Lorentzian PSD but rather exhibits the power-law form S(f ) ∼ f −β , with β ≈ 0.75.This is a persistent process that exhibits long-range memory and is called a fractional Gaussian noise (fGn).These features are also found in control runs in the CMIP5 models (Østvand et al., 2014a) and in CMIP5 simulations with discontinuous jumps in atmospheric CO 2 concentration, one observes relaxation to equilibrium where a fast response with a time constant of 1-2 years is followed by a slow decay that lasts for centuries (Geoffroy et al., 2013).Rypdal and Rypdal (2014) demonstrated that all this can be modelled by replacing the exponential response function by a power law G(t) = ct β/2−1 in Eq. ( 3).It can be shown that this corresponds to replacing the time derivative in Eq. ( 2) with a fractional derivative; hence we name it the fractional EBM.
Equation (3) suggests that the standard, as well as the fractional, EBM can be viewed as a linear filter that transforms the forcing signal into temperature signal.In the Fourier domain, the equation takes the form T (f ) = G(f ) [ F (f ) + σ w(t)], and for the PSD we get In the absence of deterministic forcing ( F (f ) = 0), we have will be a Lorentzian distribution, and the resulting stochastic process is the OU process.If G(t) is the power law , and the process is an fGn.In the absence of stochastic forcing, the filter represented by |G(f )| 2 will suppress only fluctuations on timescales smaller than τ if G(t) is exponential, while the power-law filter will systematically suppress small scales and enhance large scales.Examples were shown by Rypdal and Rypdal (2014) where a time series for the total forcing throughout the last 130 years is run through an exponential filter with τ = 4.3 years and a power-law filter with β = 0.75 (long-memory response).One observes that only the latter is able to reproduce a realistic response to the negative forcing due to volcanic eruptions (the negative spikes in the forcing signal).It also provides a better (although not perfect) fit to the large-scale trends in the observed temperature signal.
The long-memory response has important implications for the prediction of future global temperature on a century timescale.In Fig. 1 it is shown that in a moderately pessimistic forcing scenario for the next 100 years, the frac-tional, long-memory model predicts a temperature almost 1 • C than the zero response time model.The latter projection does not change much with an exponential response as long as τ is less than 1 decade.
The purpose of this paper is to assess the sensitivity of the attribution to the assumption of long, versus short, memory in the computation of the fingerprints associated with volcanic, solar, and anthropogenic forcing.Section 2 describes briefly the multiple-regression method and the regression diagnostics used, although these are very standard.Section 3.1 presents results based on instrumental surface temperature data and forcing reconstruction for the period 1880-2010 CE (Common Era), and Sect.3.2 presents the same analysis using a millennium-long multi-proxy reconstruction of Northern Hemisphere temperature and radiative forcing.Section 4 presents the paper's conclusions and discusses their implications.

Data and methods
The forcing data in this paper are given as the annual and global mean of the radiative forcing measured in W m −2 .The data from the instrumental period 1880-2010 CE are those used by Hansen et al. (2005Hansen et al. ( , 2011) ) and those for the reconstruction period 1000-1979 CE are from Crowley (2000).The instrumental temperature data are given as annual and global mean surface temperature anomalies relative to 1880 CE (the HadCRUT3 data set; Brohan et al., 2006) and the reconstructed temperature data as Northern Hemisphere annual means relative to 1000 CE (Moberg et al., 2005).The forcing data are split up in solar, volcanic, and anthropogenic components.There are more recent instrumental data sets, but for the analysis in the present paper they will only provide unimportant corrections.The reason for employing these older data sets is that it allows the use of the parameters estimated, and comparison to the results obtained, in a recent paper (Rypdal and Rypdal, 2014).
In this paper we shall compare the effects of two different response filters: the zero response time filter F ZRT and the long-range memory filter F LRM .Mathematically they represent two extremes, although we shall see that the LRM (long-range memory) filter is quite an accurate representation of the actual response.If the total forcing is written as F (t) = F sun (t) + F volc (t) + F anthr (t), and F is the filter operator, then we construct the response function (predictand), and determine the regression coefficients c 0 , c 1 by a simple least-square fit.The response function Q(t) is the fitted, filtered response to the total forcing F (t) and can be considered as the best model we can make for the temperature signal with the filter F, without allowing for different weights of the individual fingerprints.These fingerprints are defined as Blue curve in all panels is the instrumental GMST (global mean surface temperature) recorded in the period 1880-2010 CE.Panel (a): the ZRT (zero response time) regressed signal Q(t) defined in Eq. ( 5), with F being the ZRT (identity) filter and F (t) the total forcing.Panel (b): the ZRT 3P regressed signal according to Eq. ( 1) without AMO and ENSO as predictors.Panel (c): the LRM (long-range memory) regressed signal Q(t) defined in Eq. ( 5), with F being the LRM filter.Panel (d): the LRM 3P regressed signal according to Eq. ( 1) without AMO and ENSO as predictors. follows: The responses Q(t) are plotted for the two filter models in Figs.2a, c and 6a, c to provide an indication of what the filtered response will be like if we do not allow for individual feedbacks to the different forcing components.The next step is to allow for such individual weights and determine them by the construction of the linear predictand shown in Eq. (1).Our first choice is to leave out the AMO and ENSO predictors, leaving us with solar, volcanic, and anthropogenic forcing as predictors.With zero response time filter and these three predictors we have the ZRT 3P (zero response time with three predictors) regression model.The corresponding case with an LRM filter is the LRM 3P model.Including AMO and ENSO (five predictors) gives us the ZRT 5P and LRM 5P models.The weighted responses is our best estimate of climate footprints imposed by the forcing or internal modes characterized by the corresponding fingerprints.
The estimation of the regression coefficients and some diagnostics are done by the command "LinearModelFit" in "Mathematica".For each predictand Q(t) we provide the R 2 diagnostic (coefficient of determination), which measures the fraction of the total variance in the observed record that is explained by the predictand.As we move from one to three, and then to five, predictors (and the same number of fitting parameters), we increase model complexity and will increase the explained variance.In model selection assessments we have model selection criteria based on information theory where the likelihood function is used as a measure of the goodness of the fit, which is subject to a penalty for model complexity.The most commonly used of these are the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) (for an introduction to the concepts, see Burnham and Anderson, 2004).Each of these criteria produces a real number that can be positive or negative, and the model giving the smaller number is preferable in this information-theoretic sense.

Attribution from instrumental data
In this section we use the same instrumental global temperature data and the forcing data as employed by Hansen et al. (2011).The analysis is based on the annual mean time series.

Zero response time model
Figure 2a (red curve) shows the predicted signal obtained by fitting the unfiltered forcing (more precisely, by fitting Q(t) given by Eq. ( 5) subjected to the zero response time filter F ZRT ) to the instrumental GMST (global mean surface temperature; blue curve).The fit is quite poor (R 2 ≈ 0.53), and the response to the volcanic eruptions are obviously much stronger than observed.If we include an exponentially decaying response exp(−t/τ ), we will need a time constant τ larger than 1 decade in order to obtain realistic short-time responses to these eruptions (see Rypdal, 2012), provided we do not reduce the weight of the volcanic forcing.Another way of obtaining a better fit is to employ a multiple regression by using Eqs.( 1) and ( 6).The result is shown in Fig. 2b.The fit is much better (R 2 ≈ 0.80), but there is a rather strong decadal oscillation attributable to the solar cycle.The redistribution of weights is apparent from Fig. 3a and b.The three fingerprints given in Fig. 3a are just a rescaling of the three forcing components by the same factor c 1 given by Eq. ( 6) with F = F ZRT , and the red curve in Fig. 2a is just the unweighted superposition of these fingerprints plus the additive constant a 0 .The multiple three-component regression ZRT 3P is the superposition of the weighted fingerprints (i.e., the footprints) shown in Fig. 3b.The regression amplifies the solar fingerprint S(t) by a factor f sun ≈ 2.10 and the anthropogenic fingerprint by f anthr ≈ 1.58, while the volcanic fingerprint is strongly attenuated, with f volc ≈ 0.22.This strong attenuation is provoked by the unrealistically large short-time responses enforced by the zero response time model, and the suppression of the volcanic cooling is what has to be compensated for by amplified solar and anthropogenic warming.Thus, for the ZRT response model the strongly altered weights are most probably caused by an incorrect (too spiky) representation of the volcanic fingerprint.

The long-range memory response model
The ZRT response model is given by the delta function G(t) = cδ(t) and is obviously unrealistic.Next, we explore the effect of an LRM response function of the form G(t) = (t/µ) β/2−1 .In Rypdal and Rypdal (2014) a maximumlikelihood approach was applied to estimate µ = 0.84 × 10 −2 years and β = 0.75 from the same instrumental temperature data and forcing data as used in the present paper.Figure 2c shows the response variable given by Eq. ( 5), with FT (t) representing the LRM filter, and c 0 = 0.15 × 10 −2 K and c 1 = 0.92 determined by fitting Eq. ( 5) to the instrumental observation data.The fact that c 0 is close to zero and c 1 is close to unity shows that the least-square fit for these data gives results compatible with the more general maximum-likelihood approach employed in Rypdal and Rypdal (2014).Compared to the ZRT-filtered response the explained variance R 2 is increased from 0.53 to 0.81.This is partly due to a better representation of the largescale variability and a smaller immediate response to the volcanic eruptions due to the memory effects.The explained variance is only slightly increased by introducing variable weights on the solar, volcanic, and anthropogenic fingerprints (R 2 ≈ 0.83), and the improvement is mostly caused by a suppression of the volcanic response.Compared to the fingerprints shown in Fig. 3c, the volcanic footprint shown in Fig. 3d is reduced by a factor f volc ≈ 0.53, while the solar footprint is only slightly amplified by f sun ≈ 1.18 and the human footprint slightly attenuated by f anthr ≈ 0.90.The AIC and BIC are somewhat reduced, so this model is preferable compared to the unweighted LRM model, but the difference is not very large.With respect to explained variance and the information-theoretic selection criteria, the ZRT 3P, LRM, and LRM 3P models are quite similar.However, visual inspection of the shape of the responses and footprints suggests that the ZRT 3P model results in the suppression of volcanic footprint and in an amplification of the solar footprint that is unrealistically large.Similarly, the reduction in the volcanic footprint in the LRM 3P model by a factor of approximately 0.5 seems to give a much better fit to the short-time temperature around the large volcanic eruptions and suggests that the volcanic forcing signal in the forcing data may have been exaggerated.ity is interpreted as residual noise.However, some variability manifest in the global temperature is not adequately described as long-memory or short-memory noise.The ENSO signal is easily detected in the global temperature records, and even though El Niño or La Niña events are unpredictable, the PSD of ENSO indices peaks in the frequency range corresponding to periods of between 2 and 7 years.It is therefore common to include ENSO in attribution analyses (Lean andRind, 2008, 2009;Foster and Rahmstorf, 2011).Another feature that appears impossible to explain with only forcing predictors is the low temperatures in the first decades of the nineteenth century and the high temperatures in the decades after World War II.These anomalies may be compatible with an oscillation with a period of 60-70 years, as discussed extensively by Canty et al. (2013).The statistical significance of this oscillation with respect to a long-memory null hypothesis for the noise background was discussed by Østvand et al. (2014b) but has also been studied extensively by a number of authors with short-memory null models (Ghil and Vautard, 1991;Schlesinger and Ramankutty, 1994;Plaut et al., 1995;Polonski, 2008).The mode has the same period and phase as the most prominent period in the AMO index, and thus it seems reasonable to introduce the AMO index as a predictor variable in addition to the Niño3 index if one wants to increase the explained variance.One could object that the inclusion of temperature observations as predictor variables looks like circular reasoning.However, as mentioned in the introduction, regression is not really about attributing observed variance to external causes but rather about attributing global temperature variability to a set of signatures (fingerprints).These may signify responses to forcing (causation) but also the footprint in the global temperature of observed climate signals like the North Atlantic sea surface temperature or pressure differences across the tropical Pacific.However, one should also bear in mind that climate forcing is a problematic concept, since the separation of forced from internal dynamics depends on what one defines as the "system" that is subject to external forcing.The reasoning above was based on thinking about ENSO and AMO as internal modes and not as forcing.However, one can also define the Earth surface-ocean mixed layer as the system, and this system can be forced by modes involving energy exchange between the surface-mixed layer and atmospheric systems (ENSO) and between the surface-mixed layer and the deep ocean.AMO is an example of the latter and can be considered as the fingerprint of the forcing exerted on the ocean mixed layer from the Atlantic Meridional Overturning Circulation (AMOC) (DelSole et al., 2013;Medhaug and Furevik, 2004;Canty et al., 2013).

Inclusion of internal modes
The AMO index and an index for ENSO (Niño 3.4) are shown in Fig. 4. Using the LRM fingerprints for S(t), V (t), and H (t), and the AMO index for A(t) in Eq. ( 1 1) including the three forcings and AMO (but not ENSO) as predictors.Panel (b): the footprints f sun S(t) (yellow), f volc V (t) (magenta), f anthr H (t) (green), and f AMO A(t) of LRM 4PA regressed signal.Panel (c): the LRM 5P signal, i.e., the same as in (a), but with the ENSO signal added in the regression.Panel (d): the LRM 5P footprints, i.e., the same as in (b), but with the ENSO signal added in the regression.The ENSO footprint is the orange dotted curve.
prints relative to the LRM3 model shown in Fig. 3d, apart from a notable reduction in the volcanic footprint.A similar effect of using AMO as a predictor was found by Canty et al. (2013).Hence, the effect of including AMO as a predictor is mainly to increase the explained variance, but we also note a hiatus in the first decade of the twenty-first century.In Fig. 5c and d we show the effect of adding the Niño3 index as a predictor, in addition to AMO.The explained variance is raised to R 2 ≈ 0.89, and the AIC and BIC are further reduced, suggesting that both AMO and ENSO are relevant explanatory variables and that including both contributes to a better statistical model.The hiatus post 2000 CE is even more pronounced when ENSO is included, due to the strong 1998 El Niño.
The total natural footprint (the sum of solar, volcanic, AMO, and ENSO footprints) is dominated by the multidecadal oscillation with a weak growing trend caused by the growing trend in solar activity in the period 1880-1960 CE.From Fig. 5d we observe that this trend in the solar footprint is very close to the trend in the anthropogenic footprint up to t ≈ 90 (1970 CE), but after this time the solar footprint has no significant trend, while the trend in the anthropogenic footprint is approximately 0.13 K per decade.The anthropogenic footprint turns out to be very robust and quite insensitive to the inclusion of natural modes in the regression analysis.

Attribution from multiproxy data
A similar analysis is made using the Northern Hemisphere multiproxy temperature reconstruction of Moberg et al. (2005) and the forcing reconstruction of Crowley (2000) for the period 1000-1979 CE.The data are given with annual resolution, but since the temperature data are effectively smoothed on timescales shorter than 5 years, it seems unrea-   5), with F being the SRM filter and F (t) the total forcing.Panel (b): the SRM 3P regressed signal according to Eq. ( 1) without AMO and ENSO as predictors.Panel (c): the LRM regressed signal Q(t) defined in Eq. ( 5), with F being the LRM filter.Panel (d): the LRM 3P regressed signal according to Eq. ( 1) without AMO and ENSO as predictors.
sonable to use a zero response time model.Instead a shortrange memory (SRM) response model with an exponential response function G(t) = c exp(−t/τ ) is employed.The parameters c = 0.37 K yr −1 and τ = 4.3 years were estimated by Rypdal and Rypdal (2014) using the instrumental data over the period 1880-2010 CE.As for the instrumental data we will also use the LRM model with parameters estimated from the instrumental data.By employing the models with these parameters we can examine how well the SRM model works versus the LRM model for a longer data set.This is of interest because the SRM model employed with the instrumental data explains almost as large a fraction of the variance as the LRM model; hence, from these data the LRM is not strongly preferable to the SRM model based on the R 2 and AIC and BIC criteria only.
The SRM response is shown in Fig. 6a and the corresponding fingerprints in Fig. 7a.The response function does not show a good fit (R 2 ≈ 0.17) and AIC and BIC are large.The introduction of weighted fingerprints increases the explained variance (R 2 ≈ 0.26) and lowers AIC and BIC as shown in Fig. 6b.As shown in Fig. 3a and b, this improvement comes about by a considerable reduction in the volcanic footprint (f volc ≈ 0.45) and from a very strong amplification of the solar footprint (f sun ≈ 2.32).The volcanic footprint is reduced to lower the variance due to the sharp spikes in the SRM response to the volcanic forcing, and the solar footprint is amplified to reduce the unexplained variance from the cooling between the Medieval Warm Anomaly (MWA) and the Little Ice Age (LIA).The anthropogenic footprint is also amplified (f anthr ≈ 1.48).However, the LRM model in-creases the explained variance to R 2 ≈ 0.39 and drastically reduces AIC and BIC, even without introducing weighted fingerprints.This is shown in Fig. 6c and demonstrates that the LRM model is strongly preferable to the SRM model when we consider timescales up to 1 millennium.The consistency of the LRM model is supported by the observation that the introduction of weighted fingerprints introduces weights moderately different from unity.The main change is an enhancement of the solar footprint (f sun ≈ 1.44) at the expense of the volcanic one (f volc ≈ 0.78).The anthropogenic footprint is virtually unchanged (f anthr ≈ 0.99).This tendency to an enhanced solar, reduced volcanic, and only slightly affected anthropogenic footprints is consistent with what was observed when the LRM model was applied to the instrumental data.
Figure 7d suggests that the temperature difference between the maximum of the MWA (1000 CE) and the minimum of the LIA (1700 CE) can be mainly attributed to volcanic cooling, while the warming from the LIA until 1970 CE is attributed to solar and anthropogenic influence.The latter is also consistent with what we observe from Fig. 7c.
For all response models the explained variance is considerably lower for the reconstruction data than for the instrumental data.This is mainly due to the strong anthropogenic trend in the instrumental period.This trend dominates the variance and is very well predicted; hence, it increases the predicted variance.6), with F being the SRM filter to the individual forcing components (F sun (t) (yellow), F volc (t) (magenta), F anthr (t) (green)) to produce the fingerprints S(t), V (t), and H (t) for the SRM filter.Panel (b): the footprints f sun S(t) (yellow), f volc V (t) (magenta), and f anthr H (t) (green) of SRM 3P regressed signal.Panel (c): the same as in (a) but with the LRM filter.Panel (d): the same as in (b) but with the LRM filter.

Effect of initial state and prehistory
By defining the fingerprints as integrals over the time interval (0, t), we implicitly assume that there is no influence of past forcing from the interval (−∞, 0), i.e., we effectively assume zero forcing in pre-history.For the exponential (SRM) response function this has no consequence because this response function corresponds to the simple EBM which is just a first-order ordinary differential equation whose solution only depends on the initial temperature (see discussion in Rypdal and Rypdal, 2014).For the power-law response, prehistory matters in principle, since the corresponding differential equation contains a fractional derivative.However, even for the simple SRM response, we cannot assume with certainty that the initial forcing is zero, since this would be the same as assuming that the climate system is in equilibrium at time t = 0.As this may have some surprising implications, some detail may be appropriate.Consider as an illustration the simple zero-dimensional EBM where T is surface temperature in Kelvin, C is an effective heat capacity per area of the Earth's surface, σ S is the Stefan-Boltzmann constant, is the effective emissivity of the atmosphere, and I (t) is the incoming radiative flux density at the top of the atmosphere.Let T 0 = T (t = 0), I 0 = I (t = 0), T = T 0 + T , and I = I 0 + F .Note that F here is the perturbation of the radiative flux with respect to the initial flux I 0 and not with respect to the flux that would be in equilibrium with the initial temperature T 0 .The linearized equation for the temperature change relative to the temperature T 0 is The quantity I (eq) 0 ≡ σ S T 4 0 represents the incoming flux required to balance the outgoing long-wave radiation (OLR) from the top of the atmosphere when the surface temperature is T 0 .This is not necessarily equal to the actual incoming flux at time t = 0, so the difference F 0 = I 0 −σ S T 4 0 represents the initial forcing (or the initial imbalance of radiative flux density).By definition F (0) = I (0) − I 0 = 0 and represents the sum of various forcing components that we have used to establish the forcing fingerprints, as they are all defined to be zero at t = 0.The solution to Eq. ( 9) takes the form To understand the implications let us look at the case where F (t) is a stationary stochastic process, implying that the expectation value E[F (t)] is independent of t.Eq. ( 10) then describes realizations of the response to this forcing under the condition that the initial radiative imbalance is F 0 .The expectation (ensemble average) of this response is then www.earth-syst-dynam.net/6/719/2015/Earth Syst.Dynam., 6, 719-730, 2015 If we assume that the response is also stationary, this equation implies that E[F ] = −F 0 .What this means is that if the initial imbalance is a fluctuation around a stationary climate state but everything is measured as perturbations relative to the initial state, then the mean of the forcing F in the future will balance the initial forcing F 0 to render T finite.On the other hand, if F and T are not stationary processes, this is no longer true.Consider for instance that the imbalance F 0 is the result of a step in radiative influx just prior to t = 0.If the radiative flux is held constant after this step, we would have that F (t) = 0 for t > 0 and the evolution would be given by the first integral in Eq. ( 10).
For an exponential response function, this contribution converges to a constant for t τ , but for a power-law response, it takes the form T ∼ t β/2 .The divergence as t → ∞ is of course unphysical.It reminds us that the power-law response is an idealized representation of the response of a system with a large range of response times and that it must be cut-off on some timescale (Rypdal and Rypdal, 2014).However, it illustrates that if parts of the climate system respond very slowly, there may be a strong influence of an initial energy imbalance throughout the entire temperature record under consideration.
The effect of past forcing is a less serious problem.It was shown in Rypdal and Rypdal (2014) to be negligible over the instrumental period, using information about forcing and temperature over the past millennium.We have not carried out a similar computation for the millennium period, since reliable global-scale reconstructions for the previous millennia are not available.However, for past forcing to have a long-term effect, the climate system must have been driven strongly away from radiative equilibrium over an extended period.This is the case in the Anthropocene but is not believed to have occurred throughout earlier millennia in the Holocene.
We do not have direct physical information about the radiative flux imbalance in the year 1000 CE, but the high temperatures during the MWA might suggest that the OLR was higher than the incoming flux at the start of the subsequent cooling.What we can do by means of attribution techniques is to include T F 0 = F 0 t 0 G(t − t ) dt as an extra fingerprint and estimate F 0 along with the other regression coefficients.The results are shown in Fig. 8.The total response in Fig. 8a explains more variance than the model that does not include T F 0 , and the AIC and BIC prefer this model.In particular, the large discrepancy between explained and observed variability during the first century of the record (during the MWA, 1000-1100 CE) in the other models has disappeared in this long-memory four-predictor (LRM 4P) model.Figure 8b shows a strong reduction in the volcanic footprint because the long-term trend imposed by F 0 provides the cooling previously attributed to volcanic activity.It is quite apparent from Fig. 8a that the estimated response exhibits a weak short-term response to volcanic eruptions, but the estimated f volc ≈ 0.28 is only 30 % lower than what was esti-mated from the LRM 5P model applied to the instrumental data (Fig. 5d).The volcanic footprint may have been somewhat underestimated in Fig. 8a, simply because the shortterm response does not contribute very much to the total variance.However, recent work on detection and attribution which compares multiproxy reconstructions with paleoclimatic simulations to general circulation models show that the models seem more sensitive on short timescales to volcanic eruptions than observed in the reconstructions (Schurer et al., 2013).Many explanations can be offered for this observation, and one could be that the volcanic forcing used in the models, or its efficacy, has been systematically overestimated.Hence, it is difficult to rule out that the tendency shown in Fig. 8 could be more than an analysis artifact.

Conclusions
Standard linear, multiple regression has been applied to instrumental and multiproxy reconstructed global and northern hemispheric temperatures, using fingerprints derived from reconstructed forcing and internal mode indices as predictor variables.The fingerprints have been derived from simple short-memory and long-memory response models.The regression coefficient for the volcanic fingerprint will be strongly suppressed by zero response time and short-memory response models, but the explained variance is still around 80 %.The modelling of Lean and Rind (2008) is similar to our ZRT 3P but with the inclusion of finite time delays and ENSO as an additional predictor.Their results shown in their Fig. 2 are quite similar to the ZRT 3P results shown in Figs.2b and 3b in this paper, with an explained variance of 76 %.Hence, the inclusion of ENSO and finite time delay without the memory smoothing of the response does not seem to improve the explained variance, while the increased model complexity necessarily will increase the AIC and BIC scores and make the model less preferable.The model of Lean and Rind (2008) suffers from the same problem as the ZRT 3P model in that it overestimates the 11-year solar cycle response by not taking into account the attenuating effect of long-memory response to oscillatory forcing on a decadal timescale.In Rypdal (2012) it was shown that the large solar cycle response of 0.2 K peak to peak detected by Camp and Tung (2007) in global surface temperatures in the period 1959-2004 CE are largely attributable to three volcanic eruptions taking place incidentally in the descending phase of solar cycles.By correcting for the responses to these eruptions, there will be a considerably weaker response to the solar cycle in the global temperature series.
Multiple regression based on fingerprints derived from long-memory response models, and in particular with AMO and ENSO included as predictors, yields a response variable that explains 89 % of the total variance of the instrumental data set for 1880-2010 CE.Relative to the forcing data set employed for this period, the solar footprint is mod- ified by a factor f sun ≈ 1.23, the volcanic footprint by a factor f volc ≈ 0.41, and the anthropogenic footprint by a factor f anthr ≈ 0.77.In the instrumental period the natural variability is dominated by an internal oscillation with a period of 60-70 years, and this oscillation dominates over the forced trend up to 1970 CE.The forced trend before 1970 CE is shared in equal proportion between solar and anthropogenic footprints.After 1970 CE the trend in the anthropogenic footprint is approximately 0.13 K per decade, but the trend in the total response has been amplified by the upward phase of the AMO footprint and the strong El Niño in 1998.The combination of these footprints and that of the Mount Pinatubo eruption in 1991 CE yields a total response function showing a hiatus in the years 2002-2010 CE.Lean andRind (2008, 2009) attribute much of this hiatus to the descending solar cycle, while in the present analysis shown in Fig. 5 the maximal phase of the AMO in 2010 CE and the 1998 El Niño give more important contributions.A recent update of the sea surface temperature (SST) has cast doubt on the reality of the hiatus in global temperature (Karl et al., 2015).This is consistent with the present results, since these corrections to the SST also pertain to the AMO and ENSO fingerprints.Correction of these fingerprints will probably eliminate the hiatus in the LRM 5P response shown in Fig. 5c.The solar-cycle fingerprint, however, is unaffected by these corrections, so in the model of Lean andRind (2008, 2009), the hiatus will persist in their modelled response despite these corrections of the observed temperature.By including the AMO as an additional predictor in the LRM model response (but leaving out ENSO), the volcanic response is reduced by a factor of 2. This somewhat surprising result is due to the structure of the volcanic forcing over the instrumental period, with a 5-decade-long period of low volcanic forcing prior to the Mount Agung eruption in 1963 CE, and a series of major disruptions including El Chichón (1982 CE) andPinatubo (1991 CE).Without the AMO the post World War II cooling will be attributed exclusively to volcanic cooling, while the inclusion of AMO will lead to about half of this cooling being attributed to the low phase of the AMO.This finding is in close agreement with those of Canty et al. (2013), and in our paper it depends on the LRM character of the response, which we believe is strongly related to the overturning circulation.Without this delayed response (as illustrated by the results for the ZRT 3P model in Fig. 3), the volcanic response would be anomalously low both with and without AMO.
For the millennium reconstruction the short-memory response with a time constant of 4.3 years is unable to reproduce the reconstructed long timescale variability.The longmemory response offers two viable models for the large-scale variability: one where most of the cooling from the MWA to the LIA is attributed to volcanic activity; the other attributes more of this cooling to a negative radiative imbalance at the end of the MWA, represented as a negative initial forcing at 1000 CE, and giving rise to a downward temperature trend throughout the last millennium.Both explanations require there to be a significant long-memory impact up to millennium timescales.
The regression examples shown in this paper demonstrate that the results of attribution studies based on multiple linear regression depend strongly on the memory properties of the models employed to define the fingerprints.Models including long-term memory in the response tend to explain more of the observed variance and have better scores on information-theoretic model selection tests.Results also vary with the number and nature of the fingerprints used as predictors.Nevertheless, there are some tendencies that seem to be robust throughout.The weight of the anthropogenic footprint is not systematically changed by treating the individual forcing components as independent predictors, and almost all of the global warming since 1970 CE can be attributed to it.The solar footprint is enhanced by a factor of approximately 2 with a short-memory response but is not changed a lot with a long-memory response.The volcanic footprint is strongly suppressed with a short-memory response and is also somewhat weaker with a long-memory response.Even though the solar footprint is enhanced in all models, none of them attributes the Little Ice Age primarily to solar variability.

Figure 1 .
Figure 1.The light blue curves in both panels are observed global temperature in the period 1880-2010 CE.The orange curve on the left is the historic anthropogenic forcing extended with an exponential growth in atmospheric CO 2 concentration, ending at around 700 ppm 100 years from now.The red curve on the right shows the projected temperature from a standard EBM with τ = 0, and the blue curve shows this from a fractional EBM with β = 0.75.The difference between the two projections in the year 2100 CE is almost 1 • C.

Figure 3 .
Figure 3. Fingerprints and footprints for the instrumental temperature 1880-2010 CE.Panel (a): the application of Eq. (6) with F being the ZRT filter for the individual forcing components (F sun (t) (yellow), F volc (t) (magenta), F anthr (t) (green)) to produce the fingerprints S(t), V (t), and H (t) for the ZRT filter.Panel (b): the footprints f sun S(t) (yellow), f volc V (t) (magenta), and f anthr H (t) (green) of ZRT 3P regressed signal.Panel (c): the same as in (a) but with the LRM filter.Panel (d): the same as in (b) but with the LRM filter.

Figure 5 .
Figure 5. Panel (a): the LRM 4PA signal, i.e., the regressed signal according to Eq. (1) including the three forcings and AMO (but not ENSO) as predictors.Panel (b): the footprints f sun S(t) (yellow), f volc V (t) (magenta), f anthr H (t) (green), and f AMO A(t) of LRM 4PA regressed signal.Panel (c): the LRM 5P signal, i.e., the same as in (a), but with the ENSO signal added in the regression.Panel (d): the LRM 5P footprints, i.e., the same as in (b), but with the ENSO signal added in the regression.The ENSO footprint is the orange dotted curve.

Figure 6 .
Figure 6.Blue curve in all panels is the Moberg reconstructed temperature for the Northern Hemisphere plotted for the interval 1000-1979 CE.Panel (a): the SRM regressed signal Q(t) defined in Eq. (5), with F being the SRM filter and F (t) the total forcing.Panel (b): the SRM 3P regressed signal according to Eq. (1) without AMO and ENSO as predictors.Panel (c): the LRM regressed signal Q(t) defined in Eq. (5), with F being the LRM filter.Panel (d): the LRM 3P regressed signal according to Eq. (1) without AMO and ENSO as predictors.

Figure 7 .
Figure 7. Fingerprints and footprints for the Northern Hemisphere temperature for 1000-1979 CE.Panel (a): the application of Eq. (6), with F being the SRM filter to the individual forcing components (F sun (t) (yellow), F volc (t) (magenta), F anthr (t) (green)) to produce the fingerprints S(t), V (t), and H (t) for the SRM filter.Panel (b): the footprints f sun S(t) (yellow), f volc V (t) (magenta), and f anthr H (t) (green) of SRM 3P regressed signal.Panel (c): the same as in (a) but with the LRM filter.Panel (d): the same as in (b) but with the LRM filter.

Figure 8 .
Figure 8. Panel (a): the same as in Fig. 6d but with inclusion of the LRM response F 0 to the initial forcing F 0 as a predictor.Panel (b): the corresponding footprints.The smooth brown curve is T F 0 .