Causes of Greenland temperature variability over the past 4000 yr : implications for northern hemispheric temperature changes

Precise understanding of Greenland temperature variability is important in two ways. First, Greenland ice sheet melting associated with rising temperature is a major global sea level forcing, potentially affecting large populations in coming centuries. Second, Greenland temperatures are highly affected by North Atlantic Oscillation/Arctic Oscillation (NAO/AO) and Atlantic multidecadal oscillation (AMO). In our earlier study, we found that Greenland temperature deviated negatively (positively) from northern hemispheric (NH) temperature trend during stronger (weaker) solar activity owing to changes in atmospheric/oceanic changes (e.g. NAO/AO) over the past 800 yr (Kobashi et al., 2013). Therefore, a precise Greenland temperature record can provide important constraints on the past atmospheric/oceanic circulation in the region and beyond. Here, we investigated Greenland temperature variability over the past 4000 yr reconstructed from argon and nitrogen isotopes from trapped air in a GISP2 ice core, using a one-dimensional energy balance model with orbital, solar, volcanic, greenhouse gas, and aerosol forcings. The modelled northern Northern Hemisphere (NH) temperature exhibits a cooling trend over the past 4000 yr as observed for the reconstructed Greenland temperature through decreasing annual average insolation. With consideration of the negative influence of solar variability, the modelled and observed Greenland temperatures agree with correlation coefficients of r = 0.34–0.36 (p = 0.1–0.04) in 21 yr running means (RMs) and r = 0.38– 0.45 (p = 0.1–0.05) on a centennial timescale (101 yr RMs). Thus, the model can explain 14 to 20 % of variance of the observed Greenland temperature in multidecadal to centennial timescales with a 90–96 % confidence interval, suggesting that a weak but persistent negative solar influence on Greenland temperature continued over the past 4000 yr. Then, we estimated the distribution of multidecadal NH and northern high-latitude temperatures over the past 4000 yr constrained by the climate model and Greenland temperatures. Estimated northern NH temperature and NH average temperature from the model and the Greenland temperature agree with published multi-proxy temperature records with r = 0.35–0.60 in a 92–99 % confidence interval over the past 2000 yr. We found that greenhouse gases played two important roles over the past 4000 yr for the rapid warming during the 20th century and slightly cooler temperature during the early period of the past 4000 yr. Lastly, our analysis indicated that the current average temperature (1990–2010) or higher temperatures occurred at a frequency of 1.3 times per 1000 yr for northern high latitudes and 0.36 times per 4000 yr for NH temperatures, respectively, indicating that the current multidecadal NH temperature (1990–2010) is more likely unprecedented than not ( p = 0.36) for the past 4000 yr.


Introduction
Greenland temperature change is one of the most important factors for determining the impacts of future climate change on our society, as it is expected to undergo large changes owing to polar amplification (Meehl et al., 2007), and the T. Kobashi et al.: Causes of Greenland temperature variability over the past 4000 yr substantial influence that the melting of the Greenland ice sheet has on global coastal areas due to sea level rise (Jacob et al., 2012).However, Greenland temperature variability in multidecadal to centennial scales is poorly understood owing to the lack of appropriate long temperature records.In this study, we utilized a recently developed Greenland temperature record for the past 4000 yr, which was estimated from argon and nitrogen isotopes from trapped air in the GISP2 ice core (Kobashi et al., 2011), and investigated causes of Greenland temperature variability over the past 4000 yr with a 1-D energy balance model with reconstructed climate forcings.The temperature estimates are seasonally unbiased decadal averages (10-20 yr in central Greenland; Spahni et al., 2003) as the heat, gas diffusion, and bubbleclose-off process in the firn layer (unconsolidated permeable snow layer) smoothes large seasonal and inter-annual temperature variations (Severinghaus et al., 1998).
Greenland is located in the northern North Atlantic, where climatic variation is considerably affected by the North Atlantic Oscillation/Arctic Oscillation (NAO/AO).Atmospheric dipoles of high and low pressure at the surface are located near the Azores and Iceland, respectively, and fluctuation in pressure differences between the two centres of action induces variations in westerlies in the North Atlantic and northerly winds in Greenland (Hurrell, 1995;Wanner et al., 2001;Hurrell et al., 2003).Therefore, Greenland undergoes cooling when northern Europe experiences warming (positive NAO), and vice versa (Hurrell, 1995).Climate modelling and observations suggest that solar variability induces changes in atmospheric circulation similar to NAO/AO through ozone feedback in the stratosphere (Shindell et al., 2001;Kodera and Kuroda, 2002;Gray et al., 2010;Scaife et al., 2013).When solar activity is stronger (weaker), changes in the positive (negative) NAO/AO-like atmospheric circulation are induced (Shindell et al., 2001;Lean and Rind, 2008).Therefore, it can be expected that stronger (weaker) solar activity induces warming (cooling) in NH temperature, and relative cooling (warming) in Greenland through positive (negative) NAO.Consistent with this theory, Greenland temperatures have deviated negatively (positively) from the NH temperature trend when solar activity was stronger (weaker) over the past 800 yr (Kobashi et al., 2013).Climate modelling also indicates that the Atlantic meridional overturning circulation (AMOC) reduces (increases) during weaker (stronger) solar activity (Cubasch et al., 1997;Waple et al., 2002), contributing to negative Greenland temperature responses to solar variability (Kobashi et al., 2013).Currently, the past variations of NAO and/or AMOC (e.g. the Medieval Climate Anomaly to Little Ice Age) are highly debated (Mann et al., 2009;Trouet et al., 2009Trouet et al., , 2012;;Olsen et al., 2012).Therefore, proper understanding of Greenland temperature variability over the past 4000 yr will provide important constraints for these debates as well as future NAO and AMOC behaviour in a warmer climate.
Many palaeotemperature proxies are biologically mediated, generally recording spring to summer temperatures (Kaufman et al., 2004(Kaufman et al., , 2009)), and often the relationships between climatic and biological parameters are not stationary (Huntley, 2011).The story becomes even more complicated for these temperature proxies on a longer timescale as seasonal variation of solar insolation has changed radically over the Holocene (Laskar et al., 2004).To derive hemispheric to global average temperatures on a multidecadal to centennial timescale, temperature proxies with larger uncertainties in temperature and chronology are compiled such that the multidecadal to centennial trends are frequently smoothed out, with the trends increasing as they go further back in time (Wanner et al., 2008).Therefore, we explored another way to constrain multidecadal to millennial hemispheric temperature variability using the model and the Greenland temperature record.
In this study, we first introduce methodology and the Greenland temperature in Sect. 2.Then, in Sect. 3 we investigate spatial patterns of temperature change associated with the Greenland Summit temperature.In Sect.4, we examine climate forcings (orbital, solar, volcanic, greenhouse gas, aerosol forcings), and employ a 1-D energy balance model to explore the impacts of climate forcing on temperatures over the past 4000 yr.In Sect.5, we investigate causes of Greenland temperature changes over the past 4000 yr.In Sect.6, we explore corresponding implications for NH temperature changes based upon the relationship between the NH and Greenland temperatures (Kobashi et al., 2013) and the climate model, followed by conclusions in Sect.7.

Ice core chronology
Greenland ice cores provide one of the best chronologies to study the Holocene climate owing to the existence of annual layers (Alley et al., 1997a;Vinther et al., 2006).For the data generated from the GISP2 ice core, we employed the "visual stratigraphy" chronology (Alley et al., 1997b), which allowed us to compare temperature changes to other climatic information in the ice core with minimum age uncertainties (e.g.GISP2 δ 18 O ice and sulfate for volcanic signals).The uncertainty of the chronology has been estimated to be 1 % (Alley et al., 1997b).The ages of gas and ice differ in ice cores because the age of the gas is determined at the time of bubble closure at the bottom of the snow layer.Differences between gas age and ice age were estimated using a firn-densification/heat-diffusion model (Goujon et al., 2003;Kobashi et al., 2011).Uncertainty in the differences is approximately 10 %, or approximately 20 yr, as the gas age to ice age differences are approximately 200 yr in central Greenland during the past 4000 yr.The uncertainty estimate is likely conservative as timings of the changes in δ 18 O ice and temperature agree well (Kobashi et al., 2010).Therefore, absolute age uncertainties are less than 45 yr at approximately 2000 BCE (= √ 40 2 + 20 2 ; from the theory of error propagation) and less than 20 yr at approximately 1950 CE.Comparison with the GICCO05 chronology (S.O. Vinther and Rasmussen, B. M., personal communication, 2012) indicates that these are appropriate estimates.The use of the GICO05 chronology had little effect on the following analyses, as the differences between GICO05 and visual stratigraphy are relatively small over the past 4000 yr.

Correlation coefficients and significance
The Pearson product-moment correlation coefficient was calculated to evaluate the strength of the correlations.Because of the autocorrelations of the time series, it is important to consider the reduced degree of freedom for the calculation of significance.The number of effective degrees of freedom (df e ) for p statistics was estimated from the effective sample size (N e ), which depends on the effective decorrelation time (T e ) (Ito and Minobe, 2010;Zie ¸ba, 2010).The two relevant time series are represented by X and Y , and the sample size is N. Here, the time interval ( t) is 1 yr.Therefore, N e = N t/T e and df e = N e − 2. For the correlation coefficients, T e is estimated as , where τ is a lag or lead in years, and ρ is the correlation coefficient of the 2 time series indicted by the subscripts (XX and YY are autocorrelations).The range of (−N, N) was defined with τ starting from zero (no lead or lag; ρ = 1) and moving to an increased lag or lead with a 1 yr step until ρ drops to zero (Ito and Minobe, 2010).Then, p values of the correlation coefficients were calculated from a two-sided Student t distribution with the effective samples size (N e ) described above.We considered a 90 % confidence level (p = 0.1) as a significant correlation but reported all p values where p is less than 0.1 to allow the assessment of significance.

Greenland temperature reconstruction over the past 4000 yr
There have been several attempts to reconstruct the variability of Greenland temperature, which have included the use of classic methods of application of oxygen isotopes of ice (δ 18 O ice ) (Dansgaard, 1964;Jouzel et al., 1997;White et al., 1997;Johnsen et al., 2001;Vinther et al., 2010), borehole thermometry (Alley and Koci, 1990;Dahl-Jensen et al., 1998), and combinations of the two techniques (Cuffey et al., 1995;Cuffey and Clow, 1997).However, δ 18 O ice is known to be affected by factors other than local temperature change, such as changes in storm tracks, seasonal accumulation, vapour-source temperature, and the location of evaporative origins (Charles et al., 1994;Masson-Delmotte et al., 2005;LeGrande and Schmidt, 2009; Steen-Larsen et al., 2011).Physically constrained borehole thermometry is more accurate, but it loses information on shorter term temperature variation rather quickly (Dahl-Jensen et al., 1998).
Recently, Vinther et al. (2009) reconstructed Greenland temperature variations from δ 18 O ice and borehole temperature profiles while considering changes in ice sheet elevation.
A Greenland temperature record for the past 4000 yr (Fig. 1) has been developed (Kobashi et al., 2010(Kobashi et al., , 2011) ) using nitrogen and argon isotopes from trapped air in the GISP2 ice core (Kobashi et al., 2008b).Gases in the firn layer (unconsolidated permeable snow layer) fractionate according to the depth of the layer (gravitational fraction) (Craig et al., 1988;Schwander, 1989) and the temperature gradient ( T ; thermal fractionation) in the layer (Severinghaus et al., 1998).Therefore, previous information on firn depths and temperature gradients can be obtained by measuring two isotope ratios from trapped air in ice cores (e.g.nitrogen and argon) as deviations from present values in the atmosphere (Severinghaus et al., 1998;Kobashi et al., 2008b).This method has been applied to quantify the magnitudes of abrupt climate changes by calibrating δ 18 O ice with firndensification/heat-diffusion models to match measured nitrogen and argon isotopes (Severinghaus et al., 1998;Lang et al., 1999;Severinghaus and Brook, 1999;Landais et al., 2004;Huber et al., 2006;Kobashi et al., 2007).
As δ 18 O ice may not reflect local temperature changes, especially for periods such as the relatively stable Holocene, it was necessary to develop a new method to calculate surface temperatures directly from nitrogen and argon isotopes without relying on δ 18 O ice .For this purpose, we developed a method (Kobashi et al., 2008a(Kobashi et al., , 2010(Kobashi et al., , 2011) ) to calculate surface temperature change directly by integrating T (derived from nitrogen and argon isotopic data) with a firndensification/heat-diffusion model (Goujon et al., 2003).We note that the reconstructed temperature records for the past 1000 yr reported in Kobashi et al. (2010Kobashi et al. ( , 2011) ) are nearly identical, attesting that slight differences in the calculation do not affect the outcome.The reconstructed Summit temperatures over the past 160 yr match well with independently reconstructed Summit temperatures from observation and climate models (Box et al., 2009) within uncertainties (Kobashi et al., 2011).The average standard error of the mean over the past 4000 yr is estimated to be 0.8 • C (1σ ) as in Fig. 1 (Kobashi et al., 2011).For the past 1000 yr, the uncertainty is smaller (1σ = 0.5 • C) because sample density is about 6 times higher than the earlier period (Kobashi et al., 2008b(Kobashi et al., , 2011)).
The temperature record shows higher correlation (r = 0.49, p = 0.05) with stacked δ 18 O ice records (GISP2, Dye-3, GISP2, GRIP, NGRIP, and Agassiz; Fig. 2, top panel) than individual δ 18 O ice records from Greenland region likely because stacking reduces glaciological noises as suggested by earlier studies for shorter duration (Fisher et al., 1996;White et al., 1997).The potential sources of the error of the technique using argon-nitrogen isotopes include measurement errors for argon and nitrogen isotopes, uncertainty in temperature calculation using firn-densification/heat-diffusion models, potential changes in the depth of the convective zone in the firn layer that should be minimum especially during the Holocene in Greenland, and argon isotope correction (Kobashi et al., 2010(Kobashi et al., , 2011)).To understand these uncertainties as well as regional variability of Greenland temperatures, it is necessary to conduct temperature reconstructions from other Greenland ice cores.An important consideration for the temperatures derived from the Greenland Summit is an elevation change of the Greenland ice sheet.However, the change occurred likely at a minimum over the past 4000 yr (e.g.< 100 m) (Vinther et al., 2009).

Spatial correlation between Summit temperatures and other areas of Greenland and the world
The Greenland temperature record (Kobashi et al., 2011) was derived from the GISP2 ice core from the Summit region of central Greenland.To make inferences about the past 4000 yr of climate from a temperature record, it is critical to understand how the Summit temperature is connected with other parts of Greenland and the world for the observational period of the past 160 yr.Box et al. (2009) produced Greenland grid temperatures from 1840 to 2007 (now updated to 2011).Using the grid temperatures, a correlation map was created for the period of 1852-2010 (a period for which continuous data are available) for the Summit record from a fusion of observations and a climate model by Box et al. (2009).East Greenland exhibits higher correlations, but most of Greenland temperatures are highly correlated with the Summit temperatures as the 5th percentile of all the grid correlation coefficients is r = 0.73 (Fig. 2), indicating that 95 % of the grid correlation coefficients are higher than r = 0.73.This is in line with the fact that the observation-based Summit Ice core sites -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Greenland Summit (Box et al., 2009) vs global grid temperature (GISS) for 1880-2010.

Fig. 2.
Correlation maps between Greenland Summit temperatures and grid temperatures over the period.Top panel: correlation map in the period between 1852 and 2010 with the observation-modelbased Summit and grid temperatures (annual resolution) by Box et al. (2009).Parentheses indicate correlation coefficients between the Summit temperatures (Kobashi et al., 2011) and δ 18 O ice from six ice cores (Vinther et al., 2009) over the past 4000 yr with p = 0.04 for Dye3, p = 0.02 for GISP2, p > 0.1 for GRIP, p = 0.01 for NGRIP, p = 0.05 for Renland, and p > 0.1 for Agassiz.Middle panel: correlation map between Greenland Summit temperature (Box et al., 2009) and global grid temperatures for the period from 1880 to 2010 (Hansen et al., 2010).Only grids with > 90 % confidence are shown (two-sided tests).Bottom panel: correlation map between Greenland Summit temperature (Kobashi et al., 2011) and global grid temperatures (Mann et al., 2009) for the period from 1200 to 2010.Only grids with > 90 % confidence are shown (twosided tests).
temperature is highly correlated (r = 0.92; p < 0.01) with the average Greenland ice sheet surface air temperature calculated by Box et al. (2009).Therefore, the Summit temperature record should be a useful indicator of the melting of the ice sheet in coastal regions in the past, and consequently changes in the Greenland ice sheet.
To further extend our inferences of the Summit temperature, we examined the correlation map between the Greenland Summit temperatures (Box et al., 2009) and global grid surface air temperatures (GISS) (Hansen et al., 2010) for the period from 1880 to 2010 (Fig. 2).The correlation map indicates a significantly high positive correlation with the northern North Atlantic, western Mediterranean, Middle East, and North Africa, resembling one of the two regimes of the NAO/AO dipole pattern (Hurrell, 1995;Marshall et al., 2001).Some tropical regions of the Atlantic, Indian, and Pacific oceans, where the NAO has been linked with tropical dynamics (Hurrell et al., 2003), also exhibit high correlation.Interestingly, rather high correlations (r = 0.7 to 0.80) are found in geographically distant regions such as the tropical Atlantic, eastern Mediterranean, and Middle East (Red Sea and Persian Gulf), with correlations as high as near Iceland, indicating strong teleconnection patterns for these regions with Greenland over the past 130 yr (Fig. 2).
One exception from the NAO pattern is the location of the mid-latitudes of North America, which is usually identified as part of the other side of the NAO dipole.This may be explained by the fact that the climate of the region is also affected by the Atlantic multidecadal oscillation (AMO) (Enfield et al., 2001;McCabe et al., 2004).The general spatial patterns of the correlation over the past 130 yr were also identified in a correlation map between the grid proxy temperatures (Mann et al., 2009) and Greenland temperatures (Kobashi et al., 2011) over the past 800 yr, indicating that the general atmospheric and oceanic circulation pattern connecting Greenland and other parts of the world for the past 130 yr has been active over the past 800 yr (Fig. 2).
The temperature trend in west Iceland (Reykjavik) exhibits a similar pattern to East Greenland (Tasiilaq) over the past 130 yr, both showing a negative correlation with the NAO index (Hanna et al., 2004).The strong association of the temperature trends between East Greenland and west Iceland can be inferred from the Greenland and global correlation maps (Fig. 2).As large areas of the North Atlantic and the Summit temperatures exhibit high correlations except around the storm track region (Fig. 2), the North Atlantic average temperatures (Enfield et al., 2001) and the Summit temperatures (Box et al., 2009) are highly correlated (r = 0.58, p = 0.03 as in Fig. 3).Multidecadal variation (21 yr running means, RMs) of the Summit temperature (Box et al., 2009) can explain 70 % of the variance in the North Atlantic average temperatures (Enfield et al., 2001).Thus, the Summit temperature is highly influenced by the Atlantic multidecadal oscillation (AMO) (Enfield et al., 2001), which is a climate mode of variability with an approximate periodicity of 70 yr and has Greeland Summit temperature (°C) Summit temp.(Kobashi et al., 2011) Summit temp.(Box et al., 2009) North Atlantic average temp.(Enfield et al. 2010) Fig. 3. North Atlantic average temperatures and Greenland Summit temperatures.Error band (red) is 1σ (Kobashi et al., 2011).The Summit temperature by Box et al. (2009) was shifted by −1.75 • C to account for the surface-air temperature difference (Kobashi et al., 2011).
global impacts (Dima and Lohmann, 2007).The causes of the oscillation have been associated with the internal oscillation of the AMOC possibly influenced by solar variation and volcanic eruptions (Otterå et al., 2010;Knudsen et al., 2011).Greenland temperature variability over the past 4000 yr contains a fairly persistent periodicity of approximately 70 yr (Kobashi et al., 2011), indicating the consistent nature of the AMO over the past 4000 yr (Knudsen et al., 2011;Chylek et al., 2012).Notably, in the following sections, we identified parts of the causes of Greenland temperature over the past 4000 yr, which must be linked with the causes of variability in the AMO.

Climate forcings over the past 4000 yr
Because the residual Laurentide ice sheet from the last glacial period nearly disappeared by approximately 6000 yr BP (before present; present is defined as 1950), climate boundary conditions (e.g.sea level and areas covered by ice sheets) over the past 4000 yr were largely the same as those of the preindustrial period (Lambeck and Chappell, 2001;Kaufman et al., 2004;Wanner et al., 2008).Therefore, climate change over the past 4000 yr was primarily caused by variations in orbital, solar, volcanic, and greenhouse gas forcing as well as natural variability (Wanner et al., 2008).It is noted that for the past 4000 yr orbital forcing has become more important than during the past millennium especially in high-latitude regions.

Orbital forcing
As Greenland is located in the high latitudes of the NH, changes in the seasonal variation of insolation owing to changes in Earth's axial tilt have had a large influence on Greenland climate over the past 4000 yr (Laskar et al., 2004) (Fig. 4).Mean annual insolation at the Greenland Summit decreased by 1.4 % from 195.4 W m −2 in 2000 BCE to 192.8 W m −2 in 2000 CE (Fig. 4).At present, annual mean insolation is still decreasing at a rate of 0.44 W m −2 per 100 yr in Greenland (Laskar et al., 2004).As Greenland has no sunlight during the winter, most of the insolation change occurs in the summer season (Fig. 4) (Laskar et al., 2004).In contrast, global mean annual insolation through orbital configuration has been nearly constant over the past 4000 yr, indicating that differential changes in insolation at different latitudes compensated for each other (Laskar et al., 2004).The decrease at high-latitude insolation contrasts with that at low latitudes, which has exhibited a slight increase over the past 4000 yr (Fig. 4).It is noted that July insolation has decreased over the entire NH over the past 4000 yr, but the decrease has been compensated by an increase in winter insoation in the case of low latitudes.

Solar forcing
Changes in past solar activity can be estimated from the variation in cosmogenic nuclides, such as 14 C in tree rings (Vieira et al., 2011) and/or 10 Be concentration in ice cores (Steinhilber et al., 2009) (Fig. 5).The estimation of these changes is possible because during periods of stronger (weaker) solar activity, the intensity of the galactic cosmic rays that produce cosmogenic nuclides such as 10 Be and 14 C from nitrogen and oxygen is reduced (increased) due to stronger (weaker) solar magnetic fields (Lal and Peters, 1968;Gray et al., 2010).As both 10 Be-and 14 C-based estimates have different advantages and disadvantages owing to different chemical characters (Usoskin et al., 2009), a recent reconstruction of total solar irradiance (TSI) (Steinhilber et al., 2012) utilized available 10 Be (from Greenland and Antarctic ice cores) and 14 C data to produce "combined" estimates of TSI.Three TSI estimates from " 14 C-based" (Vieira et al., 2011), " 10 Be-based" (from Greenland ice cores) (Steinhilber et al., 2009), and "combined" methods (Steinhilber et al., 2012) exhibit a similar trend over the past 4000 yr (Fig. 5), but the combined estimate exhibits the least multidecadal variation as it captures only common variance in 14 C and 10 Be, which likely provide the most robust estimate of past solar activity (Steinhilber et al., 2012).We note that the correlation between the "combined" and " 14 C-based" estimates is higher (r = 0.87, p < 0.01) than that (r = 0.67, p < 0.01) between the "combined" and " 10 Bebased" estimates, likely because the local climate (e.g.high accumulation rate) in Greenland affected 10 Be deposition.Therefore, we primarily used the combined TSI estimates (Steinhilber et al., 2012) for the following analyses and explored uncertainties using 14 C-and 10 Be-based estimates.
During the last 4000 yr, total solar irradiance has exhibited variability in a range (maximum-minimum) of 1.23 W m −2 or approximately 0.1 % and had a slight decreasing trend of approximately 0.1 W m −2 per 1000 yr (Steinhilber et al., 2012) (Fig. 5).There are several periods (approximately 1500, 800, 400 BCE and 600, 1000, 1300, 1500, 1600, 1850 CE) of weaker solar activity similar to the Maunder Minimum, some of which reportedly caused cooling events in widespread areas (Mayewski et al., 2004;Wanner et al., 2008Wanner et al., , 2011)).The Little Ice Age and the 2.8 ka event are well-known cooling periods of the past 4000 yr (Renssen et al., 2006;Wanner et al., 2008).The standard deviation of radiative forcing ( F ) by TSI variation from 2000 BCE to 1750 CE is calculated to be 0.07 (W m −2 ) from the following equation: F = TSI/4 × 0.7 (Forster et al., 2007).We note that, for the 1-D energy balance model (EBM) calculation, average values of albedo in each latitudinal band were applied to calculate solar forcing.

Greenhouse gas forcing
Important long-lived greenhouse gases in the atmosphere are CO 2 , CH 4 , and N 2 O (Forster et al., 2007).Although concentrations of these gases in the atmosphere have  (Mayewski et al., 1997).Black and red lines are raw data and 51 yr RMs for reconstructed volcanic signals, respectively.Blue shaded areas are periods when volcanic forcings were lower than the standard deviation (−0.15 W m −2 ; dotted line) of the 51 yr RMs and lasted for more than 50 yr.Seventeen volcanically active periods were identified (Table 1).been rapidly increasing due to human activity since the industrialization (∼ 1850 CE), they can also vary naturally, but with much smaller amplitudes (Forster et al., 2007).Past variation in these gases can be obtained from ice cores (Fig. 5).We used the data for the past 2000 yr compiled by F. Joos for PMIP3 (Schmidt et al., 2011) (https://pmip3.lsce.ipsl.fr/wiki/doku.php/pmip3:design:lm:final).For older periods, lower resolution data were used for CO 2 (Monnin et al., 2004), CH 4 (GISP2 data) (Brook, 2009), and N 2 O (Flückiger et al., 2002).Slight adjustments were made to avoid gaps between CE and BCE by subtracting the differences in average concentrations in the overlapping period from 1 to 12 CE from the values used for BCE.The concentration data were calibrated to radiative forcing, F (W m −2 ), according to the following equations (Ramaswamy et al., 2001).For CO 2 , F = α ln(C/C 0 ), where α is 5.35, and , where α is 0.12 and N is N 2 O in ppb.The function f is described as f (M, N ) = 0.47 ln[1 + 2.01 × 10 −5 (MN) 0.75 + 5.31 × 10 −15 M(MN ) 1.52 ].The subscript 0 denotes the concentrations in the year 1750.Standard deviations of F for CO 2 , CH 4 , and N 2 O for the period from 2000 BCE to 1750 CE are 0.06 (W m −2 ), 0.02 (W m −2 ), and 0.01 (W m −2 ), respectively, indicating that CO 2 variation had the largest control on temperature changes among the greenhouse gases during this period.The standard deviation of total F from the three greenhouse gases from 2000 BCE to 1750 CE is 0.08 (W m −2 ).There is a slight increasing trend in total greenhouse gas forcing, F of 0.1 W m −2 per 1000 yr from 2000 BCE to 1750 CE (Fig. 5).

Volcanic forcing
Large volcanic eruptions are known to cause climate change by injecting chemically and microphysically active gases and aerosols into the atmosphere, which affect Earth's radiative balance (Robock, 2000).Sulfate data in ice cores provide the best record of past volcanic events; as a result, there have been many attempts to estimate stratospheric sulfate loading and optical depth (Zielinski et al., 1994;Zielinski, 1995;Cole-Dai et al., 2000;Gao et al., 2008).We used the GISP2 sulfate record for reconstructions of volcanic forcing (Appendix A) as it is the only currently available record that has high resolution and a good chronology with continuous biannual resolution (averages of two years) over the past 4000 yr (Zielinski et al., 1994;Mayewski et al., 1997).Multi-core analysis indicated that volcanic debris was fairly evenly distributed in the atmosphere after volcanic eruptions (Gao et al., 2007), which is an important factor in inferring the magnitudes of past volcanic eruptions from a single ice core.The standard deviation of the volcanic forcing over the period from 2000 BCE to 1750 CE is 0.15 W m −2 , which is approximately 2 times larger than the standard deviations of greenhouse gas and solar forcings over the same period (Fig. 5).In addition, we included aerosol forcing from 1892 CE onward (Crowley, 2000) (Fig. 5).
Over the past 4000 yr, we identified 17 volcanically active periods represented by a single large eruption or smaller but numerous eruptions (Fig. 5 and Table 1), and it is clear that volcanic events were most frequent during the Little Ice Age.In the NH, eight volcanically active periods can be identified during the first two millennial BCE, during which one or several large volcanic events had impacts on the climate (Fig. 5, Table 1; see later discussion).In contrast, Southern Hemisphere records exhibit very few volcanic events during the first two millennia BCE (Cole-Dai et al., 2000;Wanner et al., 2008).In terms of volcanic events, 0-600 CE was a rather calm period (Fig. 5).The largest volcanic event during the past 4000 yr occurred at 1258-1259 CE.Climate responses are somewhat ambiguous for this event as tree rings (Mann et al., 2012) and δ 18 O ice in Greenland (Zielinski, 1995) did not show clear signals of climate change, perhaps owing to the character of proxy (Mann et al., 2012) or lower climate sensitivity to volcanic events during the warmer Medieval Warm Period (Zielinski, 1995).

Integrating climate forcings and application of 1-D EBM
As five climate forcings (orbital, solar, volcanic, greenhouse gases, and aerosols) were quantified as F from 2000 BCE to 1950 CE, it is possible to integrate them and calculate the corresponding surface temperatures using climate models.Total forcing, excluding orbital forcing, shows a slight increasing trend from 2000 BCE to approximately 200 CE, due to increasing greenhouse gas forcing and fewer volcanic eruptions (Fig. 5), and remains relatively high during the first millennium CE.Then, total forcing decreases, owing to more volcanic eruptions and weaker solar activity toward 1850 CE.Finally, largely owing to increasing human-induced greenhouse gases, total forcing increases rapidly to the highest level in the past 4000 yr (Fig. 5).According to the calculated forcings, the multidecadal to centennial variations were caused primarily by volcanic forcing and secondarily by solar and greenhouse forcing during the pre-industrial period (Fig. 5).
To investigate the impacts of orbital forcing at different latitudes, we utilized a one-dimensional energy balance model (1-D EBM; Budyko-Sellers model) (Budyko, 1969;Sellers, 1969) that can be resolved for 18 latitudinal bands (10 • for each band) over the globe (McGuffie and Henderson-Sellers, 2005).The heat transports into and out of the bands were calculated from the deviations of zonal mean temperatures from global mean temperatures.Observed average albedo for each latitudinal band was used to calculate the energy balance at the latitudinal bands, which were set to be constant through time.EBMs are computationally efficient and calculate mean annual hemispheric temperatures close to the observations (Crowley, 2000;Mann et al., 2012).The 1-D EBM was tested for seasonal heat transport divergence and found to match observations, especially for high-latitude regions (Warren and Schneider, 1979).
Latitudinally variable solar insolation (for example, annual insolation at 65 • N was used for a latitudinal band of 60 • N −70 • N) was obtained from Laskar et al. (2004).Then, fractions of solar insolation from the present value for each latitudinal band were calculated for every year, and applied to all latitudinal bands for each time step (1 yr).Although better forcing estimates are available toward the more recent past, we used the climate forcings described above for the entire period to verify the fidelity of the forcing estimates by comparing them with available temperature records.The model simulation starts at 2200 BCE and calculates until 1975 CE.Calculated temperatures for a latitudinal band of 70-80 • N exhibited a secular decreasing trend over the past 4000 yr owing to orbital forcing (Fig. 5).It contrasts with calculated NH average temperatures with little slope, but the calculated NH temperature change closely resembles the total forcing variation (r = 0.97, p < 0.01) (Fig. 5e and f).This indicates that the influence of orbital forcing on NH annual average temperatures is limited over the past 4000 yr, which is consistent with equilibrium runs for the mid-Holocene (Braconnot et al., 2007) and transient runs (Renssen et al., 2006) of GCM (general circulation model) simulations.The model-derived 70-80 • N temperatures exhibit warmer temperatures than the Greenland Summit temperatures because of the high altitude of the Summit (3200 m a.s.l.).However, average temperatures of the model-derived 70-80 • N temperatures (Fig. 5) are similar to the mean annual temperature (−11.1 ± 2.3 • C) of the north-west coast of Greenland (Pituffik/Thule Air Base; 76 • 32 N, 68 • 45 S) (Box, 2002).

Causes of Greenland temperature variability over the past 4000 yr
Multidecadal to centennial Greenland temperature variations generally follow NH temperature changes with polar amplification, but deviate negatively (positively) when solar output increases (decreases) owing to atmospheric/oceanic changes (e.g.NAO/AO and AMOC) over the past 800 yr (Kobashi et al., 2013).Renssen et al. (2006) also proposed a mechanism to explain the negative responses of Greenland temperatures to solar variability, indicating that during weaker (stronger) solar activity, formation of deep water in the Nordic seas relocates toward (from) Iceland because of increased (decreased) sea ice in the Nordic seas, which warms (cools) south-east Greenland due to the large heat release from ocean in the convective areas.The Greenland temperature deviation from NH average temperature trend is defined as "Greenland temperature anomaly (GTA)" (Kobashi et al., 2013).As the background climate over the past 4000 yr is largely similar to the past 800 yr (Wanner et al., 2008), it can be assumed that the negative response of Greenland temperature to solar variability continued over the past 4000 yr.Therefore, the relationship can be written as follows: where notations are denoted as follows: k [unitless] is a scaling coefficient, GTA [unitless] the Greenland temperature anomalies, GT [ • C] Greenland temperature, NHT [ • C] NH average temperature or 70-80 • N average temperature, TSI [W m −2 ] total solar irradiance (linearly detrended), X average of a time series denoted by subscript, and s the standard deviation of a time series denoted by subscript.
First, we calculated the GTA using model estimates of 70-80 • N average temperatures (Fig. 5, bottom panel) and the ice-core-derived Greenland temperature (Kobashi et al., 2011) to account for the influence of orbital forcing.The calculated GTA over the past 4000 yr significantly correlated with detrended TSI (r = −0.41,p < 0.01), and the k value was found to be 1.16.However, it is noted that the calculated NH temperature contains solar signals as a forcing, such that the calculated GTA and TSI are not totally independent.As the variation of solar forcing is relatively small compared with volcanic forcing, we computed a GTA without solar forcing, and found a weak but significant correlation between GTA and TSI (r = −0.21,p = 0.05).The low correlation is expected as we do not know true NH average temperature with internal variation.For the past 800 yr, the natural variation was found to be rather large (Kobashi et al., 2013).Another useful finding here is that the calculated GTA had no long-term trend owing to nearly identical slopes between the modelled 70-80 • N and the ice-core-derived Greenland temperatures, indicating that solar influence on Greenland temperature is confined to multidecadal to millennial changes.Therefore, we used linearly detrended TSI for the following analysis.
We calculated a history of Greenland temperature (Fig. 6) from the model outputs of 70-80 • N average temperatures using the relationship in Eq. ( 2).The scaling coefficient k was estimated to be 1 ± 0.5 (2σ ) by calculating correlation coefficient with various k values (Figs. 6 7).The modelled and observed Greenland temperatures agree with correlation coefficients of r = 0.34-0.36(p = 0.1-0.04) in a multidecadal timescale (21 yr running means; RMs) and r = 0.38-0.45(p = 0.1-0.05) in a centennial timescale (101 yr RMs).Therefore, the model with the reconstructing climate forcings can explain 14 to 20 % of the variance of the observed temperature in multidecadal to centennial timescales with a 90-96 % confidence interval (Fig. 7).This is smaller than 25% for European winter and spring temperatures, or 35 % for annual temperatures in 11 yr smoothed data over the past five centuries (Hegerl et al., 2011).Unexplained variances could come from uncertainties in estimation of climate forcings, internal variability of climate system on global and regional scales (Kobashi et al., 2013), and uncertainties in the ice-core-derived Greenland temperatures.The result also confirms that weak but persistent solar influence on Greenland temperature has continued over the past 4000 yr.The k values of 1 ± 0.5 indicates that about 50 % of the explained multidecadal to centennial Greenland temperature variability originates from negative temperature response of Greenland temperature to solar variation, which provides important implications on current and future Greenland temperatures (Kobashi et al., 2013).We note that neither modelled 70-80 • N average temperature nor TSI had significant correlations with the ice-core-derived Greenland temperature by itself (Fig. 7).
Figure 6 shows the model-based Greenland temperature with various k values, comparing with the ice-core-derived Greenland temperature.It is rather clear that the long-term cooling trends of the ice-core-derived Greenland temperatures are well-reproduced by the model output of 70-80 • N average temperature (Fig. 6; see the panel for k = 0), and  (Kobashi et al., 2011) and modelled (blue) Greenland temperatures with 2σ error bands.We calculated Greenland temperatures from the model outputs of 70-80 • N average temperatures using the relationship in Eq. ( 2).The scaling coefficient k is estimated to be 1 ± 0.5 (2σ ).X GT and s GT are from the ice-core-derived Greenland temperature (Kobashi et al., 2011).Uncertainties for TSI in Eq. ( 2) were derived from the standard deviation of three detrended and normalized TSI records (Fig. 5c).Then, uncertainties in the modelled Greenland temperature were calculated by the theory of error propagation.Black line at the top panel is the Greenland Summit temperature reconstruction from the GRIP borehole temperature profile (Dahl-Jensen et al., 1998).Lower panels: modelled Greenland temperatures with different k values, compared to the ice-core-derived Greenland temperature (Kobashi et al., 2011).several large volcano events seems to correspond to cooling in the Greenland temperature, especially around 1600 and 700 BCE, and during the Little Ice Age (Fig. 6, k = 0).Model outputs with k values around 1 produce more multidecadal to centennial temperature variations that better fit with the multidecadal to centennial variation of the ice-core-derived Greenland temperatures than that of k = 0.This is especially true for the warming events around 1500, 1000, and 300 BCE, and 700 CE (Fig. 6).The k value around 1.5 apparently creates too large centennial variation compared to that of the ice-core-derived Greenland temperature (Fig. 6, bottom panel).In addition, Fig. 6 clearly demonstrates that without solar influence (k = 0.0) Greenland temperature should have been much warmer from the 20th century to present.We note that three independent estimates of Greenland temperatures over the past 4000 yr (borehole-based; Dahl-Jensen et al., 1998, nitrogen/argon isotopes-based;Kobashi et al., 2011, and model) exhibit a similar millennial trend of higher temperatures in 2000-1000 BP, cooler temperatures in 500 BC-1 CE, and warmer temperatures in 500-1000 CE (Medieval Warm Period) and 1500-1800 CE (the Little Ice Age).Part of the millennial variation is clearly caused by negative Greenland temperature responses to solar variability (Fig. 6).
6 Implications for NH temperature change over the past 4000 yr

Northern high-latitude temperature change
From the relationship between Greenland and NH temperatures and the understanding from the model experiments, it is plausible to gain insights on northern hemispheric temperature changes from the Greenland temperature for the past 4000 yr.A NH temperature history was calculated from the ice-core-derived Greenland temperatures and detrended TSI following Eq.(3) (Fig. 8).Uncertainty ranges (2σ ) are estimated by the same way for the Greenland temperature with the k value of 1.0 ± 0.5 (see caption for Figs.6 and 8).
As the calculated record is likely influenced by orbitally induced cooling of northern high latitudes, we call them northern high-latitude (NHL) temperatures (Fig. 8, top panel).
Long-term slopes of the NHL temperature and model outputs of 70-80 • N average temperature agree well, indicating that heat transport between latitudes are well constrained in the model (Fig. 8, top panel).The modelled 70-80 • N average temperature and NHL temperature are significantly correlated (r = 0.53, p < 0.01) over the past 4000 yr (Fig. 8, bottom panel).As the two time series are not completely independent as both estimates possess the same TSI signals, we tested the relationship by calculating a modelled 70-80 • N average temperature without solar forcing as we did for the GTA.It was found that the model-derived 70-80 • N temperature without solar forcing was significantly correlated (r = 0.33, p = 0.03) with the NHL temperature, noting that now the two time series are totally independent.Next, we compared the NHL temperature with a multiproxy polar region temperature record (60-90 • N; primary summer temperature) (Kaufman et al., 2009) for the past 2000 yr (Fig. 8, top panel).As the proxy temperatures are known to record spring to summer temperatures (Kaufman et al., 2004(Kaufman et al., , 2009)), this comparison provides an opportunity to test a hypothesis if multi-proxy temperature estimates of polar region reflect temperature changes without seasonal biases.Correlations between the NHL temperature and the polar region temperature are remarkably high with r = 0.60 (p = 0.01) for the past 2000 yr (Fig. 8, top panel).Therefore, the multi-proxy temperature reconstruction seems to reflect unbiased multidecadal average temperatures at least over the past 4000 yr.It is noted that correlation between the proxy polar-region temperatures and the ice-core-derived Greenland temperatures is lower and barely significant (r = 0.33, p = 0.11), and correlation between the proxy polar-region temperature and TSI is significant but lower (r = 0.40, p = 0.04), indicating that the Greenland temperature plus negative solar signals creates the highest and most significant correlation.In addition, the model outputs of 70-80 • N average temperatures with all forcing exhibited significant correlation (r = 0.62, p = 0.02) with the polar region temperature over the past 2000 yr, indicating that 38 % Fig. 8. Modelled and ice-core-derived northern high-latitude and NH temperature anomalies (from the base year of .Top panel: northern high latitudes (NHL) temperatures.Red line is the NHL temperature derived from the ice-core-derived Greenland temperature (Kobashi et al., 2011) with 2σ uncertainty band.Green line is a multi-proxy polar temperature record in 21 yr RMs (Kaufman et al., 2009).Blue and purple lines are the model outputs for 70-80 • N average temperatures with all forcings and with all but excluding greenhouse gas forcing over the past 4000 yr, respectively.Black line is observed NH temperature anomalies in 21 yr RMs for the period from 1860 to 2000 CE (HadCRUT3) (Brohan et al., 2006).Red and blue circles indicate the temperatures with all and with all excluding greenhouse gas forcing at 1976 CE (the end year of model calculation), respectively.Note that left and right y axes have the same scale but with different absolute values.The NHL temperature from the ice-core-derived Greenland temperature was derived from the relationship in Eq. ( 3) with k = 1 ± 0.5, and s NHT and X NHT are from the observed NH average temperature in 21 yr RMs (HadCRUT3) (Brohan et al., 2006).Bottom panel: reconstructed NH temperature anomaly.Blue line is the NH average temperature estimated from the Greenland temperature record with the 2σ uncertainty band.The record was calculated by adjusting the slope of the NHL temperature to be equal to the slope of the modelled NH temperature for the period of 2000 BCE-1850 CE.Black line is the observed NH average temperatures (HadCRUT3) (Brohan et al., 2006) in 21 yr RMs. Green and purple lines are two proxy NH average temperatures by Mann et al. (2008) and Moberg et al. (2005) in 21 yr RMs, respectively.The black dotted line is the present multidecadal temperature anomaly (0.42 • C) from 1990 to 2010 CE (HadCRUT3) (Brohan et al., 2006).
of the variance can be explained by the model (Fig. 8, top panel).
The top panel in Fig. 8 also shows the model outputs of 70-80 • N average temperatures with all forcings (Fig. 5e) and with all forcing but excluding greenhouse gas forcing to illustrate the roles of greenhouse gases over the past 4000 yr.From the two outputs, it can be concluded that greenhouse gas forcing played two critical roles over the past 4000 yr for the temperature changes.One is that it caused cooler temperatures in earlier periods of the past 4000 yr than temperatures without greenhouse gas forcing owing to slowly increasing greenhouse gas forcing from 2000 BCE to 1850 CE (Fig. 5d).The second is that it caused a rapid warming during the 20th century (Fig. 8).Without the increasing greenhouse gas forcing in the 20th century, the Little Ice Age may have continued into present (Fig. 8, top panel).This reinforces earlier findings (Solomon et al., 2007) that an increase in the NH temperature during the 20th century cannot be explained without the human-induced greenhouse gas increases.

Northern hemispheric average temperature change
The model outputs indicate that 70-80 • N and NH average temperature responses to climate forcings are different primarily in their long-term trends (Fig. 5, bottom panel).Over the past 4000 yr, low latitudes experienced increasing mean annual insolation, and high latitudes experienced decreasing mean annual insolation (Fig. 4).As low latitudes have Normally distributed random numbers (mean zero and unit standard deviation) were generated, and multiplied with 1σ uncertainties of the time series.Then, it was added to the original time series to generate 100 synthetic time series of each NH and NHL temperature anomaly for the period 2000 BCE to 1850 CE.Then, histograms were created with all the synthetic data plus the original data, and then the values for each bin were divided by 101.Future projections from IPCC AR4 (Meehl et al., 2007) were originally referenced to 1980-1999, but they were adjusted to the reference to 1961-1990 by adding 0.19 • C, i.e. the difference of the NH average temperature anomaly between 1980-1999and 1961-1990 (HadCRUT3) (HadCRUT3) (Brohan et al., 2006).
larger areas than high latitudes, NH average temperature generally follows low-latitude temperature trends under conditions of a constant background climate.As the model does a good job of calculating realistic heat transports between latitudes, we assume that the calculated slope of the NH average temperature is reasonably constrained (Fig. 5, bottom panel).Therefore, we adjusted the slope of the NHL temperature to be equal to the slope of the modelled NH average temperature for the period of 2000 BCE to 1850 CE to estimate an ice-core-derived NH average temperature change (Fig. 8, bottom panel).Then, we compared the NH temperatures with two proxy NH average temperature records (Fig. 8, bottom panel).The ice-core-derived NH temperature significantly correlated with the proxy NH temperatures (Moberg et al., 2005;Mann et al., 2009) for the overlapping periods (Fig. 8 The poorer correlations in the longer periods likely reflect poor quality of the NH proxy temperature as the two proxy NH temperature estimates diverge especially before 1000 CE.We note that two proxy records in 21 yr RMs did not show significant correlations with the detrended TSI (r = 0.22, p = 0.24 for Moberg for the period from 11 to 1969 CE; r = 0.3, p = 0.12 for Mann for the period from 310 to 1994 CE).
The reconstructed NH average temperature from the ice core is the most useful before 2000 yr ago as no proxy estimates of NH temperature in multidecadal to centennial scales reach that far back in time.We found that the NH temperature experienced two sequential cooling events at approximately 800 and 400 BCE (the 2.8 ka event), which had similar magnitudes to the Little Ice Age although the duration was shorter.The 2.8 ka event has been thought to be global (Chambers et al., 2007), and reportedly caused by decreasing solar outputs (Chambers et al., 2007;Martin-Puertas et al., 2012).However, we found that large volcanic eruptions during the events played an equally important role for the cooling as weaker solar forcing (Fig. 5d).The 2.8 ka event corresponds to the Iron Age Cold Epoch in the European Mediterranean climate zone.We note that Greenland did not experience the cooling as much as the NH temperature because of negative response to solar variability.Renssen et al. (2006) also found in their model and climate proxies that a warm inflow in the northern Icelandic shelf influenced areas near Iceland and East Greenland during the event such that Greenland did not experience severe cooling as other parts of the North Atlantic basin.
Lastly, Fig. 9 shows distribution of the ice-core-derived NH and NHL temperatures for the period from 2000 BCE to 1850 CE.It indicates that the current multidecadal average temperature (1990-2010) and/or higher temperatures occurred at a frequency of 1.3 times per 1000 yr for the northern high-latitude temperature and 0.36 times per 4000 yr for the NH temperatures.Therefore, we conclude that it is more likely than not (the 64 % significance level) that the current multidecadal NH average temperature (1990-2010) is unprecedented over the past 4000 yr.A recent study indicated that the present NH average temperature (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) has not exceeded the warm early Holocene temperature (Marcott et al., 2013), but our study indicates that the multidecadal average temperature (1990-2010) did exceed the temperatures at least for the past 4000 yr.Future model projections (Meehl et al., 2007) with increasing anthropogenic greenhouse gases show that future decadal NH average temperatures (2090-2099) of all the projections will be completely out of the range of the NH average temperature over the past 4000 yr (Fig. 9).

Conclusions
The Greenland temperatures over the past 4000 yr reconstructed from trapped air within the GISP2 ice core (Kobashi et al., 2011) provided an extraordinary opportunity to investigate the late Holocene climate changes because of several advantages: (1) the resolving of precise multidecadal to millennial temperature variation; (2) the recording of "mean" annual temperatures (many palaeotemperatures are spring to summer proxies); (3) the tight age control; (4) the understanding of regional climate (Kobashi et al., 2013); and (5) the plentiful palaeoclimate information that is available from the GISP2 ice core.With the reconstructed climate forcings (orbital, volcanic, solar, greenhouse gases, and aerosols) over the past 4000 yr, we calculated northern high-latitude and NH average temperatures with 1-D EBM.Then, modelled Greenland temperature was derived considering negative Greenland temperature responses to changes in solar output.Calculated Greenland temperatures were significantly correlated with the ice-core-derived Greenland temperatures, explaining 14 to 20 % of variance of the Greenland temperature variability in multidecadal to centennial timescales with a 90-96 % confidence interval.
Multidecadal to centennial NH temperature variations beyond the last millennium were shown to be difficult to constrain owing to various reasons such as large seasonal variations of climate variables, chronological issues, and changes in biological dependences on climate variables.In this study, we proposed a new way of investigating the multidecadal to centennial NH temperature variation with the climate model and the Greenland temperature record (Kobashi et al., 2011) that provides seasonally unbiased estimates of temperature change.It indicated that current multidecadal NH temperature (1990-2010) is more likely unprecedented than not (p = 0.36) over the past 4000 yr.

Reconstruction of volcanic forcing from GISP2 sulfate record
The GISP2 sulfate record has been used to estimate volcanic events over the past 9000 yr (Zielinski et al., 1994), as well as stratospheric loading and optical depth from volcanic eruptions over the past 2100 yr (Zielinski, 1995).More recently, volcanic events and their associated stratospheric loading over the past 1500 yr were reconstructed using multiple cores from both the polar regions of Greenland and the Antarctic (Gao et al., 2008).Although multi-core estimation is clearly more reliable, we attempted a single-core reconstruction with GISP2 due to the lack of multi-core estimation for the entire period, and evaluated the results with multicore estimates for an overlapping period.One caveat is that a single-core reconstruction cannot differentiate latitudinal variations in sulfate loading, as sulfate loading from volcanic eruptions near Greenland (e.g.Iceland) is inevitably larger.This is presumably less problematic when a single-core reconstruction is used as volcanic forcing for northern high latitudes, but care must be taken to interpret it as a hemisphericscale forcing.
To estimate stratospheric loading over the past 4000 yr, we extracted volcanic signals from the sulfate ion record of the GISP2 ice core (Mayewski et al., 1997), largely following the methodology developed by Gao et al. (2006Gao et al. ( , 2007Gao et al. ( , 2008)).To be consistent with the timescale of the temperature reconstruction (Kobashi et al., 2008b(Kobashi et al., , 2011)), we applied a GISP2 timescale, called "visual stratigraphy" (Alley et al., 1997b), to the GISP2 sulfate record.Therefore, age uncertainties between the temperature and reconstructed volcanic forcing are at a minimum, which is important as the climatic impacts of volcanic eruptions last for only a few years (Robock, 2000).It is noted that the GISP2 sulfate record is reliable, as it matches quite well with other Greenland core data (unpublished data).
First, sulfate ion data (ppb) (Mayewski et al., 1997) from the GISP2 ice core were converted to sulfate flux using accumulation rate data (Alley et al., 1997b).Then, a high-pass loess filter (Cleveland and Devlin, 1988) was applied to remove signals longer than 31 yr (Gao et al., 2008).We only extracted signals larger than the running median absolute deviation (RMAD) multiplied by 0.5.We chose 0.5 instead of 2, as in the original methodology by Gao et al. (2008), to extract slightly more volcanic events over the past 1500 yr, because peaks from a single core may miss some of the important volcanic signals detected in multi-core analysis.Although it is still possible to miss relatively minor volcanic signals or to identify incorrect signals, our results should not be affected because these signals are small and thus climatically less significant.With this procedure, 72 peaks were extracted from the GISP2 sulfate record, compared to 63 peaks GISP2 (this study) Multi-core (Gao et al. 2008) Figure A1.
Volcanic forcing (W/m 2 ) Fig. A1.Reconstructed volcanic forcing over the past 1500 yr on the multidecadal timescale.Red and blue lines represent volcanic forcing reconstructions from a single core (GISP2, this study) and multi-cores (Gao et al., 2008).51 yr RMs were applied for both data to reduce noises.
in Gao et al. (2008), for NH volcanic events for the past 1500 yr.
To calibrate the selected GISP2 volcanic signals to stratospheric volcanic sulfate aerosol loading (T ), we summed the volcanic signals from GISP2 over the past 1500 yr and then scaled it to the sum of the stratospheric volcanic sulfate aerosol loading (T ) for the NH by Gao et al. (2008).Then, the scaling factor was used to reconstruct the entire volcanic record over the past 4000 yr.It is assumed that sulfate lifetime in the stratosphere has an e-holding time of 1 yr (Crowley, 2000).As we were investigating the impacts of volcanic events on multidecadal to centennial temperature changes, a 51 yr running mean (RM) filter was applied to the time series to reduce noise (higher frequency variations agree less with the multi-core forcing estimate).Although for the investigation of multidecadal to centennial climate variation the filtering should not be a problem, one needs to take caution when comparing the volcanic forcing with higher resolution data such as those from tree rings.The reconstructed stratospheric volcanic sulfate aerosol loading data from the single core and multiple cores in 51 yr RMs agree well, with a correlation coefficient of r = 0.68 (p < 0.01) over the past 1500 yr (Fig. A1), indicating that a single-core estimate can reconstruct NH volcanic forcing reasonably well.Climate forcing ( F ) of the calculated stratospheric volcanic sulfate aerosol loading was derived by converting the sulfate loading to aerosol optical depth (AOD; = aerosol loadings/150 T ; Stothers, 1984).Then, the AOD was multiplied by −20 W m −2 to obtain F (Wigley et al., 2005).
et al.: Causes of Greenland temperature variability over the past 4000 yr

Fig. 5 .
Fig. 5. Climate forcings and 1-D EBM output over the past 4000 yr.(a) Greenhouse gas forcings by CO 2 , CH 4 , and N 2 O from 2000 BCE to 2000 CE.The values at 1750 CE were defined as zeros.See text for the data sources.(b) Volcanic forcing for NH over the past 4000 yr derived from the GISP2 sulfate record(Mayewski et al., 1997).Black and red lines are raw data and 51 yr RMs for reconstructed volcanic signals, respectively.Blue shaded areas are periods when volcanic forcings were lower than the standard deviation (−0.15 W m −2 ; dotted line) of the 51 yr RMs and lasted for more than 50 yr.Seventeen volcanically active periods were identified (Table1).(c) TSI anomaly over the past 4000 yr from the value of 1365.57W m −2 in 1986.Data are plotted on the standardized unit (mean zero and standard deviation of one) to show differences in trend but not in absolute values.The data for combined, 14 C-based, and 10 Be-based estimates are fromSteinhilber et al. (2012),Vieira et al. (2011), andSteinhilber et al. (2009), respectively.Left y axis is for the combined TSI(Steinhilber et al., 2012), which was used in the 1-D EBM.Shaded areas indicate periods of weaker solar activity.(d) Climate forcings from 2000 BCE to 1976 CE. (e) Total climate forcings.This is sum of various forcings (d).(f) Modelled NH and 70-80 • N temperatures using the 1-D EBM with the total climate forcing as in (e) and the orbital forcing in the Fig. 4 (see text).
Fig. 5. Climate forcings and 1-D EBM output over the past 4000 yr.(a) Greenhouse gas forcings by CO 2 , CH 4 , and N 2 O from 2000 BCE to 2000 CE.The values at 1750 CE were defined as zeros.See text for the data sources.(b) Volcanic forcing for NH over the past 4000 yr derived from the GISP2 sulfate record(Mayewski et al., 1997).Black and red lines are raw data and 51 yr RMs for reconstructed volcanic signals, respectively.Blue shaded areas are periods when volcanic forcings were lower than the standard deviation (−0.15 W m −2 ; dotted line) of the 51 yr RMs and lasted for more than 50 yr.Seventeen volcanically active periods were identified (Table1).(c) TSI anomaly over the past 4000 yr from the value of 1365.57W m −2 in 1986.Data are plotted on the standardized unit (mean zero and standard deviation of one) to show differences in trend but not in absolute values.The data for combined, 14 C-based, and 10 Be-based estimates are fromSteinhilber et al. (2012),Vieira et al. (2011), andSteinhilber et al. (2009), respectively.Left y axis is for the combined TSI(Steinhilber et al., 2012), which was used in the 1-D EBM.Shaded areas indicate periods of weaker solar activity.(d) Climate forcings from 2000 BCE to 1976 CE. (e) Total climate forcings.This is sum of various forcings (d).(f) Modelled NH and 70-80 • N temperatures using the 1-D EBM with the total climate forcing as in (e) and the orbital forcing in the Fig. 4 (see text).

Fig. 6 .
Fig. 6.Modelled and observed Greenland temperatures with various k values.Top panel: ice-core-derived (red)(Kobashi et al., 2011) and modelled (blue) Greenland temperatures with 2σ error bands.We calculated Greenland temperatures from the model outputs of 70-80 • N average temperatures using the relationship in Eq. (2).The scaling coefficient k is estimated to be 1 ± 0.5 (2σ ).X GT and s GT are from the ice-core-derived Greenland temperature(Kobashi et al., 2011).Uncertainties for TSI in Eq. (2) were derived from the standard deviation of three detrended and normalized TSI records (Fig.5c).Then, uncertainties in the modelled Greenland temperature were calculated by the theory of error propagation.Black line at the top panel is the Greenland Summit temperature reconstruction from the GRIP borehole temperature profile(Dahl-Jensen et al., 1998).Lower panels: modelled Greenland temperatures with different k values, compared to the ice-core-derived Greenland temperature(Kobashi et al., 2011).

Table 1 .
Seventeen periods of strong volcanic events at a multidecadal scale over the past 4000 yr.Positive and negative years indicate CE and BCE, respectively.Average forcing has a unit of W m −2 .Note that the dates were calculated from smoothed volcanic forcing such that a comparison with proxy records is effective only for multidecadal or longer timescales.