Analysis of CO2 mole fraction data: first evidence of large-scale changes in CO2 uptake at high northern latitudes

Atmospheric variations of carbon dioxide (CO 2 ) mole fraction reﬂect changes in atmospheric transport and regional patterns of surface emission and uptake. We report new estimates for changes in the phase and amplitude of observed high northern latitude CO 2 seasonal variations, indicative of biospheric changes, by spectrally decomposing 5 multi-decadal records of surface CO 2 mole fraction using a wavelet transform to isolate the changes in the observed seasonal cycle. We also perform similar analysis of the ﬁrst time derivative of CO 2 mole fraction, ∆ t CO 2 , that is a crude proxy for changes in CO 2 ﬂux. Using numerical experiments, we quantify the aliasing error associated with independently identifying trends in phase and peak uptake and release to be 10–25 %, 10 with the smallest biases in phase associated with the analysis of ∆ t CO 2 . We report our analysis from Barrow, Alaska (BRW) during 1973–2013, which is representative of the broader Arctic region. We determine an amplitude trend of 0.09 ± 0.02 ppmyr − 1 , which is consistent with previous work. Using ∆ t CO 2 we determine estimates for the timing of the onset of net uptake and release of CO 2 of − 0.14 ± 0.14 and − 0.25 ± 0.08 daysyr − 1 , 15 respectively, and a corresponding uptake period of − 0.11 ± 0.16 daysyr − 1 , which are signiﬁcantly di ﬀ erent to previously reported estimates. We ﬁnd that the wavelet transform method has signiﬁcant skill in characterizing changes in the peak uptake and release. We ﬁnd a trend of 0.65 ± 0.34 % ( p < 0.01) and 0.42 ± 0.34 % ( p < 0.05) for rates of peak uptake and release, respectively. Our analysis does not provide direct 20 evidence about the balance between uptake and release of carbon, but changes in the peak uptake and release together with an invariant growing period length provides indirect evidence that high northern latitude ecosystems are progressively taking up more carbon.


Introduction
Combustion of fossil fuel and cement production represent the dominant annual source of atmospheric CO 2 .There is also a minor source from the combustion of biomass and a diffuse source from the emissions and oxidation of reduced carbon (Suntharalingam et al., 2005).On an annual basis approximately 50 % of those emissions remain in the atmosphere with the remainder taken up by the land and ocean (Ballantyne et al., 2012).Regional changes to the net biospheric flux of CO 2 , and consequent changes in atmospheric CO 2 , are due to (a) spatial and temporal changes in climate, (b) different responses of vegetation to these changes in climate, and (c) other factors that may dominate over climate, e.g.nutrient availability.A recent study, building on extensive literature (e.g.Keeling et al., 1996), has reported substantial increases in the amplitude of the seasonal exchange of CO 2 since the 1950s, particularly at mid-to high northern latitudes (Graven et al., 2013).Here, we use the wavelet transform to isolate changes in the CO 2 seasonal cycle, revealing new insights about the growth rate, and changes in the amplitude and phase of CO 2 associated with the growing season.
Many previous studies have used atmospheric measurements of CO 2 to analyse changes in the seasonal cycle.
These studies have typically employed curve-fitting techniques (e.g.Bacastow et al., 1985;Thompson et al., 1986;Keeling et al., 1996;Piao et al., 2008;Barichivich et al., 2012Barichivich et al., , 2013) ) and spectral filtering methods such as complex demodulation (Thompson and Clark, 2008;Thompson, 2011).The most commonly used method is a combination of curve-fitting and spectral filtering (e.g.Thoning et al., 1989).We apply a wavelet transform, which uses a pre-defined wave-like oscillation that is noncontinuous in time or space to decompose a time series into time-frequency space, allowing us to investigate the dominant modes of variability and how they change with time.We use the wavelet transform as a filter to simultaneously decompose the CO 2 time series into seasonal, long-term, and residual components while retaining information about phase and amplitude in the seasonal cycle (Torrence and Compo, 1998).We show through extensive analysis of synthetic time series that using the time derivative of CO 2 , t CO 2 , provides more accurate estimates of changes in the phase of the CO 2 seasonal cycle and that the common use of the zero crossing points of the detrended CO 2 seasonal cycle is flawed.We show that we can also faithfully reproduce changes in the rates of peak uptake and peak release of CO 2 , allowing us to understand observed changes in the amplitude of the seasonal cycle.
In the next section we describe measurements of CO 2 mole fraction, the isotope ratio δ 13 C, surface temperature, and vegetation indices as well as the approach we have employed to impute these data.In Sect.3, we describe the wavelet transform that we use to spectrally decompose these data, including a characterization of the aliasing error associated with independent inference of changes in phase, amplitude, and the magnitude and timing of the peak uptake and release of CO 2 .In Sect.4, we present our analysis of CO 2 growth rates and changes in the phase and amplitude of the CO 2 seasonal cycle.We conclude our paper in Sect. 5.

CO 2 mole fraction data
Figure 1 shows the sites from the National Oceanic and Atmospheric Administration (NOAA) Global Greenhouse Gas Reference Network (GGGRN, Dlugokencky et al., 2015), which include at least 15 years of CO 2 mole fraction data.We focus on high northern latitude sites where (a) seasonal contributions of CO 2 are predominantly driven by boreal vegetation and (b) contributions to observed CO 2 from continents at these latitudes are approximately equal (Fig. 2).We report our CO 2 analysis for the site of Barrow, Alaska (BRW), because it is generally considered to be representative of the broader Arctic region and report our analysis from other sites in Appendix C.
Twin air samples are collected weekly at the sites and analysed for CO 2 at NOAA Earth System Research Labo- The NOAA/ESRL monitoring sites used in our CO 2 time series analysis.For the seasonal cycle analysis, we use the northern hemispheric high-latitude sites (blue).The sites shown in green, red, and magenta are used for growth rate analysis only.The six sites with a black border are those with the longest time series in each 30 • latitude band.The shaded regions are the temperate and boreal northern Hemispheric land regions defined in the initial Transcom study, and which we use for analysis of NDVI, temperature, and atmospheric transport.1) to the zonal mean concentrations over the high-, mid-, and low-latitude Northern Hemisphere averaged over 2004-2009.These values were determined by using the GEOS-Chem atmospheric transport model (see main text for further details).The error bars denote the 1σ of the year-to-year variability over the 6-year period.The zonal means are defined as the mean of the grid points sampled nearest to the sites shown in Fig. 1.Model output from December is missing.ratory (ESRL) in Boulder, Colorado, using a nondispersive infrared analyser.These data are suitable to study variations on weekly and longer timescales.Single measurement uncertainties are calculated based on the ability to propagate the World Meteorological Organisation (WMO) XCO 2 scale to working standards (±0.03 ppm; Zhao and Tans, 2006), the analytical repeatability of the analysers for a sample measurement (±0.03 ppm), and the agreement between pairs of samples collected within 20 s of one another (±0.1 ppm across the entire sampling network).The sample pairs are not collected simultaneously, so that the agreement contains an element of real atmospheric variability.The sum of these uncertainties is small in comparison to the magnitude of CO 2 variability observed at northern high latitudes.

Imputation of mole fraction data
The wavelet transform method (described below) requires a continuous time series that is regularly spaced in time.To fill a missing value in a time series we add a value from 7year average seasonal cycle (calculated using the 3 years on either side of and including the year of interest) to a value from a deseasonalised reference time series (Fig. 3), which accounts for large-scale anomalies in the growth rate in 30 • latitude bands.This ensures that gradual changes in the seasonal cycle amplitude/phase are preserved.Any remaining missing data points are extracted from a piecewise cubic spline curve fit.Parts of the time series that contain significant sections of missing data are likely to be unrepresentative of real changes; however, prolonged periods are rare and we find that isolated missing data points do not significantly impact the determination of long-term trends in the phase and amplitude.
Figure 4 shows an example of our imputation approach using the CO 2 mole fraction and δ 13 C time series from Cold Bay, Alaska.

δ 13 CO 2 data
We also use measurements of δ 13 C that are colocated with the CO 2 mole fraction data.The isotope samples are analysed at the Stable Isotope Laboratory at The Institute of Arctic and Alpine Research (White and Vaughn, 2011) using flasks of air provided by the NOAA GGGRN.These help us to attribute observed changes in CO 2 mole fraction to land biospheric uptake and release.The ratio δ 13 C is defined as where    vidual measurements of 12 C and 13 C are determined by isolating the CO 2 in a subsample of air from each flask and using a mass spectrometer to determine the isotopic composition.

Ancillary data
We use the University of East Anglia Climate Research Unit TS3.10 land temperature data set (Harris et al., 2013) to help interpret observed variations in the CO 2 time series.These data have a 0.5 • × 0.5 • spatial resolution and monthly time resolution.
To investigate large-scale vegetation change, we use the Global Inventory Modeling and Mapping Studies normalized difference vegetation index (GIMMS NDVI3g) data set derived from the NOAA Advanced Very High Resolution Radiometer (AVHRR) (Pinzon et al., 2005;Tucker et al., www.atmos-chem-phys.net/15/13739/2015/2005).NDVI, calculated from the visible and near-infrared light reflected by vegetation, is strongly correlated with photosynthetic activity in vegetation canopies; although we acknowledge photosynthesis may not accompany greenness (a) at high latitudes when water is frozen and (b) during drought when stomata are mostly closed.These NDVI data have a spatial resolution of approximately 8 km and a twice monthly temporal resolution from 1982 to the end of 2006.The data set has been corrected for calibration, viewing geometry, volcanic aerosols, and other effects not related to vegetation change.We remove pixels that have a time series mean NDVI value of < 0.1 to ensure that areas with bare or sparse vegetation are not included in spatial averages.

Wavelet transform
We use a wavelet transform to spectrally decompose the observed CO 2 variations into individual frequency bands that can be attributed to the responsible biological and physical processes.
In general a wavelet transform W n uses a wavelet function ψ 0 , a pre-defined wave-like oscillation that is noncontinuous in time or space, to decompose a time series into time-frequency space, allowing us to investigate the dominant modes of variability and how they change with time.This improves on the Fourier transform that determines frequency information using sine and cosine functions.
The wavelet transform of a time series x n is defined as where xk is the discrete Fourier transform of x n , N is the number of points in the time series, k = 0. ..N − 1 is the frequency index, and ψ • (sω k ) is the complex conjugate of the Fourier transform of a normalized, scaled, and translated version of ψ 0 (η), where s is the scale and ω k is the angular frequency.We use the Morlet wavelet, a plane wave modulated by a Gaussian envelope: where ω 0 is the nondimensional frequency and η is the nondimensional time parameter.We chose the Morlet wavelet because it is nonorthogonal, which is an attractive property for the analysis of smooth and continuous variations such as those exhibited by CO 2 mole fraction time series.The wavelet is comprised of a real and imaginary part, providing information about amplitude and phase respectively.We can recover the original time series from wavelet space using the corresponding inverse transform (Torrence and Compo, 1998) and summing over all frequencies from the real part of the wavelet transform or a subset of frequencies if we are interested in isolating signals: where ψ 0 (0) removes the energy scaling and s 1 2 j converts the wavelet transform to an energy density.C δ and ψ 0 (0) are constants determined for the specific wavelet function.
To minimize edge effects associated with the Fourier transform, we add synthetic data to pad the beginning and end of the time series.For our calculation we repeat the first (last) 3 years of the time series backwards (forwards) in time, accounting for continuity of the growth rate based on the following (preceding) years.The synthetic data used in the padding should be close to what we expect, but is essentially unknown, and this uncertainty penetrates the first and last year of the time series.We also "zero pad" the time series so that the number of points used is an integral power of 2, which further reduces edge effects and speeds up the transform.The padded data at the edges of the time series are removed post-wavelet decomposition and prior to analysis.
We quantify the numerical error associated with the wavelet transform by applying it to synthetic time series, which are representative of CO 2 time series with a prescribed trend.We find that the value for C δ previously reported (Torrence and Compo, 1998) introduces a small trend in the original minus reconstructed residual and find that C δ = 0.7785 results in a much smaller, unbiased residual with a typical value < 0.05 ppm for monthly data and < 0.002 ppm for weekly data (not shown).Table 1 shows the wavelet parameter values that we used in our analysis.
Additional uncertainties may arise in the long-term trend and detrended seasonal cycle as a result of spectral power being assigned to the incorrect frequency band.This could, for example, result in concentration changes caused by anthropogenic emissions being misattributed to the natural (seasonal) cycle of CO 2 and vice versa.However, this is a common weakness of any method used to decompose such time series.
We find that for atmospheric CO 2 , the wavelet power spectrum peaks at periods (reciprocal of frequency) of 6 and 12 months (Appendix A), with a spread across these periods associated with the sampling of the data.To study annual changes in phase and amplitude, we retain a period of Values for MLO, which are typically taken to be representative of the global growth rate, are highlighted with a circle.Time series of annual growth rates were determined for individual CO 2 measurement sites before calculating decadal mean growth rates and then binning the decadal mean growth rates into 20 • latitude bins.We subtract a global mean growth rate due to fossil fuel combustion from all sites.
3 to 18 months, and we assume that periods longer than 18 months are indicative of the growth rate and periods shorter than 3 months are due to local/regional sources that are unrelated to the seasonal cycle (described using an example in Appendix A).

Growth rates
Figure 5 shows how the decadal atmospheric growth rate has changed from the 1980s to the 2000s as a function of latitude.We find that in the 1980s and 1990s the growth rates are approximately the same in the Southern Hemisphere but diverge further north.The 1980-1989 growth rate rises sharply towards the northern high latitudes while there is a dip in the 1990-1999 in the same latitude band.We anticipate that this is primarily due to changes in biospheric uptake in the Northern Hemisphere.It should be noted that the number of CO 2 monitoring sites in the 1980s is considerably more sparse relative to later years, but this should not matter too much given the decadal averaging and the fact that the sites have been selected to be representative of background conditions.The 2000-2009 decadal mean growth rate is significantly higher than both of the previous decades by ∼ 0.35 ppm yr −1 and rises from the Southern Hemisphere to mid-latitude Northern Hemisphere before dropping off again in the northern high latitudes.We find that our annual CO 2 growth rates at Mauna Loa are within a fraction of a percent of NOAA values.
By subtracting anthropogenic fossil fuel emission estimates from the atmospheric CO 2 signal (Table 2) we can effectively isolate uptake by the oceans and terrestrial biosphere, acknowledging the uncertainties associated with the emission estimates and that we have not accounted for land use change emissions.The residual growth rate is negative, as expected (Ballantyne et al., 2012).We find that during the 1980s the net annual uptake by the terrestrial biosphere and ocean was typically −1.03 ± 0.11 ppm yr −1 when averaged across all sites and where the uncertainty is equal to 1σ .This rate increases dramatically in the 1990s to approximately −1.54 ± 0.06 and to −1.89 ± 0.08 ppm yr −1 in 2000s.This change in the growth rate supports the notion that the natural component of the carbon cycle is increasing the amount of carbon it takes up in response to the amount of carbon present in the atmosphere, although the last 2 decades show a smaller increase in net annual uptake.This apparent equilibrium state results in an approximate mean airborne fraction of 55.8 ± 18.2 % (including only fossil fuel) and 44.1 ± 14.4 % (including fossil fuel and land use change), consistent with previous work (Gloor et al., 2010).For the purpose of the following calculations we have removed the annual growth rate from the observed CO 2 concentrations, following the method described in Appendix A.

Phase and amplitude analysis
We use several metrics to interpret the CO 2 mole fraction time series but focus on phase changes and estimates of peak uptake (PU) and peak release (PR) from the first differential of CO 2 , t CO 2 , a proxy for the net flux of CO 2 .As we discuss below and report in Appendix B, analysis of t CO 2 leads to less biased estimates for trends in the phase of the CO 2 seasonal cycle.As part of our analysis we report 95 % confidence intervals, the Pearson correlation coefficient r, and p values that denote the probability of reproducing a result by chance; for practical purposes p values > 0.05 represent a result that is not significant.

Practical definitions and theoretical calculations
Figure 6 shows, using example data from BRW, how the detrended CO 2 and t CO 2 variations are related.The amplitude of the seasonal cycle, defined as the peak-to-peak difference (maxima minus minima) of the seasonal CO 2 mole fraction time series, has been used in previous studies as a measure of biological activity (e.g.Keeling et al., 1996;Graven et al., 2013).This metric alone cannot tell us whether net uptake or release is responsible for observed variations, so it is typically used as an indicator of overall carbon exchange.
Recent work has shown that the intense period of uptake during summer in the high northern latitudes contributes more to the seasonal amplitude than the longer period of emission in autumn.Based on t CO 2 we define three periods during an annual cycle: (1) an uptake period when t CO 2 < 0 and there is a net negative CO 2 flux to the atmosphere (photosynthesis is higher than respiration); (2) a release period when t CO 2 > 0 and there is a net source of CO 2 to the atmosphere; and (3) a dormant period, defined between the latter half of winter and the start of the next uptake period, when plant activity is very low due to frozen ground so that t CO 2 is typically small (but non-zero due to transport of CO 2 from the lower latitudes).
To look at changes in phase, previous studies have used the zero-crossing point (ZCP) of CO 2 which refers to times when the detrended seasonal cycle is equal to 0 (e.g.Piao et al., 2008).In one seasonal cycle there is a downward and upward ZCP (DZCP and UZCP respectively), where the DZCP is typically taken as a proxy for northern hemispheric spring onset of net carbon release and the UZCP is taken as a proxy for the onset of autumn net carbon release.The ZCPs of the CO 2 concentration can only be estimated from the detrended seasonal cycle.The long-term increase in CO 2 is driven by changes in net flux, and by detrending the seasonal cycle is shifted up or down relative to the zero line, such that the annually integrated flux is equal to 0. As such, an increase in net uptake in 1 year will cause a shift to the CO 2 DZCP and UZCP even if there is not a real change in phase.We refer to this error, associated with detrending the time series, as the The beginning of the period of net carbon uptake is difficult to determine accurately using the seasonal cycle at highlatitude sites because small mole fraction variations during the dormant period (which has a near-zero flux) are sufficient to bring t CO 2 below 0 before the carbon uptake period associated with the main growing season.To address this we tested a number of phase thresholds which represent the timing of when certain magnitudes of t CO 2 are reached (e.g. 25 % of PU).We find that using the 25 % of PU is a more robust indicator of the beginning of the net carbon uptake period and use this as our "spring" phase metric.In contrast, the t CO 2 UZCP is well defined and trivial to calculate and so we use this as our "autumn" phase metric.We define a carbon uptake period (CUP), which is the difference between the autumn and spring phase metrics defined above.PU and PR refer to the minima and maxima of the flux time series respectively.As we show below using theoretical calculations these peak values are related to annual release and uptake.
The ability to isolate changes in the phase and amplitude of the seasonal cycle with fidelity is critical for our analysis.We use Monte Carlo numerical experiments to characterize the errors associated with independently identifying changes in phase and amplitude that can result in the misinterpretation of these data and/or underestimation of uncertainties (Appendix B).These errors are not unique to using the wavelet transform but are more pronounced when using the detrended CO 2 seasonal cycle as opposed to t CO 2 .To our knowledge, no previous study has quantified these errors when estimating phase changes in the CO 2 seasonal cycle.We generally find that analysis of t CO 2 produces more reliable and less biased estimates than CO 2 trend estimation of either phase with an estimated 25 % systematic aliasing error (Appendix B).Unless explicitly stated all subsequent results will refer to our analysis t CO 2 .We also find that we can capture at least 80 % of independent trends in the PU and PR of the t CO 2 seasonal cycle, which has not been reported previously and allows us to study changes in characteristics more closely related to annual changes in biological release and uptake of CO 2 (Appendix B).

Analysis of NOAA CO 2 mole fraction data
Figure 7 shows that changes in spring and autumn phases determined from BRW t CO 2 are −0.14 days yr −1 (p < 0.05) and −0.25 days yr −1 (p < 0.01) respectively, with a corresponding CUP change of −0.11 days yr −1 (p > 0.1); the analysis of the other study sites is shown in Appendix D. We find no evidence using phase changes of CO 2 or t CO 2 for a significant change in CUP throughout the measurement period due to a simultaneous advance of the spring and autumn phase.
The concomitant observed changes in t CO 2 and in t δ 13 C (Fig. 9, also discussed in Appendix E) supports the idea that observed CO 2 variations are primarily due to changes in the terrestrial biosphere.Analysis of surface temperature reanalyses and space-borne observations of NDVI also corroborate the spring phase change of t CO 2 (Appendix E).We find the start of the thermal growing season (defined as the continuous period above 5 • C, Appendix E) is advancing 2 (3) times faster at latitudes > 45 • N (> 60 • N), which agrees with previous studies (e.g.Barichivich et al.,

2012
).However, we find an anti-correlation of autumn phase changes with NDVI and temperature anomalies.The NDVI anomalies during summer have not significantly increased on large spatial scales over the measurement period  compared with spring and autumn anomalies.This suggests that the increase in net exchange of carbon between vegetation and the atmosphere is likely a result of increased photosynthetic activity during spring and autumn.In contrast, our analysis of t CO 2 time series shows more uptake of CO 2 in spring and early summer and earlier onset of net release of CO 2 between mid-summer and autumn.A number of studies have linked increases in NDVI and subsequent carbon uptake with a CO 2 fertilization effect (Lim et al., 2004;Kaufmann et al., 2008;Los, 2013) which may be partly responsible for the observed increases in carbon uptake during this period.Our analysis of NDVI data shows that increases of vegetation greenness in spring and autumn have led to significant lengthening of the photosynthetic growing season over the measurement period, where autumn greening is changing in most regions at a greater rate than spring greening.The carbon uptake period, however, has not extended but shifted earlier in the year and retained its length.If photosynthesis has increased at the end of the growing season, and it is a change in the net ecosystem exchange that explains the change in phase, this implies that respiration must have increased more than photosynthesis to cause an advance of the phase at the end of the uptake period.Observed changes in amplitude at BRW (0.09 ± 0.02 ppm yr −1 ) are consistent in percentage terms with previous work over the same time period (Graven et al., 2013).We find that the observed change in amplitude at BRW is partially due to an increase in PR (0.42 ± 0.34 ppm yr −1 , p > 0.05) and a larger increase in PU (0.65 ± 0.34 ppm yr −1 , p < 0.01).Figure 8 shows that statistically significant trends (p < 0.05) in PU are observed at five of the seven highlatitude sites (ALT, BRW, CBA, ICE and ZEP, Table D1).In most of these cases, the change in PU is significantly larger than the change in PR.We show in Fig. B8 that trends in amplitude are determined mainly by changes in uptake during the CUP (Appendix D).Previous analysis of these data has shown that changes in atmospheric transport cannot explain changes in the amplitude (Graven et al., 2013).

Concluding remarks
We have used a wavelet transform to spectrally isolate changes in the seasonal cycle of atmospheric CO 2 mole fraction.The wavelet transform can simultaneously separate the long-term trend and seasonal cycle while retaining information about changes in amplitude and phase.We focused on high northern latitude sites where (a) seasonal contributions of CO 2 are predominantly driven by boreal vegetation and (b) contributions to observed CO 2 from continents at these latitudes are approximately equal.
We found that the atmospheric growth rate of CO 2 at these sites are within a few percent of reported values from NOAA.Our growth rates show large decadal changes, as expected, and once the anthropogenic signature has been removed we find strong evidence of a natural biospheric signal that is responding to increasing atmospheric CO 2 concentrations.This results in a near-constant airborne CO 2 fraction of 55.8 ± 18.2 % (including only fossil fuel) and 44.1 ± 14.4 % (including fossil fuel and land use change), consistent with previous studies.
Using the detrended CO 2 time series (original data minus growth rate) we examined the change in phase and amplitude of the seasonal cycle.Using a series of synthetic experiments we showed that using the first differential of CO 2 provided more accurate estimates of independent changes in phase and peak uptake and release of CO 2 , to within 10-25 % of the "true" values.
We reported an increase in amplitude of 0.09 ± 0.02 ppm yr −1 , consistent with previous studies, which can be crudely associated with an increase in biological activity.Using a series of Monte Carlo experiments we showed that amplitude changes are strongly correlated with trends in net carbon uptake during spring and summer but had a weak relationship with changes in net release of CO 2 in autumn and winter.We showed that in percentage terms, the rate of peak uptake has increased at a significant and faster rate when compared with the rate of peak release.
We diagnosed phase changes using thresholds of t CO 2 , taking the timing of uptake reaching 25 % of peak uptake as the beginning of the CUP, and the timing of t CO 2 switching to positive as the end of the CUP.These phase thresholds take into account that observed t CO 2 variations can introduce local maxima/minima particularly associated with the beginning of the CUP.We reported changes in the downward and upward phase of −0.14 ± 0.14 and −0.25 ± 0.08 days yr −1 respectively and a corresponding revision to the length of the net uptake period of −0.11 ± 0.16 days yr −1 .There is no evidence of a significant increase or decrease in the length of the CUP.Given that we characterized the method used to determine the change in phase, including a measure of uncertainty, and showed that analysing t CO 2 produced less biased estimates for these changes we argue that our values are a more faithful depiction of the truth.
Our analysis does not provide direct evidence about the balance between uptake and release of carbon, but changes in the peak uptake and release together with an invariant growing period length provides indirect evidence that high northern latitude ecosystems are progressively taking up more carbon in spring and early summer.The period of net carbon uptake has not lengthened but has become more intense.However, it is possible that this increase may be offset by a prolonged period of respiration due to warmer autumn temperatures.Changes in atmospheric CO 2 mole fraction tell us only part of the underlying carbon cycle story in terms of how the underlying ecosystems are changing.Clearly, additional measurements and models needs to be applied for us to understand observed changes in atmospheric CO 2 .A more frequent inspection of these data using advanced statistical tools such as the wavelet transform also has a role to play.

Appendix A: Example of spectral decomposition
Figure A1 shows, as an example, the spectral decomposition of CO 2 mole fraction measurements at Mauna Loa.The wavelet transforms decomposes the 1-D time series into a 2-D power spectrum, describing energy per unit time, as a function of frequency (the reciprocal of period) and time.The cone of influence is the boundary below which wavelet coefficients are most compromised by edge effects.We have padded the edges of the CO 2 time series with additional synthetic data so we are able to analyse the entire CO 2 time series (Sect.3).We find that most of the power is in the annual and semiannual periods, as expected, but also peaks in power at period > 1 year but this is likely a result of responses of the CO 2 growth rate to large-scale climate variability, e.g. the El Niño-Southern Oscillation (ENSO).This is supported by the global wavelet power spectra (integrated over all time).The interannual growth rate is determined by taking the value of the long-term trend (periods > 18 months) on 1 January in 1 year and subtracting the value from the previous year to leave the net change in concentration.
As discussed above, we use the spectrally decomposed data set to interpret the observed variability of CO 2 mole fraction data.Figure A1 shows two example applications: (1) as a lowpass filter to deseasonalize the CO 2 data (removing periods < 18 months) and ( 2) the associated annual growth rate (ppm yr −1 ), which we find is within < 0.1 ppm of the reported values from NOAA (not shown)., 1959, -2012. Middle row: (left) . Middle row: (left) the wavelet power spectrum where the colour scale is log (power).The black solid lines denotes the cone of influence.The power spectrum tends to emphasise very low frequency information so we have subtracted an exponential term prior to applying the wavelet transform to emphasise the high frequency variability (right) the corresponding time-integrated power spectrum.Bottom row: the inferred annual growth rate of CO 2 (ppm yr −1 ).

Appendix B: Error characterization of phase and amplitude estimates
We use synthetic CO 2 time series data, defined with specific changes in amplitude and phase, to characterize aliasing errors due to application of the wavelet transform of CO 2 concentration data or its first derivative ( t CO 2 ).Insights from this synthetic analysis are directly applied to our interpretation of NOAA mole fraction measurements in the main paper.

B1 Synthetic model framework
We use a simple box model based on the CO 2 mole fraction time series at Barrow, Alaska (Fig. B1).BRW is the most suitable site for this purpose because is has a long time series and as it is representative of high-latitude CO 2 in the Northern Hemisphere.We take the first derivative of the detrended time series at BRW to get the "flux" time series.We then take the mean seasonal cycle of the CO 2 flux and adjust it so that in its initial state, the source and sink terms are balanced.This cycle is then repeated for 40 years (equivalent to the time span of the BRW time series) and integrated to convert the flux to CO 2 concentration.For our experiments, described below, we introduce trends and variability to various aspects of t CO 2 before integrating with respect to time to recover CO 2 mole fraction.Detrending is as described in the main paper.

B2 Numerical experiments
The starting point of our numerical experiments is the detrended time series of atmospheric CO 2 mole fraction.Our analysis here as it is in the main paper does not provide direct evidence about the balance between uptake and release of carbon.The detrending process results in a seasonal cycle that integrates to 0 over a year, which can, if not properly accounted for, introduce false trends and variability in the seasonal cycle metrics.We combine the metrics defined above to provide indirect evidence of trends in the carbon balance of the Northern Hemisphere.The following three broad set of experiments are designed to identify the best metrics to describe changes in the contemporary cycle from detrended CO 2 mole fraction measurements.First, we perturb the timing of spring or autumn by adding or subtracting a smooth Gaussian curve with a flat top centred roughly about the onset of net uptake or release and increase the magnitude of the curve each year to introduce a trend across the time series.Second, we perturb the magnitude of net uptake or net release by multiplying the uptake (negative t CO 2 ) or release (positive t CO 2 ) by some factor, and increase the factor each year to introduce a trend.Finally, we add year-to-year variability (or noise) to the time series to assess the ability of our spectral method to extract trends from the data.We compare each metric by calculating the percentage difference in trend from the input time series and the wavelet detrended time series.

B2.1 Perturbing the timing of the spring and autumn phases
Figure B2 shows the results of our analysis of a time series for which we introduced a progressively earlier onset of net CO 2 uptake of 0.50 days yr −1 for t CO 2 DZCP.The t CO 2 DZCP is very sensitive to the curve we use to perturb the time series due to the relatively flat period of near-zero flux during the dormant period preceding it (it does not take much to bring this below 0).While for the synthetic example we have used a smoothed version of the BRW time series, in practice there is substantial variability in the spring shoulder so that it is often difficult to accurately define a trend in the t CO 2 DZCP.To address this we use an operational definition that is defined as 25 % from 0 to the PU, which in this case has a trend of 0.35 days yr −1 .The t CO 2 metrics were found to be better at capturing the springtime trend to within 23 and 16 % respectively when compared with the equivalent CO 2 mole fraction metric which underestimates the trend by 63 %.This has implications for using the CO 2 mole fraction ZCPs to interpret changes in the phase.There is little change in any of the UZCP metrics (typically < 0.025 days yr −1 ) as a result of aliasing.The wavelet detrending introduces a −0.01 % yr −1 trend in peak CO 2 uptake and a concurrent increase in peak CO 2 release of 0.14 % yr −1 corresponding to −0.4 and 5.6 % across the 40-year time series respectively.This is considered an aliasing error when interpreting the real data in the main paper and is relatively small considering the large trends introduced in spring uptake.
Figure B3 shows the same calculation but for introducing an earlier autumn onset of net CO 2 release of 0.30 days yr −1 .We find that the metrics for the spring phase respond to the prescribed change in autumn phase due to aliasing, where the mole fraction and t CO 2 = 0 metrics had non-zero trends up to ∼ −0.16 days yr −1 .All three UZCP phase metrics underestimate the change in the defined phase change by amounts ranging from 11 to 22 % where the CO 2 UZCP performed the  best.The earlier onset of net CO 2 release aliases into a 2.5 % increase in peak CO 2 release and a 5 % increase in peak CO 2 across the entire time series.

B2.2 Perturbing the magnitude of net uptake and release of CO 2
Figure B4 shows the results of introducing a progressive enhancement of CO 2 uptake of roughly 0.70 % yr −1 , equivalent to a 28 % increase over 40 years.We introduce the trend by multiplying the negative flux by an increasing amount each year, which does not have an effect on timing of net CO 2 uptake or release.We also introduce 2 exceptional years to emulate the effect of interannual variability such as that driven by climate phenomena like ENSO.We find that the wavelet transform attributes the 0.70 % yr −1 increased uptake as 0.59 % yr −1 uptake and 0.20 % yr −1 release.The mole fraction metrics infer non-zero UZCP and DZCP phase changes of 0.06 and 0.16 days yr −1 respectively, while the 25 % t CO 2 UZCP and DZCP metrics, our operational metrics, exhibit negligible trends as expected.The exceptional years are captured in the PU and PR metrics, while the CO 2 UZCP is the most affected out of the phase metrics.In addition, information from the exceptional years of uptake is aliased into the CO 2 UZCP and is spread over a number of years rather than just one.This is not the case for the t CO 2 metrics indicating that they are better for estimating interannual variability.

B2.3 Simultaneous variations in phase and peak uptake and release
Figure B5 shows the results from a final experiment that describes a calculation in which we simultaneously perturb the phase of the spring and autumn, as diagnosed by the t CO 2 = 0, and the PU and PR.We also superimpose Gaussian random noise within ±10 days and ±25 % to describe year-to-year changes to the phase and to the PU and PR respectively.
Despite large interannual variability, there is a negligible trend in the spring timing of CO 2 uptake (−0.02 days yr −1 ) which is captured by the t CO 2 phase metric (0.02 days yr −1 ).The CO 2 DZCP trend has the opposite sign and additionally overestimates the magnitude of the trend by a factor of 4. The trend in the autumn t CO 2 phase metric (0.05 days yr −1 ) underestimates the expected trend (0.09 days yr −1 ) by ∼ 45 %, while the CO 2 UZCP overestimates it by a factor of 2.8.The estimated trend in PU is 0.54 % yr −1 which is 80 % of the expected trend (0.68 % yr −1 ), while the estimated PR trend (0.14 % yr −1 ) is opposite in sign and double the magnitude of the expected trend (−0.07 % yr −1 ).The estimated CUP trend is positive but roughly 0, which is a little smaller than the expected trend of 0.12 days yr −1 .The increase in PU (which is a factor of 3 larger than the rise in PR) and the roughly 0 trend estimated for the CUP hints at a probable increase in annually integrated net uptake.The trend in net flux in this example is indeed negative with an increase in uptake of  We also superimpose Gaussian random noise to describe interannual variation.
We find that the analysis of synthetic time series indicates that t CO 2 metrics can reproduce prescribed phase changes to within 30 %, but trends with a magnitude of < 0.1 days yr −1 were uncertain in magnitude and sign.Strong shifts in spring and autumn phase caused changes in PU and PR of < 6 % due to aliasing.Strong trends in PU and PR were estimated to within 25 %.

B3 Monte Carlo simulations (MCSs)
We used an MCS to study the ability of the wavelet transform to simultaneously determine the PU, PR, and changes in phase.We generated 1000 synthetic time series with random trends and variability such as the one illustrated in Fig. B5, where Fig. B6 shows the probability distributions of the trends introduced in the net carbon fluxes and changes in the CUP.Trends in integrated uptake and release of carbon was in the range of −0.25 to 0.25 ppm yr −2 , while changes in the phase were within 1 day yr −1 .We then regressed the expected trends in phase, PU, and PR against the values we estimated using our analysis.The regression coefficient was used as an estimate of the mean bias, while the Pearson correlation coefficient r is indicative of consistency in the bias and the likelihood of the estimates to deviate far from the expected value.
Figure B7 shows some of the results from the MCS regression analysis where we compare expected and estimated trends.The figure also shows estimates where we detected the wrong sign of the trend and the quantity of statistically significant trends (p < 0.05) that were and were not detected in the analysis.The results of the MCS indicated not only a large mean negative bias in the CO 2 DZCP trend (−0.57±4 %) but also a large spread about the mean bias that suggests that the CO 2 DZCP is more susceptible to aliasing.However, the use of t CO 2 = 25 % PU resulted in a relatively small mean bias (−14 ± 2 %) with high consistency (r 2 = 0.94).Although the mean bias was less in the MCS for the CO 2 UZCP (−1±3 %), it was less consistent (r 2 = 0.80).The t CO 2 UZCP had a mean bias of −23±1 % (r 2 = 0.97).Differences between the spring and autumn phase biases calculated from CO 2 and t CO 2 phase metrics carry through to the respective CUP estimates, where the t CO 2 CUP had a mean bias of −28 ± 1 % (r 2 = 0.93) relative to a bias of −55±1 % (r 2 = 0.45) in the CO 2 CUP.Estimates of t CO 2 phase metrics tended to be more consistent, and while it resulted in significantly more accurate estimates of the trend in spring phase, the autumn phase was better represented by the CO 2 UZCP.We expect that this is a result of the asymmetry of the high-latitude CO 2 seasonal cycle.Analysis of peak rates of uptake and release resulted in mean biases of −18 ± 2 and −28 ± 2 % for PU and PR respectively.In general, the trend estimates from the analysis had the correct sign so long as the trend was sufficiently large (> 0.25 % yr −1 for PU and PR, and > 0.1 days yr −1 for changes in phase).The CO 2 phase metric trend estimates were the most likely to have the wrong sign compared to the t CO 2 phase metrics by a factor of 4.5, 4, and 1.5 for the DZCP, UZCP, and CUP respectively.Finally, the t CO 2 metrics were far more effective at detecting statistically significant trends where the CO 2 metrics typically missed 33-50 % of them.
Figure B8 shows a regression of the linear trend in integrated CO 2 uptake and release against the estimated seasonal amplitude from the individual MCS runs.We find that the linear trends in annually integrated CO 2 uptake (ppm yr −2 ) are highly correlated with the amplitude trend (ppm yr −1 ), but the amplitude trends are poorly correlated with changes in integrated release of CO 2 .Previous work has shown that this is due to the rapid temporal change in CO 2 during the period of net carbon uptake relative to the more gradual release of CO 2 outside of this period (Graven et al., 2013).

Appendix C: Analysis of detrended CO 2 seasonal cycle
Our analysis of phase changes in the CO 2 seasonal cycle at BRW shows a much tighter coupling between the timing of the downward and upward zero crossing points with values of −0.20 days yr −1 (p < 0.01) and −0.18 days yr −1 (p < 0.05) respectively.This results in a more conserved carbon uptake period, with a coefficient of 0.02 days yr −1 (p > 0.1), which is consistent with the ecosystem having an intrinsic or fixed uptake period (not shown).Recent work using changes in CO 2 has reported a change of −0.17 days yr −1 for the downward phase over a similar time period (Graven et al., 2013).Although we have included these values for completeness, we have already shown that there are significant weaknesses in using detrended CO 2 as opposed to t CO 2 for this analysis.
Table D1.Estimated trends of downward and upward zero crossing points (DZCP and UZCP respectively), peak uptake and release (PU and PR respectively), and carbon uptake period (CUP) calculated from CO 2 and t CO 2 data for seven high-latitude measurement sites (Fig. 1).The 95 % confidence intervals and p values are calculated for each trend estimate.in TGS LEN suggests that the potential period during which plant growth is not hindered by low temperatures has been significantly extended by approximately 11 days (> 45 • N) and 20 days (> 60 • N) since 1970, consistent with previous findings (Linderholm, 2006;Barichivich et al., 2012).Table E2 shows the relationship between northern high-latitude land surface temperature anomalies with the BRW CO 2 and t CO 2 phase metrics throughout 1973-2012.We find there are significantly different results depending on whether CO 2 and t CO 2 phase metrics are used.
The warming-induced earlier onset of springtime carbon uptake is also supported by observed increases in vegetation greenness described by NDVI inferred from space-borne sensors (Gong and Shi, 2003;Mao et al., 2012;Cong et al., 2013).Increases in autumn NDVI have also been observed and while this is indicative of increased photosynthetic activity it is not necessarily inconsistent with the observed early onset of net carbon release.This is because it does not provide information about respiration processes.Our analysis of NDVI data (not shown) finds an increases of vegetation greenness in spring and autumn have led to significant lengthening of the photosynthetic growing season over the measurement period, where autumn greening is changing in most regions at a greater rate than spring greening.

E2 δ 13 C data
Figure 4 shows δ 13 C data over CBA, with the corresponding CO 2 mole fraction data.Measurements of δ 13 C show a strong seasonal variation, which is anti-correlated with CO 2 .Plants preferentially take the lighter carbon 12 C isotope out of the atmosphere through photosynthesis during spring and summer resulting in an increase in δ 13 C and release more 12 C than 13 C during autumn and winter resulting in a decrease in δ 13 C.
Figure 9 shows a similar phase analysis for (−1) × δ 13 C and (−1) × δ 13 C, comparing it with variability and trends with the corresponding CO 2 values.Table E3 shows regression coefficients and mean statistics for the spring and autumn phase and the CUP.We find that at least 68 % of the observed trend in CO 2 DZCP and UZCP can be explained by variations in colocated measurements of δ 13 C.This suggests that the terrestrial biosphere is largely responsible for observed CO 2 variability with the remainder due to atmospheric transport and other minor source variations.This result is consistent with previous work (Graven et al., 2013) that showed using an atmospheric transport model that atmospheric transport variations contributed < 7 % of the observed variation in CO 2 seasonal amplitudes at high northern latitudes.
Figure1.The NOAA/ESRL monitoring sites used in our CO 2 time series analysis.For the seasonal cycle analysis, we use the northern hemispheric high-latitude sites (blue).The sites shown in green, red, and magenta are used for growth rate analysis only.The six sites with a black border are those with the longest time series in each 30 • latitude band.The shaded regions are the temperate and boreal northern Hemispheric land regions defined in the initial Transcom study, and which we use for analysis of NDVI, temperature, and atmospheric transport.

Figure 2 .
Figure2.The maximum CO 2 perturbations caused by biosphere carbon fluxes from five Transcom land regions (Fig.1) to the zonal mean concentrations over the high-, mid-, and low-latitude Northern Hemisphere averaged over 2004-2009.These values were determined by using the GEOS-Chem atmospheric transport model (see main text for further details).The error bars denote the 1σ of the year-to-year variability over the 6-year period.The zonal means are defined as the mean of the grid points sampled nearest to the sites shown in Fig.1.Model output from December is missing.

Figure 4 .
Figure 4. Weekly (top) CO 2 mole fraction (ppm) measurements (black) and (bottom) δ 13 C values (‰) at Cold Bay, Alaska (CBA; 55.2 • N, 162.7 • W), from 1980 to 2012.Imputed values, shown in red, are inferred from a locally averaged seasonal cycle adjusted for anomalies in growth rate.Any remaining missing values are extracted from a fitted piecewise cubic spline curve (magenta).

Figure 5 .
Figure5.Decadal mean CO 2 growth rates inferred from individual site measurements and averaged in 20 • latitude bins having retained (left) and subtracted (right) the decadal mean global fossil fuel emissions (CDIAC).The solid line with error bars represents the decadal mean growth rate in each latitude bin with ±1σ representing the standard deviation between individual sites in that latitude bin.The global decadal mean growth rate is indicated by the dashed lines and mean values with ±1σ representing the standard deviation between all sites.Values for MLO, which are typically taken to be representative of the global growth rate, are highlighted with a circle.Time series of annual growth rates were determined for individual CO 2 measurement sites before calculating decadal mean growth rates and then binning the decadal mean growth rates into 20 • latitude bins.We subtract a global mean growth rate due to fossil fuel combustion from all sites.

Figure 6 .
Figure6.A schematic describing the metrics we use to characterize changes in the amplitude and phase of atmospheric CO 2 (ppm).In this example we use detrended annual and semiannual components of CO 2 data from Barrow, Alaska.

Figure 7 .
Figure7.Scatterplots of the t CO 2 =25 % PU (spring phase) and t CO 2 UZCP (autumn phase) (days) at four high northern latitude sites (see main text).The coloured lines show the trajectory of the 2-year running mean of the scatterplot, where colours represent the year of measurement.

Figure 8 .
Figure 8.Time series of the percentage change of peak uptake and release at four high northern latitude sites (see main text).Each panel shows the data as blue closed circles and the 25 % uncertainty interval.The dashed black line is the fitted linear trend that is reported inset of each panel.

Figure 9 .
Figure9.Time series of phase changes and the corresponding change to the carbon uptake period of δ 13 C and CO 2 (left) and t δ 13 C and t CO 2 (right), expressed as days.The red line is the 2-year running mean.

Figure A1 .
Figure A1.Top row: weekly mean (black) and low-pass filtered (red, periods > 18 months) CO 2 mole fraction time series (ppm) at MaunaLoa, 1959-2012.Middle row: (left)  the wavelet power spectrum where the colour scale is log (power).The black solid lines denotes the cone of influence.The power spectrum tends to emphasise very low frequency information so we have subtracted an exponential term prior to applying the wavelet transform to emphasise the high frequency variability (right) the corresponding time-integrated power spectrum.Bottom row: the inferred annual growth rate of CO 2 (ppm yr −1 ).

Figure B1 .
Figure B1.Synthetic CO 2 "flux" (left), expressed as ppm week −1 over an annual cycle, and (right) the corresponding mole fraction (ppm) time series repeated over successive years.The CO 2 annual cycle is based on the observed cycle at Barrow Alaska.

Figure B2 .
Figure B2.analysis of t CO 2 flux time series including a prescribed earlier onset of net CO 2 uptake.Top left panel: the defined flux time series and the associated detrended time series.Top right panel: the expected (defined) and actual change in peak uptake and release of CO 2 .Bottom panels: the expected (defined) and actual change in (left) DZCP and (right) UZCP, including an operational version of the phase metric as described in the main text.

Figure B3 .
Figure B3.Same as Fig. B2 but including an earlier autumn onset of net CO 2 release.

Figure B4 .Figure B5 .
Figure B4.Same as Fig.B2but introducing a trend of 0.75 % yr −1 trend in the peak uptake and years of anomalously high and low uptake respectively.

Figure B6 .
Figure B6.Probability densities of trends introduced as part of a 1000-member ensemble of synthetic time series generated for the Monte Carlo experiment where the black line is the fitted probability distribution.

Table 1 .
Parameters used by the control wavelet transform for monthly and weekly spectral decomposition of CO 2 mole fraction.

Table 2 .
Global decadal mean growth rates (GR) and the corresponding growth rate due to fossil fuel combustion (FF) and natural sources (GR-FF).Units are ppm yr −1 .

Table E1 .
Temperature linear trend analysis (1970Temperature linear trend analysis ( -2011) )for the beginning and end of the thermal growing season.TGS BEG (days decade −1 ) Spring T ( • C decade −1 )

Table E2 .
Linear regression coefficients that describe the relationship between changes in CO 2 , t CO 2 , and temperature phase metrics at different latitude bands in the high northern latitudes.

Table E3 .
Linear regression coefficients that describe the relationship between changes in CO 2 , t CO 2 , and δ 13 C at BRW during the overlapping time span of the data.