On the linearity of the temperature response in Holocene : the spatial and temporal dependence

Previous studies show that the evolution of global mean temperature forced by the total forcing is almost the same as the sum of those forced by individual forcing in the last 21,000 years as simulated in three independent climate models, CCSM3, FAMOUS and LOVECLIM. But how does the linearity of the climate response depend on the spatial and temporal scales quantitatively and in what regions the linear response tends to dominate remain unknown. Here, based on the 15 TraCE-21ka climate simulation outputs, the spatial and temporal dependence of the linear response of temperature evolution in the Holocene is studied using correlation coefficient and a linear error index. The results show that the linear response of global mean temperature is strong on the orbital, millennial and centennial scales throughout the Holocene, but is poor on the decadal scale. The linear response differs significantly between the Northern Hemisphere (NH) and Southern Hemisphere (SH). In the NH, the linear response is strong on millennial scale, while in the SH the linear response is strong on orbital 20 scale, such that the linear response of the global mean orbital variability is dominated by the SH response, but that of the global millennial variability seems dominated by the NH response. Furthermore, at the regional scales, the linear responses differ substantially between the orbital, millennial, centennial and decadal timescales. On the orbital scale, the linear response is dominant for most regions, even at the small area of about a mid-size country like Germany. On the millennial scale, the linear response is still strong in the NH over many regions, albeit weaker than on the orbital scale. However, the 25 millennial variability shows relatively poor linear response over most regions in the SH. On the centennial and decadal timescales, the linear response is no longer significant in almost all the regions. The regions of strong linear response on the millennial scale are mostly consistent with those on the orbital scale, notably western Eurasian, North Africa, subtropical North Pacific, tropical Atlantic and Indian Ocean, because of a large signal-to-noise ratios over these regions. This finding can improve our understanding of the regional climate response to various climate forcings in the Holocene. 30 Clim. Past Discuss., https://doi.org/10.5194/cp-2018-177 Manuscript under review for journal Clim. Past Discussion started: 7 January 2019 c © Author(s) 2019. CC BY 4.0 License.


Introduction
Long term climate evolution in the Holocene (~ last 11,000 years) is believed to be mainly driven by several external forcing, notably, the orbit forcing, the greenhouse gases (GHGs), the meltwater forcing and ice sheets forcing (COHMAP members, 1988;Liu et al., 2014.Here, we treat the climate system as the coupled ocean-atmosphere system, such that Earth orbital parameters, GHGs, meltwater discharge and continental ice sheet are considered as external forcing).Conceptually, in the simplest linear climate model,   = − ×  + , where the total forcing  consists of multiple forcing  = ∑    =1 , the response is linear such that the climate response to the total forcing  equals the linear sum of all the responses to individual forcing   .In the real world, however, climate evolution in the Holocene may exhibit nonlinear features (deMenocal et al., 2000).Furthermore, there is strong internal climate variability of a wide range of timescales, from daily variability of synoptic weather storms (Hasselmann, 1976), to interannual variability of El Nino (Cobb et al., 2013), interdecadal climate variability (Zhang et al., 1997;Kushnir, 1994;Delworth and Mann, 2000), all the way to millennial climate variability (Bond et al., 1997).Therefore, climate response to external forcing during the Holocene may not exhibit linear response everywhere over the globe, especially at the regional scale where internal climate variability of interannual to interdecadal timescales is strong in the coupled ocean-atmosphere system.However, the effect of internal variability should be suppressed by an average over space and/or time.Furthermore, one can expect that the effect of internal variability will be suppressed more as the average spatial and temporal scales increase.Therefore, it is conceivable that the linearity of the climate response depends on both the spatial and temporal scales, with the linearity more dominant for larger spatial scale and longer timescale.In case of the largest scale of the global mean surface temperature, the evolution of global mean temperature forced by the total forcing is almost the same as the sum of those forced by individual forcing in the last 21,000 years as simulated in three independent climate models, CCSM3, FAMOUS and LOVECLIM (Fig. 2 of Liu et al., 2014).This naturally leads to two questions.First, how does the linearity of the climate response depend on the spatial and temporal scales quantitatively?Second, in what regions, the linear response tends to dominate?The answer to these questions will great improve our understanding of regional climate response during the Holocene.
In this paper, we address the spatial and temporal dependence of the surface temperature response on climate forcings in the Holocene using a set of climate simulations in CCSM3 (Liu et al., 2014).We will analyse the linearity of the climate response on the timescales of orbital, millennial, centennial and decadal, and at the spatial scales of global, hemispheric and regional.The data and methodology are given in section 2. The dependence of linearity on spatial and temporal scales are analysed in section 3. A summary is given in section 4.

Data
The data is from TraCE-21ka (Liu et al., 2009 and2014), which consists of a set of synchronously coupled atmosphereocean general circulation model simulations for the last 21,000 years.The simulations are completed using the CCSM3  (Shakun et al., 2012;Marsicek et al., 2018), suggesting a reasonable climate sensitivity in CCSM3.In addition to the all forcing run (ALL), there are four individual forcing runs forced by the orbital forcing (ORB), the continental ice sheets (ICE), the GHGs (GHG) and meltwater forcing (MWF), respectively (Liu et al., 2014, Table 1).Here, we examine the surface temperature response in the Holocene (last 11, 000 years), although the simulations go back to the last glacial maximum which is about 21,000 years.Because the Holocene period poses a more stringent test of the linearity of the climate response, as it removes the deglacial warming response that is dominated by the response to increased CO2.Unless otherwise specified, each time series is averaged in 100-yr bins of the length of 110.The 100-year average filters out variability under centennial timescale.
It should be noted that, the continental ice sheet topography and sea-land distribution remains at the condition of 19 ka BP except in the all forcing run and ice sheet run, in which the continental ice sheet topography and sea-land distribution are changed with time.

Methods
We use both the temporal correlation coefficient  and a linear error index   to evaluate the goodness of the linear response.
The correlation coefficient is derived from the formula below Here,   is the linear SUM of the temperature time series of the four single forcing experiments,   is the full temperature response time series in the ALL run, both at time , and  is the length of the time series.The overbar represents the mean of the entire time series.The correlation coefficient represents the similarity of the temporal evolution between the linear SUM response and the ALL response.The linear error index   represents the magnitude of the relative error between the linear sum and the full response.It is defined as the root mean square error (RMSE) between the sum temperature response and the full temperature response divided by the standard deviation of the full temperature response in the ALL run: A higher  and a smaller   represent a better linear response.
The dependence of the linear response on spatial and temporal scales will be studied by filtering the time series in different scales.For the spatial scale, we will divide the globe in 9 succeeding cases, denoted by 9 division factors: f=0, 1,2,3,4,6,8,12 and 24 from the largest global scale to the smallest model grid scale.The f=0 case is for global average while the f=1 case is for hemispheric average in the NH and SH.Further division will be done within each hemisphere.Note that each hemisphere has 96(.) × 24(.) grid boxes, with a ratio of 4:1 between longitude span and latitude span.We divide each hemisphere into  ×  sections of equal latitude and longitude spans, with each area containing the same number of (96/) × (24/) grid boxes, maintaining the ratio of 4:1 between longitude span and latitude span.For example, f=2 is for the 2 × 2 division, with each area containing 48 × 12 grid boxes; f=24 is for the 24 × 24 division with each area containing 4×1 grid boxes, about the size of 15°(.) × 3.75°(.), like the size of a mid-size country of Spain or Germany in the mid-latitude.
To study timescale dependence of the linear response, we will decompose a full 100-yr binned temperature time series into three components of, roughly, orbital, millennial and centennial timescales.Following Marsicek et al. (2018), we derive the orbital and millennial variability using the locally weighted regression fits (Loess fits) (Cleveland, 1979).First, we apply a 6500-yr Loess fit onto the temperature series of 100-yr bin to derive the orbital scale variability, which contains the trend and the slow evolution longer than ~6500 years.Second, we apply a 2500-yr Loess fit onto the 100-yr binned temperature to derive a low frequent component, which is then used to subtract the 6500-yr Loess fit to produce the millennial variability.
Finally, the 100-yr binned temperature minus the 2500-yr Loess fit gives the centennial variability.Also, we use running mean on the 10-yr binned temperature time series to derive the decadal scale variability.First, we apply a 100-yr running mean on the time series of the 10-yr bins (of the length of 1100).Second, we use the 10-yr binned temperature minus the 100-yr running mean to derive the decadal scale variability.
Statistical significance of the  of a particular timescale is tested using the Monte Carlo method (Kroese, 2011 and2014;Kastner, 2010;Binder, 1997) with 1,000,000 realizations on the corresponding red noise in the AR1 model (autoregressive model of order 1).As a benchmark, the significance level is tested against the full global mean temperature series in the ALL run.For the total, orbital, millennial, centennial and decadal temperature time series, the 95% confidence levels are 0.72 (with the AR1 coefficient 0.96), 0.76 (0.97), 0.65 (0.95), 0.21 (0.31) and 0.19 (0.06), respectively.
Statistical significance of the   of a particular timescale is tested using the bootstrap method (Efron, 1979 and1993) with 1,000,000 realizations on the corresponding time series.As a benchmark, the significance level is tested against the global mean temperature series in the ALL run.For the total, orbital, millennial, centennial and decadal temperature time series, the 95% confidence levels are 1.23, 1.23, 1.21, 1.24 and 1.36, respectively.

Linear responses at different temporal scales
The global mean temperature provides a useful example to study the dependence of the linear response on timescales.We first examine the linear response of the global mean temperature based on its components of orbital, millennial, centennial and decadal variability (Fig. 1).Fig. 1a is the total variability of global surface temperature derived from the ALL run and the SUM.The global temperature evolution shows a strong linear response on the orbital and millennial scales throughout the Holocene (Fig. 1b and 1c).The orbital scale evolution is characterized by a warming trend of about 1℃ from 11ka to ~4.5ka before decreasing slightly afterwards.This feature is captured in the linear sum albeit with a slightly smaller magnitude and an additional local minimum around 3ka (Fig. 1b).The total variability is very similar to the orbital variability (the correlation coefficient is 0.99, Fig. 1a vs Fig. 1b).The millennial variability shows 5 major peaks around ~9.8ka, 7.8ka, 4.7ka, 3.7ka and 1.8ka.All these peaks seem to be captured in the linear SUM response albeit with a slightly larger amplitude (Fig. 1c).For the orbital and millennial variability, the correlation coefficients between the sum and the full responses are =0.83and 0.71, respectively, both significant at the 95% confidence level and explaining over 50% of the variance; the linear errors are   =0.63 and 0.92, respectively, also significant at the 95% confidence level.The global centennial variability also appears to be a modest linear response (Fig. 1d), with a correlation of =0.44 and linear error of   =1.21.But the linear response of the global decadal variability is poor (Fig. 1e), with a correlation of =-0.02 and linear error of   =1.99, which is not statistically significant at the 95% confidence level.The result of the analysis on global mean temperature is qualitatively consistent with previous hypothesis that the linear response tends to degenerate at shorter temporal scale.It shows that the global linear response is largely valid at orbital and millennial timescales, weakened substantially at centennial scales and is invalid at decadal scales.

Linear responses at different spatial scales
In order to assess the linear response at smaller spatial scales, we first analyse the linear response of the NH and SH (f=1).It is interesting that the linear response significantly differs between the NH and SH (Fig. 2).Fig. 2a and 2f are the total variability of hemispheric surface temperature derived from the ALL run and the SUM for the NH and SH, respectively.
Their components on the 4 timescales (i.e.orbital, millennial, centennial and decadal) are shown in the Fig. 2b-2e and Fig. 2g-2j, respectively.In the NH, the linear response is strong at the millennial scale (=0.82,  =1.01, Fig. 2c), but not so strong on orbital (=0.55,  =0.84, Fig. 2b) and centennial (=0.32,  =2.29, Fig. 2d) scales, since only   is significant on the orbital scale, and only  is significant on the centennial scale.At the decadal scale, the linear response fails completely (= -0.04,   =1.99, Fig. 2e), as expected.In comparison, in the SH, the linear response is dominant at the orbital scale (=0.92,  =0.43, Fig. 2g), but poor on all other timescales, including millennial (=-0.12,  =2.32, Fig. 2h), centennial (=0.14,   =3.19, Fig. 2i) and decadal (=0.03,  =2.07, Fig. 2j) scales.The linear response of the global mean temperature discussed in Fig. 1 seems to be dominated by the SH response at the orbital scale, but by the NH response on the millennial scale.This suggests that the linear responses of the NH and SH are complex and should be treated with caution.This further highlight the need to study the linear response at regional scales.
The linear response at different regions on the orbital, millennial, centennial and decadal timescales are summarized in Fig. 3 and Fig. 4 in the correlation coefficient and linear error index, respectively.Fig. 3a  features.First, as expected, the correlation coefficient tends to decrease towards smaller area (larger f).Quantitatively, however, the correlation coefficient does not decrease much, such that even at the smallest area (f=24), the correlation in most regions are still above 0.8, statistically significant at 95% confidence level.This suggests that the orbital evolution in most regions exhibit a strong linear response, even at the smallest scale of about a mid-size country like Germany (f=24, 15°(.) × 3.75°(.)).Second, the linear response in the NH is slightly better than SH (for f≥3).The reason will be discussed later.Third, subareas in both hemispheres show comparable linear response across all the spatial scales, with the median correlation all above ~0.8,except that of the NH mean temperature (f=1.The linear response of NH is not better than that of regional variability at smaller spatial scales f≥2).This is a special feature which should be treated with caution.The correlation coefficient therefore shows that, for orbital scale evolution, temperature response is dominated by the linear response over most of the globe, even at regional scales.These features are also consistent with the linear error analysis in Fig. 4a.
The millennial variability also shows a weaker linear response for smaller scales, in both the correlation (Fig. 3b) and error (Fig. 4b).Quantitatively, the millennial variability still exhibits strong linear response in the NH over many regions, albeit weaker than at the orbital scale.The correlation coefficients remain above 0.6 across most regions even at the smallest division area (f=24), contributing to ~40% of the variance.In contrast to the orbital variability, where regions in both hemispheres show comparable linear response, the millennial variability shows poor linear response in most regions in the SH, with the median no longer significant at 95% confidence level.As for orbital scale, the linear responses of all the spatial scales in the NH are better than those in the SH.
In contrast to the orbital and millennial variability, there is almost no strong linear response for the centennial variability.
The median linear response in centennial timescale in either hemisphere across spatial scales (f>3, Fig. 3c and Fig. 4c) is no longer significant, with few correlation coefficients larger than 0.3 and contributing less than 10% of the variance.Finally, the decadal variability exhibits absolutely no linear response over any spatial scales in both of the NH and SH (Fig. 3d and Fig. 4d).
The strong linear response at the orbital and millennial scales suggest that these two groups of variability are generated predominantly by the external forcing.In contrast, the poor linear response of centennial and decadal variability suggest that these two groups of variability are caused mainly by the internal coupled ocean-atmosphere processes.

Pattern of the Linear Responses
We now further study the pattern of the linear response.Fig. 5 shows the spatial patterns of the correlation coefficients at orbital (Fig. 5a1-a3) and millennial (Fig. 5b1-b3) for three representative spatial scales, f=3, 6 and 24 (the other factors are not shown because they are similar to the above mentioned three representative spatial scales).For orbital variability (Fig. 5a1-a3), the linear response is strong in most regions in the NH in all three spatial scales, with the correlation coefficient above 0.8.In the SH, the linear response is also strong over the continents, but is poor over the ocean.This leads to the significantly reduced linear response in the SH as discussed in Fig. 3a-4a.This suggests that orbital variability is likely forced For millennial variability, in the NH, the linear response shows a similar feature to that of the orbital variability, but that is poor in almost all the SH.Fig. 5b1-b3 show that the linear response is strong in most regions in the NH at the three spatial scales, with the correlation coefficient above 0.6.At the regional scale, e.g.f=6 and 24 (Fig. 5b2 and 5b3), the linear response is strong over the northwestern Eurasian continent, Northern Africa, northern North America, northern Pacific Ocean, southern North Atlantic Ocean and western Arctic Ocean, as seen in the Fig. 5b2 and 5b3 (only those significant correlation coefficients are shown).But the linear response is poor over the southern North America, the eastern and southern Eurasian continent.While in the SH, the linear relationship is poor over almost the entire SH as seen in the correlation map (Fig. 5b1-b3).Interestingly, over the NH, the regions of strong linear response for millennial variability are mostly consistent with the regions of strong linear response for orbital variability, notably western Eurasian, North Africa, subtropical North Pacific, tropical Atlantic and Indian Ocean.These preferred region of linear response suggests some common mechanism of the climate response in these regions.
In order to understand the reason for the preferred regions of linear response, we examine the signal-to-noise ratio.Since our study above shows that the linear response is largely valid for orbital and millennial variability, but not for centennial and decadal variability, we define the variance of the orbital and millennial variability crudely as the linear signals, while define the variance of the sum of the centennial and decadal variability, which is dominated by internal variability, as the linear noise.The ratio of the signal to noise will be used for shedding light on the regional preference of linear response.Figure 6 shows the signal-to-noise ratio for orbital and millennial variability for the three representative spatial scales (f=3, 6 and 24).
For orbital variability (Fig. 6a1-a3), the signal-to-noise ratio is large (above 10 (0.6) , the log base 10 is taken on the signal to noise ratios on the orbital scale) in most regions in the NH.In the SH, the signal-to-noise ratio is also large over the continents, but is small over the ocean.At different regional scales, e.g.f=6 and 24 (Fig. 6a2 and 6a3), the signal-to-noise ratio is large over the western half of the Eurasian continent, Northern Africa, Central and South America, most NH oceans and SH tropical oceans, and the Antarctic Continent.But the signal-to-noise ratio is small over the North American continent and the eastern Eurasian continent, and the entire Southern Ocean.These features of spatial pattern of signal-to-noise ratio on the orbital scale (Fig. 6a1-a3) are similar to those in the correlation map (Fig. 5a1-a3).The spatial correlation between the map of the signal-over-noise ratio in Fig. 6 and the corresponding correlation coefficient in Fig. 5 is 0.53, 0.58 and 0.49 for f=3, 6 and 24 respectively (Fig. 7), all significant at 95% confidence level.For millennial variability, the signal-to-noise ratio also shows a similar feature to that of orbital variability although overall somewhat smaller (note the different color scales!).Fig. 6b1-b3 show the signal-to-noise ratio is large in most regions in the NH in all three spatial scales, with the signal-to-noise ratio above 0.15.At the regional scale, e.g.f=6 and 24 (Fig. 6b2 and 6b3), the signal-to-noise ratio is large over the northwestern Eurasian continent, Northern Africa, central North America, northern North Pacific Ocean, Northern Atlantic Ocean and the Arctic Ocean.But the signal-to-noise ratio is poor over the southern and northern North America continent, the South America and the eastern and southern Eurasian continent.While the signal-to-noise ratio is poor over almost the entire SH.Over the NH, interestingly, the regions of large signal-to-noise ratio for millennial variability, are mostly consistent with the regions of strong signal-to-noise ratio for orbital variability, notably northwestern Eurasian, North Africa, subtropical North Pacific, Northern Atlantic and tropical Northern Indian Ocean.These features of spatial pattern of signal-to-noise ratio on the millennial scale (Fig. 6b1-b3) are similar to those in the correlation map (Fig. 5b1-b3).The correlation coefficient between the maps of signal-to-noise ratio and correlation is 0.73, 0.43 and 0.27 for f=3, 6 and 24, respectively (Fig. 7), again all significant at 95% confidence level.

Conclusions
It is found that the linear response of the temperature tends to degenerate at shorter temporal scale and smaller regional scales.The global mean temperature evolution shows a linear response on the orbital, millennial, and centennial scales throughout the Holocene, but shows no linearity for decadal variability (Fig. 1).Further analysis on regional scale suggests that the linear response is strong on the orbital and millennial scale for most continental regions over the NH and SH.
However, the linear response is poor over much of the ocean, especially over the ocean in the SH.There are specific regions where linear response tends to be dominant, notably the western Eurasian continent, North Africa, the central and South America, the Antarctica continent, and the North Pacific.The strong linear response is interpreted as the region of large signal-to-noise ratio.That is, in these regions, either the orbital and millennial response signal is large, or the influence of the centennial and decadal variability noise is small, or both.This suggests that the orbital and millennial variability in these regions are relatively easy to understand.This finding can help paleoclimatological scientists to understanding the effects of climate forcings.
Much further work needs to be done.The result here represents the first such assessment and would need to be substantiated with other models.Furthermore, the physical mechanism for the preferred region of linear response also needs to be better understood.
shows the correlation coefficients of the orbital variability in each region for the 9 division factors.The cases of global mean (f=0) and hemispheric mean (f=1) have been discussed in details before.The correlation coefficients for succeeding division factors (f=2, 3, …, 24) show several Clim.Past Discuss., https://doi.org/10.5194/cp-2018-177Manuscript under review for journal Clim.Past Discussion started: 7 January 2019 c Author(s) 2019.CC BY 4.0 License.
Clim.Past Discuss., https://doi.org/10.5194/cp-2018-177Manuscript under review for journal Clim.Past Discussion started: 7 January 2019 c Author(s) 2019.CC BY 4.0 License.predominantly by the external forcing over continents.Over the ocean, especially the SH oceans, strong internal variability may manifest itself all the way to affect the orbital scale variability.More specifically, at the regional scales, e.g.f=6 and 24 (Fig.5a2 and 5a3), the linear response is strong over the western half of the Eurasian continent, Northern Africa, Central and South America, most NH oceans and SH tropical oceans, and the Antarctica Continent, as seen in the Fig.5a2 and 5a3(only those significant correlation coefficients are shown).But the linear response is poor over the North America continent and the eastern Eurasian continent, and the entire Southern Ocean.Similar features can also been seen in the map of linear error (not shown).

Figure 1 :
Figure 1: The global annual mean surface temperature time series derived from the ALL run (black) and the SUM (the sum of four single forcing run, red).In the (a), the thin line is the 100-yr binned time series, the thick line is 2500-yr loess fitted time series.The orbital scale variability (b) is represented by the 6500-yr loess fitted series.The millennial variability (c) is represented by the 2500-yr loess fitted data subtracting the 6500-yr loess fitted data.The centennial variability (d) is represented by the 100-yr binned data subtracting the 2500-yr loess fitted data.The decadal variability (e) is represented by the 10-yr binned origin time series subtracting its 100-yr running mean.The correlation coefficient is given at the upper left corner and the linear error is given at the lower left corner of each panel.The x axis is year (ka, 0 is AD 1950, negative is before AD 1950), and the y axis is temperature anomaly (℃, relative to 11ka-0ka).10

Figure 2 :
Figure 2: The surface temperature time series derived from the ALL run (black) and the SUM (red).The (a-e) is similar to the Fig.1a-e but for NH, and the (f-j) is similar to the Fig1.a-e but for SH.The x axis is year (ka, 0 is AD 1950, negative is before AD 1950), and the y axis is temperature anomaly (℃, relative to 11ka-0ka).

Figure 3 :
Figure 3: Correlation coefficient of mean surface temperature between the ALL and SUM outputs in different regions (red dots for global, blue dots for NH and cyan dots for SH) and timescales (a, orbital; b, millennial; c, centennial; d, decadal).The red line is connecting the medians of all blue and cyan dots at all division factors.The blue (cyan) line is connecting the medians of all blue (cyan) dots at all division factors.The black thick dashed line is 95% confidence level.The black thin solid line is zero.The x axis 5

Figure 4 :
Figure 4: Same as Figure 3, but for linear error (  , y axis).The black thick dashed line is 95% confidence level of linear error.There are only about 10 points with the Le larger than 6 so they are neglected in a-c.

Figure 5 :
Figure 5: The correlation coefficient of mean surface temperature between the ALL and SUM outputs of the two timescales (a1-a3, orbital; b1-b3, millennial) for three representative spatial scales, f=3, 6 and 24 (the other factors are not shown because they are similar to the abovementioned three representative spatial scales).Only those significant at 95% confidence level are shaded in colors.

Figure 6 :
Figure 6: The signal-to-noise ratios on the orbital (a1-a3) and millennial (b1-b3) time scales derived from the ALL run.Here the signal is defined by the orbital (millennial) variability variance, and the noise is defined by the sum of the centennial and decadal variabilities variance.The numbers of 1, 2 and 3 are for the three representative spatial scales, f=3, 6 and 24, respectively.In order to show the signal clearly, the log base 10 is taken on the signal to noise ratios on the orbital scale (a1-a3).

Figure 7 :
Figure 7: Scatter diagram of the correlation coefficient and the signal to noise ratio on the orbital (a) and millennial (b) time scale.The log base 10 is taken on the signal to noise ratio on the orbital scale, as indicated in Figure 6.Only those for the three representative spatial factors (f=3, 6 and 24) are shown in the panels.The blue stars and linear fitting line are for factor 3, the red 5