Coherence among the Northern Hemisphere land , cryosphere , and ocean responses to natural variability and anthropogenic forcing during the satellite era

A lack of long-term measurements across Earth’s biological and physical systems has made observation-based detection and attribution of climate change impacts to anthropogenic forcing and natural variability difficult. Here we explore coherence among land, cryosphere and ocean responses to recent climate change using 3 decades (1980–2012) of observational satellite and field data throughout the Northern Hemisphere. Our results show coherent interannual variability among snow cover, spring phenology, solar radiation, Scandinavian Pattern, and North Atlantic Oscillation. The interannual variability of the atmospheric peak-totrough CO2 amplitude is mostly impacted by temperature-mediated effects of El Niño/Southern Oscillation (ENSO) and Pacific/North American Pattern (PNA), whereas CO2 concentration is affected by Polar Pattern control on sea ice extent dynamics. This is assuming the trend in anthropogenic CO2 emission remains constant, or the interannual changes in the trends are negligible. Our analysis suggests that sea ice decline-related CO2 release may outweigh increased CO2 uptake through longer growing seasons and higher temperatures. The direct effects of variation in solar radiation and leading teleconnections, at least in part via their impacts on temperature, dominate the interannual variability of land, cryosphere and ocean indicators. Our results reveal a coherent long-term changes in multiple physical and biological systems that are consistent with anthropogenic forcing of Earth’s climate and inconsistent with natural drivers.


718
A. Gonsamo et al.: Coherence among the Northern Hemisphere land, cryosphere, and ocean responses ity between models and the stochastic climate system with limited knowledge on processes and mechanisms involved in external forcing and climate model response, make disentangling the relative roles of natural variability and anthropogenic forcing challenging (e.g., Fyfe et al., 2013;Hegerl and Zwiers, 2011).
Observation-based causal attribution analysis of recent biosphere responses to climate change (e.g., Parmesan and Yohe, 2003;Parmesan, 2006;Walther et al., 2002;Kelly and Goulden, 2008;Post et al., 2013;Wu et al., 2011;Menzel et al., 2006;Montes-Hugo et al., 2009) is further complicated because of a lack of long-term observational data across life supporting natural systems to attribute the detected climate change impacts to natural variability and anthropogenic forcing. Although the attribution of climatic conditions and detection of the anthropogenic signal is now a mature discipline going back several decades (e.g., National Research Council, 1983;Wigley and Barnett, 1990;IPCC, 2007IPCC, , 2014b, detection and attribution of climate change impacts on human and natural systems started only recently (e.g., attribution of physical and biological impacts to warming, Rosenzweig et al., 2008; methodological framework for climate change impact attribution in conservation and ecological research, Parmesan et al., 2013; conceptual framework to detect and attribute effects of climate change, Stone et al., 2013; climate change impacts on marine life, Poloczanska et al., 2013). For detection and attribution of climate change impact assessments to work, understanding of numerous drivers, which may not often have additive interactions, are required along with observational data across a broader range of Earth's systems (Stone et al., 2013;Parmesan et al., 2013). Previous studies on change detection and attribution of both climatic condition and its impact on physical and biological systems often focused on a single or few climatic (e.g., National Research Council, 1983;Wigley and Barnett, 1990;IPCC, 2007IPCC, , 2014b or response indicator variables (e.g., Rosenzweig et al., 2008;Parmesan et al., 2013;Stone et al., 2013;Poloczanska et al., 2013) with the analysis of mechanisms. However, the test for a coherent detection and attribution of impacts of climate change requires observations of consistent patterns across natural systems, including the need to synthesize information from a much broader range of disciplines. Therefore, with accumulation of satellite observations over the last 3 decades, we here synthesize data sets across several physical and biological systems to test the relative roles of natural variability and anthropogenic forcing on impacts of recent climate change. Our aim is not to explain the entire observed variances, but rather to study the relative contribution of each of the examined driving variables on interannual changes and long-term trends of physical and biological systems. Furthermore, we study the auto-association among the response and driving variables in a way that best explains the observed variances. Unexplained variances remain, and may be attributed to missing drivers, errors in data, or methodological difficulties in capturing feedbacks and short-term adaptations in Earth's natural systems. Due to a lack of highquality observational time series, the key missing drivers in our analysis include volcanism and aerosol.
Variations in incoming solar radiation have the potential to influence global climate trends (Rind, 2002;Hameed and Gong, 1994). The rates of incoming solar radiation changes due to different solar cycles and apparently produce corresponding changes in weather and climate on the Earth's surface (Lamb, 1972;Burroughs, 1992). The changes in the amount of incoming solar radiation affects climate, both directly through changes in the rate of solar heating of the Earth and atmosphere, and indirectly, by changing cloud formation (Udelhofen et al., 1999) and ozone patterns since changes in solar cycle affect the solar ultraviolet radiation the most (Hood, 1997). But the extent of responses of climatic, biological, and physical systems to solar variability remains largely untested due to a lack of long time series measurements. It has been suggested that cosmic rays induce low-level cloud formation through atmospheric ionization and during periods of low solar activity, more cosmic rays enter the Earth's atmosphere affecting the Earth's climate system (Svensmark and FriisChristensen, 1997;Svensmark, 1998;Carslaw et al., 2002). Solar and cosmic ray activities have now been monitored for the past 3 decades, and those data sets can be tested for associations with climatic dynamics and biosphere responses.
Here we present the relationships of land, cryosphere and ocean indicators to recent changes in surface temperature, greenhouse gases, internal climatic variability, solar radiation, sunspots, and cosmic rays. While longer-term analyses (e.g. century scale) can sometimes yield better statistics, we focus on the satellite era when data quality far exceeds that of earlier years. The satellite era begins in the year 1980, the starting year for continuous full global coverage observations of atmosphere, and land and ocean surfaces. All satellite data sets are spatially averaged time series to partially diminish the effects of stochastic noise. We have selected several indicators for which high-quality, long time series satellite observations, with high retrieval accuracy, covering most or all of the Northern Hemisphere are available, and relate to temperature. This is because temperature fulfils the key assumption of detection and attribution studies where the response to external forcing is a deterministic change and to first order, and signals and noise superimpose linearly (Meehl et al., 2003). We also analyse the relative roles of natural and anthropogenic driving variables on the changes of response variables with and without temperature mediation. The term "temperature mediation" refers to the impacts of natural and anthropogenic driving variables on the response variables primarily through changes in temperature (as opposed to changes in precipitation, radiation and associated variables such as cloudiness and humidity).  Barnston and Livezey (1987), and Wallace and Gutzler (1981) East Atlantic Pattern (EA) 1980-2012 Northern Hemisphere (> 20 • N) Barnston and Livezey (1987), and Wallace and Gutzler (1981) West Pacific Pattern (WP) 1980-2012 Northern Hemisphere (> 20 • N) Barnston and Livezey (1987), and Wallace and Gutzler (1981) Pacific/North American Pattern (PNA)

Observations
We use several observational satellite and field data throughout the Northern Hemisphere including temperature, cryosphere indicators, land indicators, ocean indicators, external anthropogenic forcing indicators, external natural forcing indicators, and internal climatic variability indicators (Table 1). In total we use 13 driving (external anthropogenic forcing indicators, external natural forcing indicators, and internal climatic variability indicators) and 9 response (temperature, cryosphere indicators, land indicators, and ocean indicators) variables (Table 1). The key land indicators include satellite-measured spring thaw (ST) (Barichivich et al., 2013) and start of growing season (SOS) (Barichivich et al., 2013) of the Northern Hemisphere (> 45 • N), the Point Barrow station atmospheric CO 2 concentration peak-to-trough amplitude (AMP), and field-measured first flower bloom (FFB) day of Canada (Gonsamo et al., 2013). The cryosphere indicators include satellite-measured sea ice concentration (SIC) and extent (SIE) of the Northern Hemisphere (> 31 • N), and snow cover (SC) of the Northern Hemisphere (0-90 • N). Satellite-measured sea level (SL) of the Northern Hemisphere (0-90 • N) was also included as ocean indicator to assess the response of open water bodies to climate change. The forcing and natural variability indicators include the following: Point Barrow station atmospheric CO 2 concentration (PPM); sunspot number (SP), satellite-measured solar irradiance (RAD); cosmic ray (CR) counts at Kiel sta-720 A. Gonsamo et al.: Coherence among the Northern Hemisphere land, cryosphere, and ocean responses tion; and eight leading National Oceanic and Atmospheric Administration (NOAA) Northern Hemisphere teleconnection indices: the North Atlantic Oscillation (NAO), East Atlantic Pattern (EA), West Pacific Pattern (WP), Pacific/North American Pattern (PNA), East Atlantic/West Russia Pattern (WR), Scandinavia Pattern (SCA), Polar/Eurasia Pattern (POL), and El Niño/Southern Oscillation (ENSO)-Niño 3.4 index (NINO). We have also included sets of well-mixed global greenhouse gases (WMGHG), GISS global stratospheric aerosol optical thickness at 550 nm, Atlantic Multidecadal Oscillation (AMO), and total solar irradiance to identify the relative roles of large-scale and long-term forcing and decadal internal climate variability during the satellite era. For surface temperature, we use the Goddard Institute for Space Studies (GISS) analysis of Northern Hemisphere (0-90 • N) (Hansen et al., 2010). Details and summary of each variable are given below and in Table 1, respectively.

Spring thaw
The spring thaw (ST) day of year was estimated from the daily combined freeze-thaw dynamics as the first day when at least 12 out of 15 consecutive days were classified as non-frozen (a.m. and p.m. thawed) between January and June (Barichivich et al., 2013;Kim et al., 2012;Xu et al., 2013) based on global microwave observations from morning (a.m.) and afternoon (p.m.) equatorial crossings of the Special Sensor Microwave Imager (SSM/I). The ST data set is for Northern Hemisphere (> 45 • N) for the period of 1988-2007.

Start of growing season
The start of season (SOS) day of year is calculated from the biweekly 8 km third generation (NDVI3g) normalized difference vegetation index (NDVI) data set produced from Advanced Very High Resolution Radiometer (AVHRR) observations by the Global Inventory Modeling and Mapping Studies (GIMMS) group at NASA Goddard Space Flight Center to characterize the photosynthetic growing season of the Northern Hemisphere (> 45 • N) for 1982-2011 (Xu et al., 2013). SOS was calculated from maximum rate (inflection point) of green-up as determined by the first derivative of the seasonal curve of smoothed NDVI data (Barichivich et al., 2013).

First flower bloom day of Canada
The first flower bloom (FFB) day of Canada is obtained from phenology records of PlantWatch Canada Citizen Science network. FFB is defined as a plant stage at which the first flower buds have opened in an observed tree or shrub or in a patch of smaller plants. We have selected only the FFB day records observed by at least five observers at a minimum of five different locations in order to remove observer bias for 19 Canadian plant species recorded by several observers across Canada between years 2001 and 2012 totalling 1784 unique site-year data points (Gonsamo et al., 2013).

Sea ice extent and concentration
Annual means calculated from the daily sea ice extent (SIE) and concentration (SIC) observations are obtained from the Scanning Multichannel Microwave Radiometer (SMMR;1980-August 1987 and the Special Sensor Microwave/Imager (SSM/I; July 1987 to present) onboard the Nimbus-7 satellite and Defense Meteorological Satellite Program, respectively. The data are provided by the National Snow and Ice Data Center (NSIDC) (Fetterer et al., 2009). SIC is the fraction of Ocean area covered by sea ice whereas SIE is the total area covered by at least 15 percent of ice. The SIE and SIC data sets are for the Northern Hemisphere (> 31 • N) for the period of 1980-2012.

Snow cover
Annual means calculated from the monthly mean snow cover (SC) extent (Robinson et al., 2012) are obtained from the Rutgers University Global Snow Lab (available at http:// climate.rutgers.edu/snowcover). The SC extent is based on AVHRR satellite observations. The SC data set is for the entire Northern Hemisphere for the period of 1980-2012.

Sea level
Annual means calculated from 10-day estimates of sea level (SL) are obtained from University of Colorado Sea Level Research Group Research Group (available at http://sealevel. colorado.edu). The SL estimate was derived from the TO-Pography EXperiment (TOPEX) and Jason series of satellite radar altimeters calibrated against a network of tide gauges (Nerem et al., 2010). The SL data set is for entire Northern Hemisphere for the period of 1993-2012.

Surface temperature
The annual mean Northern Hemisphere surface temperature was obtained from the GISS (Hansen et al., 2010) data set (available at http://data.giss.nasa.gov/gistemp) for the period of 1980-2012. The global GISS and the 2013 reconstruction of Cowtan and Way (Cowtan and Robert, 2014) Had-CRUT4 hybrid UAH temperature anomalies are shown in Fig. 1b. The Cowtan and Way reconstruction of HadCRUT4 temperature data corrects for the incomplete global coverage, thereby alleviating the cool bias in recent decades (Cowtan and Robert, 2014). The global mean annual surface temperature anomalies, total solar irradiance, and global well-mixed greenhouse gases (WMGHG) patterns during the satellite era. The WMGHG includes CO 2 , major gases (CH 4 , N 2 O, CFC12 and CFC11) and a set of 15 minor longlived halogenated gases (CFC-113, CCl 4 , CH 3 CCl 3 , HCFCs 22, 141b and 142b, HFCs 134a, 152a, 23, 143a, and 125, SF 6 , and halons 1211, 1301 and 2402). Global annual mean surface temperatures relative to 1951-1980 base period from GISS (black) and Cowtan and Way HadCRUT4 hybrid UAH reconstruction (grey), and total solar irradiance are also shown in (b). (c) 11-year moving averages of global annual mean surface temperatures from GISS and Cowtan and Way HadCRUT4 hybrid UAH reconstruction, total solar irradiance, Atlantic Multidecadal Oscillation (AMO) index, and global stratospheric aerosol optical thickness at 550 nm anomalies relative to 1951-1980 base period. The zoom-in arrows in (c) show the time period covered in this study when the temperature anomalies diverged from solar forcing and internal climatic variability.

Atmospheric CO 2 measurements at Point Barrow station
Monthly averaged atmospheric CO 2 concentrations at Point Barrow station (Alaska, USA, 71.3 • N, 156.6 • W), based on continuous in situ observations, are obtained from the Earth System Research Laboratory (ESRL) of the National Oceanic and Atmospheric Administration (NOAA) (available at http://www.esrl.noaa.gov/gmd/dv/data). Atmospheric CO 2 concentration measurements from in situ stations cover the period of 1980-2012. The peak-to-though amplitude (AMP) for an individual year is calculated as a difference between maximum and minimum of monthly means to avoid influences of different curve fitting and data smoothing meth-ods. The annual means of parts per million (PPM) of atmospheric CO 2 concentration were also used in this study.

Sunspot measurements
Annual means calculated from the international daily mean sunspot (SP) number were obtained from the SIDC (Solar Influence Data Center) at World Data Center for the Sunspot Index, Royal Observatory of Belgium (available at http://www.sidc.be/sunspot-data). The sunspot number data used in this study cover the period of 1980(SIDCteam, 1980.

Solar irradiance
Annual means calculated from version d41_62_1302 daily averaged solar irradiance (W m −2 ) are obtained from Physikalisch-Meteorologisches Observatorium Davos World Radiation Centre (PMODWRC) (available at ftp://ftp. pmodwrc.ch/pub/data/irradiance). The composite algorithm, corrections for the radiometers other than VIRGO (Variability of solar IRradiance and Gravity Oscillations), and detailed methodologies are given in Frohlich and Lean (1998) and Frohlich (2000Frohlich ( , 2003Frohlich ( , 2006. The Active Cavity Radiometer Irradiance Monitor (ACRIM) II data are used to fill the gap during the Solar and Heliospheric Observatory (SOHO) Vacations in 1998 and early 1999. The solar irradiance (RAD) data used in this study cover the period of 1980-2012. The long-term solar irradiance data shown in Fig. 1c are obtained from (Krivova et al., 2007) for the  period and were merged with the PMODWRC data for 1980-2012.

Cosmic ray measurements at Kiel station, Germany
The hourly pressure corrected cosmic ray (CR) neutron monitor data of Kiel neutron monitoring station (Kiel, Germany, 54.3 • N, 10.2 • E) is obtained from the National Geophysical Data Center (NGDC) of the National Oceanic and Atmospheric Administration (NOAA) (available at ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/ COSMIC_RAYS/STATION_DATA/Kiel). The annual mean of hourly CR count calculated from hourly data for the period of 1980-2007 was used in this study.

Northern Hemisphere teleconnection indices
We restricted the teleconnection indices to those that dominate the interannual variability of climatic oscillations in phase and amplitude with continental to global scale implications accounting for the most spatial variance of the observed standardized anomaly (Quadrelli and Wallace, 2004;IPCC, 2007). The eight teleconnection indices (Barnston and Livezey, 1987;Wallace and Gutzler, 1981 (Gonsamo and Chen, 2015;Gonsamo et al., 2016). We then removed trends from the resulting winter teleconnection index by detrending the time series for the 1982-2011 base period.

Decadal teleconnection index
The Atlantic Multidecadal Oscillation (AMO) index is calculated from Kaplan sea surface temperature (SST) data set and is obtained from NOAA Earth System Research Laboratory (available at www.esrl.noaa.gov/psd/data/timeseries/AMO).

Stratospheric aerosol optical thickness
Annual means calculated from monthly mean stratospheric aerosol optical thickness at 550 nm are obtained from the GISS (Sato et al., 1993) data set (available at http://data.giss. nasa.gov/modelforce/strataer) for the period of 1850-2011. Background aerosols with an optical thickness 0.0001 were added as a lower limit for aerosol amount at all times. The effective radius (R) of the aerosol particles is defined as R = 0.20 (in the state with small optical thickness); other- is a function of time derived from the observed R for Pinatubo, while keeping the observed values for El Chichon and Pinatubo.

Analysis
For both response and driving variables, we have included observational data sets within the Northern Hemisphere (0-90 • N, 180 • E-180 • W) (Fig. 1a) available between 1980 and 2012. We present interannual variability analysis from detrended data to examine correlations at various timescales and minimize the risk of detecting spurious correlations. All interannual variability assessments were done based on detrended time series at annual timescale using the common base period of each pair of analysis. Trend analysis based on the raw data are presented only in Figs. 2 and 3, to show the long-term trends in both response and driving variables. Our interpretation of the results starts with basic correlation analysis (Table 2) comprising percent explained variance (coefficient of determination) among all variables from both raw and detrended data sets for trend and interannual covariability analyses, respectively. To investigate if the correlations are worthy of interpretation and to illustrate the coherence by analysing which variables are similar (or different), we use the squared loading of the variables (Abdi and Williams, 2010) from orthogonally projected driving and response variables, using principal component analysis (PCA). The squared loading of variables used in this study is alternatively called "squared cosines". We used a PCA algorithm with Pearson correlation coefficient as index of similarity to remove the effect of scale. Alternatively, this is called normalized PCA. The squared loading of variables can be interpreted numerically as the coefficient of determination between a PCA axis and a given variable, and reveals the internal structure and auto-association among the response and driving variables in a way that best explains the total variance. We use the 95 % confidence level from a two-tailed Student's t test to identify variables contributing significantly to each PCA axis. Generally speaking, the squared loading of the variables help interpret which variables are significantly coherent from the point of view of total variance analysis. The squared loading of variables that are small and not significant are interpreted as likely an artifact of projection into a low dimensional subspace, or an indication that the observed changes in response and driving variables are not coherent. Furthermore, we use several sets of PCAs (including natural and anthropogenic drivers together and separately with or without temperature mediation) to show the relative contribution of natural variability and anthropogenic forcing, and to predict each of the indicator variables with and without temperature mediation. We follow stepwise regression with Akaike Information Criterion (AIC) using different sets of PCA coordinates as regressors to reduce the effects of multicollinearity. All trend slopes in this study are calculated using a simple least squares linear regression.  The numbers indicate the statistically significant (p <

Results and discussion
This section starts with the Northern Hemisphere temperature trend analysis followed by the results and detailed interpretations of four groups of variables, i.e., spring phenology indicators, snow cover, sea level, and finally the atmospheric CO 2 dynamics in response to sea ice decline and climate variability.

Trends in the Northern Hemisphere surface temperature
The Northern Hemisphere experienced increases in surface temperature during the last 3 decades that are unprecedented in the anthropocene era ( Fig. 1a) with climate extremes during the period 1990-2010, which included the warmest decades since the start of modern measurements around 1850 (Fig. 1c). Warming in these recent decades is larger over land than over ocean (Fig. 1a), in part because the ocean responds more slowly than the land due to the ocean's large thermal inertia (Hansen et al., 2010). Warming during the past 3 decades is enhanced in Eurasia and the Arctic (Fig. 1a).
Warming of the ocean surface has been largest over the Arctic Ocean, and smallest and even slightly cooler over the North Pacific Ocean (Fig. 1a), partly due to the La-Niñalike cooling in the tropics affecting the extratropics (Kosaka and Xie, 2013). During the study period, the radiative effects from the increased WMGHG concentrations follow the rise in global surface temperatures (all p < 0.1 × 10 −7 ), whereas the solar irradiance is not and has an overall declining trend (Fig. 1b). The timeframe covered in this study coincides with the period when the global temperature anomalies diverged from trends in solar forcing and the internal climatic oscillation indicated by AMO (Fig. 1c).

Spring phenology of vegetation and soil thaw
We quantitatively assessed the spring anomalies of the start of growing season (SOS), spring thaw (ST), and first flower bloom (FFB) days. FFB, unlike hemispheric averaged spring indices such as SOS and ST, shows the contrasting roles of NAO and SCA on North American and Eurasian part of the Northern Hemisphere (see Fig. 4). Ob- Long-term trends of ST and SOS are overwhelmingly correlated with changes in annual mean surface temperature (p < 0.001) (Fig. 2b). ST and SOS are also significantly correlated with temperature after data detrending (Table 2) indicating both long-term and interannual covariability (p < 0.01). Changes in Canadian FFB are moderately explained by changes in annual mean surface temperature of the Northern Hemisphere (Fig. 2b), although there was, unsurprisingly, a stronger association with Canada's annual mean temperature (R = −0.85, p < 0.001) (Gonsamo et al., 2013). Associated with post-1998 slow down in surface temperature increase, SOS and ST advanced more slowly, even with   ( slight delays (Fig. 2b). Although the Canadian FFB shows strong correlation with NAO and SCA (Table 2), the stepwise regression selects 9 PCAs for FFB prediction (Fig. 5g, o) and the degrees of freedom becomes zero -indicating that there is no way to affirm or reject the prediction model for FFB. The interannual changes in solar radiation, NAO, and SCA are the most covarying variables with spring phenology anomalies (Tables 2 and 3). The interannual variability of spring phenology indices are explained more by natural forcing (i.e., solar radiation) and teleconnections than greenhouse gases (GHG) ( Table 4). Temperature mediation on interannual variability of spring phenology indices is only apparent with GHG, and less relevant with natural forcing (Table 4). Solar radiation and teleconnections may have non-temperature mediated effects on spring phenology through their impacts on incident solar radiation, cloudiness, precipitation and snowfall. The timing of spring events in many plant life cycles is advancing in response to climate warming (Parmesan and Yohe, 2003;Parmesan, 2006;Walther et al., 2002;Menzel et al., 2006;Post et al., 2009;Fitter and Fitter, 2002;Barichivich et al., 2013). The observed earlier spring activities (Fig. 2b) of terrestrial ecosystem increase the length of the growing season and consequently the primary productivity of vegetation. This same condition may also increase soil and plant respiration (Piao et al., 2008) and expose plants to widespread late spring frost damage (Hufkens et al., 2012), leading to carbon loss. However, the tradeoffs between increased primary productivity and enhanced ecosystem respiration and soil carbon release related to advancing spring activity remain poorly understood.
NAO and SCA are two of the most dominant teleconnections related to dynamics in terrestrial ecosystems of the Northern Hemisphere ( Fig. 4 and Tables 2 and 3) (Gonsamo and Chen, 2015;Gonsamo et al., 2016). NAO has a strong negative relation with SCA (Table 2), affecting much of Canada and Eurasia, with SCA dominant in Midwestern Europe (Fig. 4). Therefore, we only discuss the NAO results in relation to spring activity. The detrended NAO index is negatively correlated with fluctuations in snow cover   Table presents the percent explained variance derived from a PCA and stepwise regression analysis of land, cryosphere and ocean indicators by interannual changes in greenhouse gases and temperature (GHG-T ); solar radiation and eight leading teleconnections (SR-TEL); and solar radiation, eight leading teleconnections and temperature (SR-TEL-T ). Temperature in GHT-T and in SR-TEL-T are included to show the temperature mediated natural variability and anthropogenic forcing. Snow cover = SC, sea ice extent = SIE, sea ice concentration = SIC, start of growing season = SOS, peak-to-trough amplitude of CO 2 = AMP, spring thaw = ST, first flower bloom day = FFB, and sea level = SL. GHG includes concentration (PPM) of atmospheric CO 2 at Point Barrow and 20 global well-mixed greenhouse gases (WMGHG) radiative forcing. SR-TEL includes solar radiation, sunspot numbers and eight leading Northern Hemisphere teleconnection indices. All three sets (columns) of PCA analyses were done separately each only for the selected predicting variables. a p < 0.05, b p < 0.01, c p < 0.001, two tailed Student's t test. The explained variances by combination of all drivers with and without temperature mediation are given in Fig. 5. (p < 0.01), and positively correlated with changes in the FFB days (p < 0.01) (see lower left in Table 2). A steeper atmospheric pressure gradient (the high or positive NAO index phase), indicating an intensified Icelandic low, is associated with warmer Northern Hemisphere (mostly Europe and Asia) winter temperatures. This explains the negative relationship between the detrended NAO index and snow cover observed in our analysis (Table 2). Under steeper atmospheric pressure gradient or positive NAO (i.e., negative SCA index) phase when the westerlies in the North Atlantic are shifted poleward, there is enhanced advection of warm air across North-ern Europe and Asia, increasing vegetation productivity on this region (Gonsamo and Chen, 2015) (Fig. 4). Continental winter temperatures to the east are raised as a consequence.
To the west in northern Canada and Greenland the winters are colder and drier, delaying the Canadian first flower bloom days (Table 2) and overall vegetation productivity (Fig. 4).

Snow cover
Although the response of snow cover (SC) to global warming is complicated, as snow formation and melt are closely related to a temperature threshold of 0 • C (Brown and Mote, 2009), SC is the most predictable indicator (Fig. 5a, i) among the studied variables while teleconnections and solar radiation alone explain more than 74 % of the interannual variability (Table 4). The SC decline over the Northern Hemisphere of 0.14 × 10 6 km 2 decade −1 for 1980-2012 was not statistically significant for the region (p = 0.26), and showed less interannual variability after the 1998 global warming slowdown (Fig. 2c). Figure 5a, i and Table 4 show that SC is well explained by teleconnections and solar radiation whereas temperature mediation has only a marginal effect. The two leading Northern Hemisphere teleconnections (i.e., NAO and SCA) contribute the biggest natural climatic contributions to the interannual SC variability (see Table 2 and PCA1 column in Table 3). Temperature mediation on interannual variability of SC is only conspicuous with GHG, and less relevant with internal climatic variability (Table 4).

Sea level
For timescales relevant to anthropogenic warming, the rate of sea level (SL) rise is roughly proportional to the magnitude of warming above the temperatures of the pre-industrial age, with a proportionality constant of 3.4 mm year −1 • C −1 Each PCA for the two sets (i.e. left and right panels) were conducted separately. (Rahmstorf, 2007). Our analysis of simple linear trend shows that over the past 20 years, Northern Hemisphere SL has risen at a rate of 27 mm decade −1 for 1993-2012. Natural factors (solar radiation and teleconnections) impacting temperature explains 63 % of SL interannual variability (Table 4), while the combination and interactions of all stud-ied driving variables together with temperature explain 78 % of the observed variability (Fig. 5h). Although long-term SL rise is related to temperature rise (Fig. 2d) (Rahmstorf, 2007), the interannual variability is mostly controlled by temperature mediated (Fig. 5h, p and Table 4) changes in PNA and WR teleconnections (Table 2). PNA and WR modulate the location and strength of jet streams and fluxes of heat, moisture and momentum and can thus directly warm and expand, or cool and contract large areas of Northern Hemisphere water. Locally, the possible link between SL and teleconnections could be through changes in the surface atmospheric pressure via the inverse barometer effect, and water balance and density changes in response to temperature. Both PNA and WR are highly related to NINO (Table 2), and PNA phases are related to warm and cold Pacific episodes and sea level (Bromirski et al., 2011).
3.5 Atmospheric CO 2 variation in response to sea ice and climate variability The concentration of Northern Hemisphere atmospheric CO 2 decreases in spring as vegetation grows, and increases in fall when vegetation senesces resulting in an annual peakto-trough amplitude (AMP) of CO 2 concentration. The seasonal cycles of the Point Barrow CO 2 concentration is mainly explained by dynamics of growing and shrinking extratropical land ecosystems (e.g., Graven et al., 2013;Barichivich et al., 2013). The monthly Point Barrow measurements show that the CO 2 AMP has increased over the last 3 decades at a rate of 0.96 ppm decade −1 (p = 1.2 × 10 −6 ) for 1980-2012. Both CO 2 AMP and concentration (PPM) increases are significantly (p < 0.001) correlated with the long-term temperature increases (Fig. 2e) but changes in temperature do not directly explain the interannual variability (see lower left Table 2). The interannual variability of CO 2 AMP is explained by large-scale teleconnections such as EA, PNA and NINO and their temperature mediation (Tables 3-4), although the direct explanatory power of temperature on CO 2 AMP is negligible (Fig. 5e, m and Table 2). Our results show that there is no direct interannual link between CO 2 AMP and PPM in the Northern Hemisphere -the former is controlled by EA, PNA and NINO and their temperature mediation while the later is controlled by the influence of POL on sea ice dynamics (Table 3). Warm ENSO phases (i.e. positive NINO), coincides with lower CO 2 AMP (Table 2) indicating decreased CO 2 sink capacity which is in agreement with previous finding (Miralles et al., 2013). This decrease of CO 2 sink during positive NINO phase is due to reduced CO 2 uptakes by northern biosphere and may not be linked to sea ice dynamics (Table 3). In other words, the interannual variability of seasonal dynamics of CO 2 concentration is mostly controlled by EA, PNA and NINO influence on temperature while the absolute interannual variability in PPM is controlled by the POL influence on sea ice dynamics (see details below). This is assuming the trend in anthropogenic CO 2 emission remains constant, or the interannual changes in the trends are negligible. Following a decade with nine of the lowest minima on record, sea ice concentration (SIC) and extent (SIE) have received increased attention in light of climate warming (Post et al., 2013). Over the last 3 decades there have been rapid declines in both SIE (0.53 × 10 6 km 2 decade −1 ) and SIC (1.8 % decade −1 ) (both p < 0.5×10 −11 ) for 1980-2012. The decline rate in SIE is much faster than that of the SIC (Fig. 2a), indicating that sea ice is diminishing more rapidly in areas with thinner ice cover. The rapid decline of SIE is highly correlated with temperature rise (R = −0.8, p < 0.2 × 10 −7 ) (Fig. 2a). Our results show that the interannual variability of SIE and SIC are less controlled by temperature (Table 4), the least predictable indicators (Fig. 5b, c, j, k), more affected by POL teleconnection, and have the biggest direct control on atmospheric CO 2 concentration at Point Barrow but not globally ( Table 3). The interannual changes in CO 2 concentration are negatively related to changes in sea ice extent (p < 0.01) and concentration (p < 0.001) (see lower left Table 2).
The rapid changes in the Arctic are a consequence of the enhanced warming that the Arctic experiences compared with the rest of the world both on land and in the ocean, caused by a complex interaction of forcing and feedbacks, known as Arctic amplification. Inferring causality between correlated time series is difficult but may be supported when the sea ice response and feedback displays the expected physical understanding. There could be several explanations for the negative relationship between sea ice extent and atmospheric CO 2 concentration. First, water column stratification due to ocean freshening from melting sea ice restrain nutrient availability in Arctic primary productivity (Wassmann et al., 2011). Second, sea ice decline may indirectly contribute to periodic massive pulses of terrestrial carbon release as shown by the link between ice loss and the annual extent of tundra fires in Alaska (Post et al., 2013;Hu et al., 2010). Third, sea ice algae and sub-ice phytoplankton account for more than half of the total annual primary production in the Arctic Ocean (Gosselin et al., 1997), thus the decline in sea ice contributes to substantial loss of habitat for the primary producers. Forth, parallel to changes in the oceanic cryosphere, the lengthening of the growing season and a reduction in snow cover have also been observed in terrestrial ecosystems across the Arctic, which may induce large releases of carbon due to permafrost thaw (Schuur et al., 2011). On the other hand, sea ice decline can also contribute to increased carbon uptake. First, large phytoplankton blooms in the Arctic, where light transmission has increased in recent decades due to the thinning ice cover and proliferation of melt ponds can increase carbon uptake (Arrigo et al., 2012). Several studies show that rapid decline in sea ice related to climate warming is responsible for the increased sub-ice primary production (e.g., Post et al., 2013;Parmentier et al., 2013). Second, solubility-driven sea carbon uptake increases with increased ice-free sea surface. Third, sea ice decline is also strongly linked to longer growing season and increased vegetation productivity of the circumpolar north (> 45 • N) terrestrial ecosystems (Gonsamo and Chen, 2016) which indicates enhanced carbon uptake by northern plants. Generally speaking, through temperature and sea ice dynamics (Fig. 2a), the ocean may have a large impact on the terrestrial greenhouse gas balance of the Northern Hemisphere: earlier snowmelt and higher land surface temperatures -leading to longer growing seasons -can potentially increase plant uptake of atmospheric CO 2 , and these same conditions also increase respiration, permafrost thaw, wildfire, and droughts. Overall, our analysis strongly suggests that the increased carbon loss due to sea ice declinerelated processes may outweigh the carbon uptake enhancement through the parallel and concomitant processes, at least during the current climate regime.
The Polar/Eurasia Pattern (POL), which enhances the strength of the circumpolar vortex during its positive phase, is related to gradients in total mass of the atmosphere between polar and continental regions. The ice-albedo feedback due to declining sea ice results in warmer Arctic sea surface temperature, which increases ocean heat content and evaporation in polar region, further decreasing the temperature gradient of polar and continental regions. This may, in turn, result in strong negative POL phase events that lead to a weaker circumpolar vortex, and the resulting cold air spill will delay the spring vegetation activity of continental areas, reducing the CO 2 sequestration by terrestrial ecosystem. The negative phase of POL is strongly linked to decreased vegetation productivity of the circumpolar north (> 45 • N) terrestrial ecosystems (Gonsamo and Chen, 2016). With disproportionately accelerating warming of the polar region, the negative phase of POL will be prevalent resulting in less sea ice extent, and colder winters in continental areas. Currently, our results suggest that the sea ice loss is linked to net increase in atmospheric CO 2 concentration (Table 2). Given the above explanations, the anticipated sea ice decline in the future may lead to increased atmospheric CO 2 concentration, further strengthening the vicious circle of Arctic amplification.

Summary and concluding remarks
Prior to data detrending, our results reveal strong long-term relationships between temperature and several land, ocean and cryosphere indicators (Figs. 2 and 3). From Fig. 1, it seems that the atmospheric CO 2 forcing time series has less interannual variability but shows strong long-term relationships with rising temperature. When both the response and driving variables are detrended, the relationship between the long-term trend of temperature and land, ocean and cryosphere, and CO 2 forcing on temperature is partially removed. Consequently, the effects of the rapidly adjusting interannual variability of solar output and teleconnections become evident on several indicator variables. Unlike a single climatic variable such as temperature or precipitation, teleconnections control the entirety of heat, moisture and momentum fluxes, and incidence radiation through their effects on cloudiness (IPCC, 2007). This makes solar output and teleconnections the main drivers of the interannual variability of land, cryosphere and ocean indicators. However, new evidence is emerging regarding external forcing precursors on teleconnections (Fowler et al., 2012;Risbey et al., 2014;Collins, 2005), which may intensify in the midst of long-term climate changes. There is no relationship between solar irradiance and sunspot numbers with key land, cryosphere and ocean indicators (Fig. 3) if the trends are not removed from the data sets. This suggests that the recent trend in solar output has no discernible influence on the trends of the physical and biological systems indicators studied in the current work.
We found several coherent interannual patterns among related detrended response and driving variables. There are three sets of statistically strong auto-associations of driving and response variables which have at least 30-years of observations in the Northern Hemisphere (Table 3): (i) start of season, snow cover, Sun outputs, global well-mixed greenhouse gases (WMGHG), North Atlantic Oscillation (NAO), and Scandinavia Pattern (SCA); (ii) temperature, peak-to-trough CO 2 amplitude, East Atlantic Pattern (EA), Pacific/North American Pattern (PNA), and ENSO; and (iii) sea ice extent and concentration; CO 2 concentration (PPM), and Polar/Eurasia Pattern (POL). Overall, our results show that key land, cryosphere and ocean indicators are behaving as expected if they are responding to rising annual mean surface temperature and atmospheric CO 2 concentration in the Northern Hemisphere, and global well-mixed greenhouse gases (WMGHG), over the last 3 decades (Figs. 1b and 2).
The long-term trend analysis indicates that changes in surface temperature in the last 3 decades are strongly correlated (p < 0.05) with sea ice and sea level, spring phenology and thaw, and atmospheric CO 2 concentration (Fig. 2). Globally, rising temperature is also on a par with increasing radiative forcing of WMGHG (Fig. 1b). Recent changes in the Sun's output, decadal climatic oscillations, sunspot number, and cosmic ray counts have little or no relationship with long-term trends of Northern Hemisphere warming and its effect on land, cryosphere or ocean indicators (Figs. 1 and 2). During the last 3 decades, the Sun's energy output followed its historical 11-year cycle, with a slight overall decrease (Fig. 1b), temperature anomalies diverged from solar forcing, stratospheric Aerosol, and internal climatic oscillation indicated by the Atlantic Multidecadal Oscillation (Fig. 1c), and the two major volcanic eruptions of the last 3 decades have had only short cooling effects on climate (Gillett et al., 2012). Therefore, the combination of solar and volcanic activity should actually have led to a slight cooling if they were the primary drivers of long-term trends (Gillett et al., 2012). The recent multidecadal warming of Northern Hemisphere surface temperature cannot be explained by natural variability, or by any known mode of internal variability (Santer et al., 2013a, b) (Fig. 1c). Slow changes in the Earth's tilt and orbit around the Sun are only relevant in timescales of several thousands of years and cannot explain the recent rapid warming. Therefore, the observed rapid climate warming and its impacts on land, cryosphere, and ocean may best be attributed to anthropogenic factors, largely the radiative effects from increased concentrations of WMGHG (Fig. 1b). Despite the apparently slower rate of post-1998 global warming, a coherent pattern of changes across multiple life supporting natural systems is very likely to continue with increasing greenhouse gases.
How robust are our results? Although most of the variables were spatially averaged and multicollinearity was removed, the uncertainties from residual atmospheric effects and calibration errors in satellite data, missing drivers, errors in ground measurements, and methodological difficulty in capturing interactive effects of drivers and short-term feedbacks, are data source specific, difficult to quantify and cannot be ruled out. This work; however, contributes not only to observation-based detection and attribution of changes in climate index (i.e., here temperature), but also to the detection and attribution of impacts of climate changes on physical and biological systems following the 2014 Working Group II IPCC report (IPCC, 2014a) and other recent works (e.g., Rosenzweig et al., 2008;Parmesan et al., 2013;Stone et al., 2013;Poloczanska et al., 2013).