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

Abstract

Atmospheric variations of carbon dioxide (CO${}_{\mathrm{2}}$) mole fraction reflect changes in atmospheric transport and regional patterns of surface emission and uptake. Here we present a study of changes in the observed high northern latitude CO${}_{\mathrm{2}}$ seasonal cycle. We report new estimates for changes in the phase and amplitude of the seasonal variations, indicative of biospheric changes, by spectrally decomposing multi-decadal records of surface CO${}_{\mathrm{2}}$ mole fraction using a wavelet transform to isolate the changes in the observed seasonal cycle. We also perform similar analysis of the first derivative of CO${}_{\mathrm{2}}$ mole fraction,${\mathrm{\Delta }}_t$CO${}_{\mathrm{2}}$, that is a crude proxy for changes in CO${}_{\mathrm{2}}$ flux. Using numerical experiments, we quantify the aliasing error associated with independently identifying trends in phase and peak uptake and release to be 10–25 %, with the smallest biases in phase associated with the analysis of${\mathrm{\Delta }}_t$CO${}_{\mathrm{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\pm 0.02$ $\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{1}}$, which is consistent with previous work. Using ${\mathrm{\Delta }}_t$CO${}_{\mathrm{2}}$ we determine estimates for the timing of the onset of net uptake and release of CO${}_{\mathrm{2}}$ of $-0.14\pm 0.14$ and $-0.25\pm 0.08$ $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ respectively and a corresponding net uptake period of $-0.11\pm 0.16$ $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$, which are significantly different to previously reported estimates. We find that the wavelet transform method has significant skill in characterizing changes in the peak uptake and release. We find a trend of $0.65\pm 0.34$ $\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$ ($p<0.01$) and $0.42\pm 0.34$ $\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$ ($p<0.05$) for rates of peak uptake and release respectively. Our analysis does not provide direct evidence about the balance between uptake and release of carbon when integrated throughout the year, but the increase in the seasonal amplitude of CO${}_{\mathrm{2}}$ together with an invariant net carbon uptake period provides evidence that high northern latitude ecosystems are progressively taking up more carbon during spring and early summer.

Introduction

Combustion of fossil fuel and cement production represent the dominant annual source of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$. There is also a minor source from the combustion of biomass and a diffuse source from the emissions and oxidation of reduced carbon . On an annual basis approximately 50 % of those emissions remain in the atmosphere with the remainder taken up by the land and ocean . Regional changes to the net biospheric flux of ${\mathrm{CO}}_{\mathrm{2}}$, and consequent changes in atmospheric ${\mathrm{CO}}_{\mathrm{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. ), has reported substantial increases in the amplitude of the seasonal exchange of ${\mathrm{CO}}_{\mathrm{2}}$ since the 1950s, particularly at mid- to high northern latitudes . Here, we use the wavelet transform to isolate changes in the ${\mathrm{CO}}_{\mathrm{2}}$ seasonal cycle, revealing new insights about the growth rate, and changes in the amplitude and phase of ${\mathrm{CO}}_{\mathrm{2}}$ associated with the growing season.

Many previous studies have used atmospheric measurements of ${\mathrm{CO}}_{\mathrm{2}}$ to analyse changes in the seasonal cycle. These studies have typically employed curve-fitting techniques e.g. and spectral filtering methods such as complex demodulation . The most commonly used method is a combination of curve-fitting and spectral filtering (e.g. ). 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 ${\mathrm{CO}}_{\mathrm{2}}$ time series into seasonal, long-term, and residual components while retaining information about phase and amplitude in the seasonal cycle . We show through extensive analysis of synthetic time series that using the time derivative of ${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$, provides more accurate estimates of changes in the phase of the ${\mathrm{CO}}_{\mathrm{2}}$ seasonal cycle and that the common use of the zero crossing points of the detrended ${\mathrm{CO}}_{\mathrm{2}}$ seasonal cycle is flawed. We show that we can also faithfully reproduce changes in the rates of peak uptake and peak release of ${\mathrm{CO}}_{\mathrm{2}}$, allowing us to understand observed changes in the amplitude of the seasonal cycle.

In the next section we describe measurements of ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction, the isotope ratio ${\mathit{\delta }}^{13}\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$. In Sect. 4, we present our analysis of ${\mathrm{CO}}_{\mathrm{2}}$ growth rates and changes in the phase and amplitude of the ${\mathrm{CO}}_{\mathrm{2}}$ seasonal cycle. We conclude our paper in Sect. 5.

Data

${\mathrm{CO}}_{\mathrm{2}}$ mole fraction data

Figure  shows the sites from the National Oceanic and Atmospheric Administration (NOAA) Global Greenhouse Gas Reference Network (GGGRN,), which include at least 15 $\mathrm{years}$ of ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction data. We focus on high northern latitude sites where (a) seasonal contributions of ${\mathrm{CO}}_{\mathrm{2}}$ are predominantly driven by boreal vegetation and (b) contributions to observed ${\mathrm{CO}}_{\mathrm{2}}$ from continents at these latitudes are approximately equal (Fig. ). We report our ${\mathrm{CO}}_{\mathrm{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.

The NOAA/ESRL monitoring sites used in our ${\mathrm{CO}}_{\mathrm{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${}^{\circ }$ 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 omitted. See PDF]

The maximum ${\mathrm{CO}}_{\mathrm{2}}$ perturbations caused by biosphere carbon fluxes from five Transcom land regions (Fig. ) 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$\mathit{\sigma }$ 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. . Model output from December is missing.

[Figure omitted. See PDF]

Twin air samples are collected weekly at the sites and analysed for ${\mathrm{CO}}_{\mathrm{2}}$ at NOAA Earth System Research Laboratory (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) X${\mathrm{CO}}_{\mathrm{2}}$ scale to working standards ($\pm$0.03 $\mathrm{ppm}$; ), the analytical repeatability of the analysers for a sample measurement ($\pm$0.03 $\mathrm{ppm}$), and the agreement between pairs of samples collected within 20 s of one another ($\pm$0.1 $\mathrm{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 ${\mathrm{CO}}_{\mathrm{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 7-year average seasonal cycle (calculated using the 3 $\mathrm{years}$ on either side of and including the year of interest) to a value from a deseasonalised reference time series (Fig. ), which accounts for large-scale anomalies in the growth rate in 30${}^{\circ }$ 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  shows an example of our imputation approach using the ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction and ${\mathit{\delta }}^{13}\mathrm{C}$ time series from Cold Bay, Alaska.

Reference weekly ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction measurements (ppm) covering various time spans for 30${}^{\circ }$ latitude bands used to impute missing data points. BRW, NWR, MLO, SMO, CGO, and SPO are codes to denote Barrow (71.3${}^{\circ }$ N, 156.6${}^{\circ }$ W), Niwot Ridge (40.0${}^{\circ }$ N, 105.6${}^{\circ }$ W), Mauna Loa (19.5${}^{\circ }$ N, 155.6${}^{\circ }$ W), American Samoa (14.2${}^{\circ }$ S, 170.5${}^{\circ }$ W), Cape Grim (40.7${}^{\circ }$ S, 144.7${}^{\circ }$ E), and South Pole (89.9${}^{\circ }$ S, 24.8${}^{\circ }$ W).

[Figure omitted. See PDF]

Weekly (top) ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction (ppm) measurements (black) and (bottom) ${\mathit{\delta }}^{13}\mathrm{C}$ values ($\mathrm{‰}$) at Cold Bay, Alaska (CBA; 55.2${}^{\circ }$ N, 162.7${}^{\circ }$ 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 omitted. See PDF]

${\mathit{\delta }}^{13}{\mathrm{CO}}_{\mathrm{2}}$ data

We also use measurements of ${\mathit{\delta }}^{13}\mathrm{C}$ that are colocated with the ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction data. The isotope samples are analysed at the Stable Isotope Laboratory at The Institute of Arctic and Alpine Research using flasks of air provided by the NOAA GGGRN. These help us to attribute observed changes in ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction to land biospheric uptake and release. The ratio ${\mathit{\delta }}^{13}\mathrm{C}$ is defined as $${\mathit{\delta }}^{13}\mathrm{C}=\left[{\displaystyle \frac{{\left[\frac{{}^{13}\mathrm{C}}{{}^{12}\mathrm{C}}\right]}_{\text{sample}}}{{\left[\frac{{}^{13}\mathrm{C}}{{}^{12}\mathrm{C}}\right]}_{\text{standard}}}},-,\mathrm{1}\right]\times 1000,$$ where ${\left[\frac{{}^{13}\mathrm{C}}{{}^{12}\mathrm{C}}\right]}_{\text{sample}}$ is the ratio of ${}^{13}\mathrm{C}$ to ${}^{12}\mathrm{C}$ (mol/mol) within the sample, and ${\left[\frac{{}^{13}\mathrm{C}}{{}^{12}\mathrm{C}}\right]}_{\text{sample}}$ is the ratio of ${}^{13}\mathrm{C}$ to ${}^{12}\mathrm{C}$ of the internationally accepted Pee Dee Belemnite standard. Individual measurements of ${}^{12}\mathrm{C}$ and ${}^{13}\mathrm{C}$ are determined by isolating the ${\mathrm{CO}}_{\mathrm{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 to help interpret observed variations in the ${\mathrm{CO}}_{\mathrm{2}}$ time series. These data have a ${0.5}^{\circ }\times {0.5}^{\circ }$ 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) . 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 $\mathrm{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 ${\mathrm{CO}}_{\mathrm{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 ${\mathit{\psi }}_{\mathrm{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 $$W_n\left(s\right)=\sum _{k=\mathrm{0}}^{N-\mathrm{1}}{\widehat{x}}_k\widehat{\mathit{\psi }}\cdot \left(s{\mathit{\omega }}_k\right)e^{i{\mathit{\omega }}_{k}n\mathit{\delta }t},$$ where ${\widehat{x}}_k$ is the discrete Fourier transform of $x_n$, $N$ is the number of points in the time series, $k=\mathrm{0}\mathrm{\dots }N-\mathrm{1}$ is the frequency index, and $\widehat{\mathit{\psi }}\cdot \left(s{\mathit{\omega }}_k\right)$ is the complex conjugate of the Fourier transform of a normalized, scaled, and translated version of ${\mathit{\psi }}_{\mathrm{0}}\left(\mathit{\eta }\right)$, where $s$ is the scale and ${\mathit{\omega }}_k$ is the angular frequency. We use the Morlet wavelet, a plane wave modulated by a Gaussian envelope: $${\mathit{\psi }}_{\mathrm{0}}\left(\mathit{\eta }\right)={\mathit{\pi }}^{\frac{-\mathrm{1}}{\mathrm{4}}}{e}^{i{\mathit{\omega }}_{\mathrm{0}}\mathit{\eta }}{e}^{\frac{-{\mathit{\eta }}^{\mathrm{2}}}{\mathrm{2}}},$$ where ${\mathit{\omega }}_{\mathrm{0}}$ is the nondimensional frequency and $\mathit{\eta }$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction time series. The wavelet is comprised of a real and imaginary part, providing information about amplitude and phase respectively.

Parameters used by the control wavelet transform for monthly and weekly spectral decomposition of ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction.

Parameter	$\mathit{\delta }t=\frac{\mathrm{1}}{12}$	$\mathit{\delta }t=\frac{\mathrm{1}}{52}$
$\mathit{\delta }j$	0.25	0.01
$s_{\mathrm{0}}$	2$\mathit{\delta }t$	$\mathit{\delta }t$
$C_{\mathit{\delta }}$	0.7785	0.7785
${\mathit{\psi }}_{\mathrm{0}}$	${\mathit{\pi }}^{-\frac{\mathrm{1}}{\mathrm{4}}}$	${\mathit{\pi }}^{-\frac{\mathrm{1}}{\mathrm{4}}}$

We can recover the original time series from wavelet space using the corresponding inverse transform 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:

$${\widehat{W}}_n={\displaystyle \frac{\mathit{\delta }j\mathit{\delta }{t}^{\frac{\mathrm{1}}{\mathrm{2}}}}{C_{\mathit{\delta }}{\mathit{\psi }}_{\mathrm{0}}\left(\mathrm{0}\right)}}\sum _{j=\mathrm{0}}^J{\displaystyle \frac{\mathrm{\Re }\left\{W_n,\left(,s_j,\right)\right\}}{s_j^{\frac{\mathrm{1}}{\mathrm{2}}}}},$$ where ${\mathit{\psi }}_{\mathrm{0}}\left(\mathrm{0}\right)$ removes the energy scaling and $s_j^{\frac{\mathrm{1}}{\mathrm{2}}}$ converts the wavelet transform to an energy density. $C_{\mathit{\delta }}$ and ${\mathit{\psi }}_{\mathrm{0}}\left(\mathrm{0}\right)$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ time series with a prescribed trend. We find that the value for $C_{\mathit{\delta }}$ previously reported introduces a small trend in the original minus reconstructed residual and find that $C_{\mathit{\delta }}=0.7785$ results in a much smaller, unbiased residual with a typical value $<0.05$ $\mathrm{ppm}$ for monthly data and $<0.002$ $\mathrm{ppm}$ for weekly data (not shown). Table  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 ${\mathrm{CO}}_{\mathrm{2}}$ and vice versa. However, this is a common weakness of any method used to decompose such time series.

Decadal mean ${\mathrm{CO}}_{\mathrm{2}}$ growth rates inferred from individual site measurements and averaged in 20${}^{\circ }$ 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 $\pm$1$\mathit{\sigma }$ 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 $\pm$1$\mathit{\sigma }$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ measurement sites before calculating decadal mean growth rates and then binning the decadal mean growth rates into 20${}^{\circ }$ latitude bins. We subtract a global mean growth rate due to fossil fuel combustion from all sites.

[Figure omitted. See PDF]

We find that for atmospheric ${\mathrm{CO}}_{\mathrm{2}}$, the wavelet power spectrum peaks at periods (reciprocal of frequency) of 6 and 12 months (Appendix ), 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 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 ).

Results

Growth rates

Figure  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 ${\mathrm{CO}}_{\mathrm{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 $\sim 0.35$ $\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ growth rates at Mauna Loa are within a fraction of a percent of NOAA values.

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${}^{-\mathrm{1}}$.

Decade	No. sites	FF	GR	GR 1$\mathit{\sigma }$	GR–FF
1960–1969	1	1.51	0.86	n/a	$-$0.65
1970–1979	2	2.25	1.21	0.055	$-$1.04
1980–1989	13	2.61	1.58	0.108	$-$1.03
1990–1999	38	3.02	1.48	0.056	$-$1.54
2000–2009	49	3.79	1.90	0.076	$-$1.89

By subtracting anthropogenic fossil fuel emission estimates from the atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ signal (Table ) 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 . We find that during the 1980s the net annual uptake by the terrestrial biosphere and ocean was typically $-1.03\pm 0.11$ $\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{1}}$ when averaged across all sites and where the uncertainty is equal to 1$\mathit{\sigma }$. This rate increases dramatically in the 1990s to approximately $-1.54\pm 0.06$ and to $-1.89\pm 0.08$ $\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{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\pm 18.2$ % (including only fossil fuel) and $44.1\pm 14.4$ % (including fossil fuel and land use change), consistent with previous work . For the purpose of the following calculations we have removed the annual growth rate from the observed ${\mathrm{CO}}_{\mathrm{2}}$ concentrations, following the method described in Appendix .

Phase and amplitude analysis

We use several metrics to interpret the ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$, a proxy for the net flux of ${\mathrm{CO}}_{\mathrm{2}}$. As we discuss below and report in Appendix , analysis of ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ leads to less biased estimates for trends in the phase of the ${\mathrm{CO}}_{\mathrm{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  shows, using example data from BRW, how the detrended ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ variations are related. The amplitude of the seasonal cycle, defined as the peak-to-peak difference (maxima minus minima) of the seasonal ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction time series, has been used in previous studies as a measure of biological activity e.g.. 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 ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ we define three periods during an annual cycle: (1) an uptake period when ${\mathrm{\Delta }}_t{\mathrm{CO}}_{\mathrm{2}}<\mathrm{0}$ and there is a net negative ${\mathrm{CO}}_{\mathrm{2}}$ flux to the atmosphere (photosynthesis is higher than respiration); (2) a release period when ${\mathrm{\Delta }}_t{\mathrm{CO}}_{\mathrm{2}}>\mathrm{0}$ and there is a net source of ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ is typically small (but non-zero due to transport of ${\mathrm{CO}}_{\mathrm{2}}$ from the lower latitudes).

A schematic describing the metrics we use to characterize changes in the amplitude and phase of atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ (ppm). In this example we use detrended annual and semiannual components of ${\mathrm{CO}}_{\mathrm{2}}$ data from Barrow, Alaska.

[Figure omitted. See PDF]

To look at changes in phase, previous studies have used the zero-crossing point (ZCP) of ${\mathrm{CO}}_{\mathrm{2}}$ which refers to times when the detrended seasonal cycle is equal to 0 e.g.. 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 ${\mathrm{CO}}_{\mathrm{2}}$ concentration can only be estimated from the detrended seasonal cycle. The long-term increase in ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{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 aliasing error. As the first derivative is closely related to the actual flux, it is less affected by this shifting up or down of the seasonal cycle relative to the zero-line. This is shown in Figs. , , and .

The beginning of the period of net carbon uptake is difficult to determine accurately using the seasonal cycle at high-latitude sites because small mole fraction variations during the dormant period (which has a near-zero flux) are sufficient to bring ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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 ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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 ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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 ). These errors are not unique to using the wavelet transform but are more pronounced when using the detrended ${\mathrm{CO}}_{\mathrm{2}}$ seasonal cycle as opposed to ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$. To our knowledge, no previous study has quantified these errors when estimating phase changes in the ${\mathrm{CO}}_{\mathrm{2}}$ seasonal cycle. We generally find that analysis of ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ produces more reliable and less biased estimates than ${\mathrm{CO}}_{\mathrm{2}}$ trend estimation of either phase with an estimated 25 % systematic aliasing error (Appendix ). Unless explicitly stated all subsequent results will refer to our analysis ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$. We also find that we can capture at least 80 % of independent trends in the PU and PR of the ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ (Appendix ).

Analysis of NOAA ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction data

Figure  shows that changes in spring and autumn phases determined from BRW ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ are $-$0.14 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ ($p<0.05$) and $-$0.25 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ ($p<0.01$) respectively, with a corresponding CUP change of $-$0.11 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ ($p>0.1$); the analysis of the other study sites is shown in Appendix . We find no evidence using phase changes of ${\mathrm{CO}}_{\mathrm{2}}$ or ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ for a significant change in CUP throughout the measurement period due to a simultaneous advance of the spring and autumn phase.

Scatterplots of the ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$=25 % PU (spring phase) and ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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 omitted. See PDF]

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 omitted. See PDF]

Time series of phase changes and the corresponding change to the carbon uptake period of ${\mathit{\delta }}^{13}\mathrm{C}$ and ${\mathrm{CO}}_{\mathrm{2}}$ (left) and ${\mathrm{\Delta }}_t$${\mathit{\delta }}^{13}\mathrm{C}$ and ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ (right), expressed as days. The red line is the 2-year running mean.

[Figure omitted. See PDF]

The concomitant observed changes in ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ and in ${\mathrm{\Delta }}_t$${\mathit{\delta }}^{13}\mathrm{C}$ (Fig. , also discussed in Appendix ) supports the idea that observed ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ (Appendix ). We find the start of the thermal growing season (defined as the continuous period above 5 ${}^{\circ }$C, Appendix ) is advancing 2 (3) times faster at latitudes $>45$${}^{\circ }$ N ($>60$${}^{\circ }$ N), which agrees with previous studies e.g.. 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 (1982–2006) 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 ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ time series shows more uptake of ${\mathrm{CO}}_{\mathrm{2}}$ in spring and early summer and earlier onset of net release of ${\mathrm{CO}}_{\mathrm{2}}$ between mid-summer and autumn. A number of studies have linked increases in NDVI and subsequent carbon uptake with a ${\mathrm{CO}}_{\mathrm{2}}$ fertilization effect 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\pm 0.02$ $\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{1}}$) are consistent in percentage terms with previous work over the same time period . We find that the observed change in amplitude at BRW is partially due to an increase in PR ($0.42\pm 0.34$ $\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{1}}$, $p>0.05$) and a larger increase in PU ($0.65\pm 0.34$ $\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{1}}$, $p<0.01$). Figure  shows that statistically significant trends ($p<0.05$) in PU are observed at five of the seven high-latitude sites (ALT, BRW, CBA, ICE and ZEP, Table ). In most of these cases, the change in PU is significantly larger than the change in PR. We show in Fig.  that trends in amplitude are determined mainly by changes in uptake during the CUP (Appendix ). Previous analysis of these data has shown that changes in atmospheric transport cannot explain changes in the amplitude .

Concluding remarks

We have used a wavelet transform to spectrally isolate changes in the seasonal cycle of atmospheric ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ are predominantly driven by boreal vegetation and (b) contributions to observed ${\mathrm{CO}}_{\mathrm{2}}$ from continents at these latitudes are approximately equal.

We found that the atmospheric growth rate of ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ concentrations. This results in a near-constant airborne ${\mathrm{CO}}_{\mathrm{2}}$ fraction of $55.8\pm 18.2$ % (including only fossil fuel) and $44.1\pm 14.4$ % (including fossil fuel and land use change), consistent with previous studies.

Using the detrended ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ provided more accurate estimates of independent changes in phase and peak uptake and release of ${\mathrm{CO}}_{\mathrm{2}}$, to within 10–25 % of the “true” values.

We reported an increase in amplitude of $0.09\pm 0.02$ $\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{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 ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$, taking the timing of uptake reaching 25 % of peak uptake as the beginning of the CUP, and the timing of ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ switching to positive as the end of the CUP. These phase thresholds take into account that observed ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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\pm 0.14$ and $-$0.25 $\pm$ 0.08 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ respectively and a corresponding revision to the length of the net uptake period of $-0.11\pm 0.16$ $\mathrm{days}{\mathrm{yr}}^{-\mathrm{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 ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$. A more frequent inspection of these data using advanced statistical tools such as the wavelet transform also has a role to play.

Example of spectral decomposition

Figure  shows, as an example, the spectral decomposition of ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ time series with additional synthetic data so we are able to analyse the entire ${\mathrm{CO}}_{\mathrm{2}}$ time series (Sect. ). We find that most of the power is in the annual and semiannual periods, as expected, but also peaks in power at period $>\mathrm{1}$ year but this is likely a result of responses of the ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction data. Figure  shows two example applications: (1) as a lowpass filter to deseasonalize the ${\mathrm{CO}}_{\mathrm{2}}$ data (removing periods $<18$ months) and (2) the associated annual growth rate ($\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{1}}$), which we find is within $<0.1$ $\mathrm{ppm}$ of the reported values from NOAA (not shown).

Top row: weekly mean (black) and low-pass filtered (red, periods $>18$ months) ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction time series (ppm) at Mauna Loa, 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 ${\mathrm{CO}}_{\mathrm{2}}$ ($\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{1}}$).

[Figure omitted. See PDF]

Error characterization of phase and amplitude estimates

We use synthetic ${\mathrm{CO}}_{\mathrm{2}}$ time series data, defined with specific changes in amplitude and phase, to characterize aliasing errors due to application of the wavelet transform of ${\mathrm{CO}}_{\mathrm{2}}$ concentration data or its first derivative (${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$). Insights from this synthetic analysis are directly applied to our interpretation of NOAA mole fraction measurements in the main paper.

Synthetic model framework

We use a simple box model based on the ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction time series at Barrow, Alaska (Fig. ). BRW is the most suitable site for this purpose because is has a long time series and as it is representative of high-latitude ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{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 $\mathrm{years}$ (equivalent to the time span of the BRW time series) and integrated to convert the flux to ${\mathrm{CO}}_{\mathrm{2}}$ concentration. For our experiments, described below, we introduce trends and variability to various aspects of ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ before integrating with respect to time to recover ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction. Detrending is as described in the main paper.

Synthetic ${\mathrm{CO}}_{\mathrm{2}}$ “flux” (left), expressed as $\mathrm{ppm}{\mathrm{week}}^{-\mathrm{1}}$ over an annual cycle, and (right) the corresponding mole fraction (ppm) time series repeated over successive years. The ${\mathrm{CO}}_{\mathrm{2}}$ annual cycle is based on the observed cycle at Barrow Alaska.

[Figure omitted. See PDF]

Numerical experiments

The starting point of our numerical experiments is the detrended time series of atmospheric ${\mathrm{CO}}_{\mathrm{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.

Wavelet analysis of ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ flux time series including a prescribed earlier onset of net ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{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 omitted. See PDF]

The following three broad set of experiments are designed to identify the best metrics to describe changes in the contemporary cycle from detrended ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$) or release (positive ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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.

Perturbing the timing of the spring and autumn phases

Figure  shows the results of our analysis of a time series for which we introduced a progressively earlier onset of net ${\mathrm{CO}}_{\mathrm{2}}$ uptake of 0.50 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ for ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ DZCP. The ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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 ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$. The ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ metrics were found to be better at capturing the springtime trend to within 23 and 16 % respectively when compared with the equivalent ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction metric which underestimates the trend by 63 %. This has implications for using the CO${}_{\mathrm{2}}$ mole fraction ZCPs to interpret changes in the phase. There is little change in any of the UZCP metrics (typically $<0.025$ $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$) as a result of aliasing. The wavelet detrending introduces a $-$0.01 $\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$ trend in peak ${\mathrm{CO}}_{\mathrm{2}}$ uptake and a concurrent increase in peak ${\mathrm{CO}}_{\mathrm{2}}$ release of 0.14 $\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$ corresponding to $-$0.4 and 5.6 % across the 40-$\mathrm{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.

Same as Fig.  but including an earlier autumn onset of net ${\mathrm{CO}}_{\mathrm{2}}$ release.

[Figure omitted. See PDF]

Figure  shows the same calculation but for introducing an earlier autumn onset of net ${\mathrm{CO}}_{\mathrm{2}}$ release of 0.30 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{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 ${\mathrm{\Delta }}_t{\mathrm{CO}}_{\mathrm{2}}=\mathrm{0}$ metrics had non-zero trends up to $\sim -$0.16 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$. All three UZCP phase metrics underestimate the change in the defined phase change by amounts ranging from 11 to 22 % where the ${\mathrm{CO}}_{\mathrm{2}}$ UZCP performed the best. The earlier onset of net ${\mathrm{CO}}_{\mathrm{2}}$ release aliases into a 2.5 % increase in peak ${\mathrm{CO}}_{\mathrm{2}}$ release and a 5 % increase in peak ${\mathrm{CO}}_{\mathrm{2}}$ across the entire time series.

Perturbing the magnitude of net uptake and release of ${\mathrm{CO}}_{\mathrm{2}}$

Figure  shows the results of introducing a progressive enhancement of ${\mathrm{CO}}_{\mathrm{2}}$ uptake of roughly 0.70 $\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$, equivalent to a 28 % increase over 40 $\mathrm{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 ${\mathrm{CO}}_{\mathrm{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.

Same as Fig.  but introducing a trend of 0.75 $\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$ trend in the peak uptake and years of anomalously high and low uptake respectively.

[Figure omitted. See PDF]

Same as Fig.  but introducing simultaneous trends in spring and autumn phase and in the peak amplitude and release of ${\mathrm{CO}}_{\mathrm{2}}$. We also superimpose Gaussian random noise to describe interannual variation.

[Figure omitted. See PDF]

We find that the wavelet transform attributes the 0.70 $\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$ increased uptake as 0.59 $\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$ uptake and 0.20 $\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$ release. The mole fraction metrics infer non-zero UZCP and DZCP phase changes of 0.06 and 0.16 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ respectively, while the 25 % ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ UZCP is the most affected out of the phase metrics. In addition, information from the exceptional years of uptake is aliased into the ${\mathrm{CO}}_{\mathrm{2}}$ UZCP and is spread over a number of years rather than just one. This is not the case for the ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ metrics indicating that they are better for estimating interannual variability.

Simultaneous variations in phase and peak uptake and release

Figure  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 ${\mathrm{\Delta }}_t{\mathrm{CO}}_{\mathrm{2}}=\mathrm{0}$, and the PU and PR. We also superimpose Gaussian random noise within $\pm$10 days and $\pm$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 ${\mathrm{CO}}_{\mathrm{2}}$ uptake ($-$0.02 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$) which is captured by the ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ phase metric (0.02 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$). The ${\mathrm{CO}}_{\mathrm{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 ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ phase metric (0.05 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$) underestimates the expected trend (0.09 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$) by $\sim 45$ %, while the ${\mathrm{CO}}_{\mathrm{2}}$ UZCP overestimates it by a factor of 2.8. The estimated trend in PU is 0.54 % yr${}^{-\mathrm{1}}$ which is 80 % of the expected trend (0.68 % yr${}^{-\mathrm{1}}$), while the estimated PR trend (0.14 % yr${}^{-\mathrm{1}}$) is opposite in sign and double the magnitude of the expected trend ($-$0.07 % yr${}^{-\mathrm{1}}$). The estimated CUP trend is positive but roughly 0, which is a little smaller than the expected trend of 0.12 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{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 $-$0.16 $\mathrm{ppm}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$.

We find that the analysis of synthetic time series indicates that ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ metrics can reproduce prescribed phase changes to within 30 %, but trends with a magnitude of $<0.1$ $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ were uncertain in magnitude and sign. Strong shifts in spring and autumn phase caused changes in PU and PR of $<\mathrm{6}$ % due to aliasing. Strong trends in PU and PR were estimated to within 25 %.

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. , where Fig.  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 $\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{2}}$, while changes in the phase were within 1 $\mathrm{day}{\mathrm{yr}}^{-\mathrm{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  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 ${\mathrm{CO}}_{\mathrm{2}}$ DZCP trend ($-0.57\pm \mathrm{4}$ %) but also a large spread about the mean bias that suggests that the ${\mathrm{CO}}_{\mathrm{2}}$ DZCP is more susceptible to aliasing. However, the use of ${\mathrm{\Delta }}_t{\mathrm{CO}}_{\mathrm{2}}=25$ % PU resulted in a relatively small mean bias ($-14\pm \mathrm{2}$ %) with high consistency ($r^{\mathrm{2}}=0.94$). Although the mean bias was less in the MCS for the ${\mathrm{CO}}_{\mathrm{2}}$ UZCP ($-\mathrm{1}\pm \mathrm{3}$ %), it was less consistent ($r^{\mathrm{2}}=0.80$). The ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ UZCP had a mean bias of $-23\pm \mathrm{1}$ % ($r^{\mathrm{2}}=0.97$). Differences between the spring and autumn phase biases calculated from ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ phase metrics carry through to the respective CUP estimates, where the ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ CUP had a mean bias of $-28\pm \mathrm{1}$ % ($r^{\mathrm{2}}=0.93$) relative to a bias of $-55\pm \mathrm{1}$ % ($r^{\mathrm{2}}=0.45$) in the ${\mathrm{CO}}_{\mathrm{2}}$ CUP. Estimates of ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ UZCP. We expect that this is a result of the asymmetry of the high-latitude ${\mathrm{CO}}_{\mathrm{2}}$ seasonal cycle. Analysis of peak rates of uptake and release resulted in mean biases of $-18\pm \mathrm{2}$ and $-28\pm \mathrm{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$ $\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$ for PU and PR, and $>0.1$ $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ for changes in phase). The ${\mathrm{CO}}_{\mathrm{2}}$ phase metric trend estimates were the most likely to have the wrong sign compared to the ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ phase metrics by a factor of 4.5, 4, and 1.5 for the DZCP, UZCP, and CUP respectively. Finally, the ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ metrics were far more effective at detecting statistically significant trends where the ${\mathrm{CO}}_{\mathrm{2}}$ metrics typically missed 33–50 % of them.

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.

[Figure omitted. See PDF]

Regression of expected and estimated linear trends for peak uptake (PU), peak release (PR), and the ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ phase metrics. Coloured points represent (red) trends that were not statistically significant, (black) trends where we estimated the incorrect sign of the trend, (blue) statistically significant trends that were successfully detected, and (green) statistically significant trends that were not detected in the analysis. Statistical significance is at the 5 % level. The numbers, $N_{\mathrm{1}\mathrm{\dots }n}$ are the number of points in each category and sum to 1000.

[Figure omitted. See PDF]

Scatterplot and associated linear regression coefficients of the amplitude trend ($\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{1}}$) against the trend in integrated ${\mathrm{CO}}_{\mathrm{2}}$ uptake and release ($\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{2}}$) from the 1000-member ensemble used in the Monte Carlo experiment.

[Figure omitted. See PDF]

Figure  shows a regression of the linear trend in integrated ${\mathrm{CO}}_{\mathrm{2}}$ uptake and release against the estimated seasonal amplitude from the individual MCS runs. We find that the linear trends in annually integrated ${\mathrm{CO}}_{\mathrm{2}}$ uptake ($\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{2}}$) are highly correlated with the amplitude trend ($\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{1}}$), but the amplitude trends are poorly correlated with changes in integrated release of ${\mathrm{CO}}_{\mathrm{2}}$. Previous work has shown that this is due to the rapid temporal change in ${\mathrm{CO}}_{\mathrm{2}}$ during the period of net carbon uptake relative to the more gradual release of ${\mathrm{CO}}_{\mathrm{2}}$ outside of this period .

Analysis of detrended CO${}_{\mathrm{2}}$ seasonal cycle

Our analysis of phase changes in the ${\mathrm{CO}}_{\mathrm{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 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ ($p<0.01$) and $-$0.18 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ ($p<0.05$) respectively. This results in a more conserved carbon uptake period, with a coefficient of 0.02 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ ($p>0.1$), which is consistent with the ecosystem having an intrinsic or fixed uptake period (not shown). Recent work using changes in ${\mathrm{CO}}_{\mathrm{2}}$ has reported a change of $-$0.17 $\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$ for the downward phase over a similar time period . Although we have included these values for completeness, we have already shown that there are significant weaknesses in using detrended ${\mathrm{CO}}_{\mathrm{2}}$ as opposed to ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ for this analysis.

Analysis of other sites

Table  summarizes the analysis of all the high northern latitude sites we have considered in this study.

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 ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ data for seven high-latitude measurement sites (Fig. ). The 95 % confidence intervals and $p$ values are calculated for each trend estimate.

Site info		Spring phase					Autumn phase					Uptake period			C. exchange
Site	Timespan	${\mathrm{CO}}_{\mathrm{2}}$ DZCP ($\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$)	${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ DZCP ($\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$)	${\mathrm{\Delta }}_t{\mathrm{CO}}_{\mathrm{2}}=25$ % PU ($\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$)	${\mathrm{\Delta }}_t{\mathrm{CO}}_{\mathrm{2}}=\text{PU}$ ($\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$)		${\mathrm{CO}}_{\mathrm{2}}$ UZCP ($\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$)	${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ UZCP ($\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$)	${\mathrm{\Delta }}_t{\mathrm{CO}}_{\mathrm{2}}=25$ % PR ($\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$)	${\mathrm{\Delta }}_t{\mathrm{CO}}_{\mathrm{2}}=\text{PR}$ ($\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$)		${\mathrm{CO}}_{\mathrm{2}}$ CUP ($\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$)	${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ CUP ($\mathrm{days}{\mathrm{yr}}^{-\mathrm{1}}$)		Seas. amp. ($\mathrm{ppm}{\mathrm{yr}}^{-\mathrm{1}}$)	PU ($\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$)	PR ($\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$)
ALT	1986–2013	$-0.14\pm 0.15$ ($p<0.1$)	$-0.34\pm 0.83$ ($p>0.1$)	$-0.16\pm 0.26$ ($p>0.1$)	$-0.18\pm 0.26$ ($p>0.1$)		$-0.27\pm 0.19$ ($p<0.01$)	$-0.10\pm 0.17$ ($p>0.1$)	$-0.09\pm 0.20$ ($p>0.1$)	$0.15\pm 0.59$ ($p>0.1$)		$-0.02\pm 0.20$ ($p>0.1$)	$0.05\pm 0.32$ ($p>0.1$)		$0.10\pm 0.04$ ($p<0.01$)	$0.61\pm 0.60$ ($p<0.01$)	$0.40\pm 0.60$ ($p<0.1$)
BRW	1973–2013	$-0.20\pm 0.08$ ($p<0.01$)	$-0.02\pm 0.47$ ($p>0.1$)	$-0.14\pm 0.14$ ($p<0.05$)	$-0.21\pm 0.15$ ($p<0.01$)		$-0.18\pm 0.14$ ($p<0.05$)	$-0.25\pm 0.08$ ($p<0.01$)	$-0.26\pm 0.10$ ($p<0.01$)	$-0.25\pm 0.10$ ($p<0.1$)		$0.02\pm 0.15$ ($p>0.1$)	$-0.11\pm 0.16$ ($p>0.1$)		$0.09\pm 0.02$ ($p<0.01$)	$0.65\pm 0.34$ ($p<0.01$)	$0.42\pm 0.34$ ($p<0.05$)
CBA	1979–2012	$-0.14\pm 0.15$ ($p<0.1$)	$-0.56\pm 0.34$ ($p<0.01$)	$0.06\pm 0.11$ ($p>0.1$)	$-0.24\pm 0.37$ ($p>0.1$)		$-0.27\pm 0.27$ ($p<0.05$)	$-0.16\pm 0.17$ ($p<0.1$)	$-0.17\pm 0.20$ ($p<0.1$)	$0.14\pm 0.33$ ($p>0.1$)		$-0.07\pm 0.29$ ($p>0.1$)	$-0.22\pm 0.34$ ($p>0.1$)		$0.07\pm 0.04$ ($p<0.01$)	$0.66\pm 0.48$ ($p<0.01$)	$0.58\pm 0.48$ ($p<0.05$)
ICE	1993–2013	$0.34\pm 0.27$ ($p<0.05$)	$0.62\pm 0.98$ ($p<0.01$)	$0.63\pm 0.65$ ($p<0.01$)	$0.25\pm 0.54$ ($p<0.1$)		$-0.13\pm 0.28$ ($p>0.1$)	$0.18\pm 0.25$ ($p>0.1$)	$0.22\pm 0.24$ ($p<0.1$)	$0.11\pm 0.99$ ($p>0.1$)		$-0.21\pm 0.33$ ($p>0.1$)	$-0.45\pm 0.64$ ($p>0.1$)		$0.06\pm 0.04$ ($p<0.01$)	$0.97\pm 0.94$ ($p<0.01$)	$0.92\pm 0.92$ ($p<0.05$)
SHM	1987–2012	$-0.40\pm 0.18$ ($p<0.01$)	$-0.59\pm 0.45$ ($p<0.05$)	$-0.45\pm 0.34$ ($p<0.05$)	$-0.54\pm 0.40$ ($p<0.01$)		$-0.27\pm 0.22$ ($p<0.05$)	$-0.13\pm 0.23$ ($p>0.1$)	$-0.11\pm 0.25$ ($p>0.1$)	$-0.15\pm 0.33$ ($p>0.1$)		$-0.13\pm 0.24$ ($p>0.1$)	$0.32\pm 0.44$ ($p>0.1$)		$0.06\pm 0.05$ ($p<0.05$)	$-0.24\pm 0.75$ ($p>0.1$)	$-0.05\pm 0.69$ ($p>0.1$)
STM	1981–2010	$-0.17\pm 0.14$ ($p<0.05$)	$-0.60\pm 0.74$ ($p>0.1$)	$-0.03\pm 0.65$ ($p>0.1$)	$-0.04\pm 0.27$ ($p>0.1$)		$-0.24\pm 0.25$ ($p<0.1$)	$-0.01\pm 0.15$ ($p>0.1$)	$-0.01\pm 0.18$ ($p>0.1$)	$0.16\pm 0.62$ ($p>0.1$)		$-0.07\pm 0.31$ ($p>0.1$)	$0.02\pm 0.66$ ($p>0.1$)		$0.05\pm 0.03$ ($p<0.01$)	$0.04\pm 0.63$ ($p>0.1$)	$0.72\pm 0.62$ ($p<0.05$)
ZEP	1994–2013	$-0.01\pm 0.21$ ($p>0.1$)	$-1.24\pm 1.78$ ($p>0.1$)	$0.01\pm 0.61$ ($p>0.1$)	$-0.06\pm 0.40$ ($p>0.1$)		$0.40\pm 0.52$ ($p>0.1$)	$-0.16\pm 0.33$ ($p>0.1$)	$-0.24\pm 0.38$ ($p>0.1$)	$0.43\pm 1.20$ ($p>0.1$)		$0.12\pm 0.25$ ($p>0.1$)	$-0.18\pm 0.76$ ($p>0.1$)		$0.14\pm 0.05$ ($p<0.01$)	$1.00\pm 1.06$ ($p<0.05$)	$-0.30\pm 1.07$ ($p>0.1$)

Analysis of ancillary data

Surface temperature and NDVI

Table  shows that mean surface land temperature has warmed significantly at high latitudes since 1970. We define a thermal growing season (TGS) with a threshold temperature of 5 ${}^{\circ }$C, the minimal temperature typically required for the onset of photosynthesis, following . We find that an earlier onset of the mean temperature reaching 5 ${}^{\circ }$C in spring, ${\text{TGS}}_{\text{BEG}}$, and a delay in the temperature dropping below 5 ${}^{\circ }$C in autumn, ${\text{TGS}}_{\text{END}}$, results in a significant lengthening of the thermal growing season, ${\text{TGS}}_{\text{LEN}}$ since 1970 for a number of high-latitude regions. Of the Transcom regions, we find that Europe exhibits the largest change in ${\text{TGS}}_{\text{LEN}}$ of $\sim 3.41\pm 0.9$ $\mathrm{days}{\mathrm{decade}}^{-\mathrm{1}}$, resulting from equal shifts in ${\text{TGS}}_{\text{BEG}}$ and ${\text{TGS}}_{\text{END}}$. Europe is followed by roughly equal changes in boreal North America and Asia, however these regions exhibit different changes in spring and autumn temperature. The largest overall changes are seen $>60$${}^{\circ }$ N where ${\text{TGS}}_{\text{LEN}}$ has increased by up to $\mathrm{5}\pm 1.7$ $\mathrm{days}{\mathrm{decade}}^{-\mathrm{1}}$ where a larger proportion of this change is due to autumn warming. This increase in ${\text{TGS}}_{\text{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$${}^{\circ }$ N) and 20 days ($>60$${}^{\circ }$ N) since 1970, consistent with previous findings . Table  shows the relationship between northern high-latitude land surface temperature anomalies with the BRW ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$ phase metrics throughout 1973–2012. We find there are significantly different results depending on whether ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{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 . 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.

${\mathit{\delta }}^{13}\mathrm{C}$ data

Figure  shows ${\mathit{\delta }}^{13}\mathrm{C}$ data over CBA, with the corresponding ${\mathrm{CO}}_{\mathrm{2}}$ mole fraction data. Measurements of ${\mathit{\delta }}^{13}\mathrm{C}$ show a strong seasonal variation, which is anti-correlated with ${\mathrm{CO}}_{\mathrm{2}}$. Plants preferentially take the lighter carbon ${}^{12}\mathrm{C}$ isotope out of the atmosphere through photosynthesis during spring and summer resulting in an increase in ${\mathit{\delta }}^{13}\mathrm{C}$ and release more ${}^{12}\mathrm{C}$ than ${}^{13}\mathrm{C}$ during autumn and winter resulting in a decrease in ${\mathit{\delta }}^{13}\mathrm{C}$.

Figure  shows a similar phase analysis for $(-\mathrm{1})\times {\mathit{\delta }}^{13}\mathrm{C}$ and $(-\mathrm{1})\times \mathrm{\Delta }{\mathit{\delta }}^{13}\mathrm{C}$, comparing it with variability and trends with the corresponding ${\mathrm{CO}}_{\mathrm{2}}$ values. Table  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 ${\mathrm{CO}}_{\mathrm{2}}$ DZCP and UZCP can be explained by variations in colocated measurements of ${\mathit{\delta }}^{13}\mathrm{C}$. This suggests that the terrestrial biosphere is largely responsible for observed ${\mathrm{CO}}_{\mathrm{2}}$ variability with the remainder due to atmospheric transport and other minor source variations. This result is consistent with previous work that showed using an atmospheric transport model that atmospheric transport variations contributed $<\mathrm{7}$ % of the observed variation in ${\mathrm{CO}}_{\mathrm{2}}$ seasonal amplitudes at high northern latitudes.

Temperature linear trend analysis (1970–2011) for the beginning and end of the thermal growing season.

	${\text{TGS}}_{\text{BEG}}$ ($\mathrm{days}{\mathrm{decade}}^{-\mathrm{1}}$)					Spring $T$ (${}^{\circ }\mathrm{C}{\mathrm{decade}}^{-\mathrm{1}}$)
Region	Trend	unc	$r^{\mathrm{2}}$	$p$ value		Trend	unc	$r^{\mathrm{2}}$	$p$ value
ASBor	$-$1.39	$\pm$0.49	0.45	$<0.01$		0.58	0.27	0.32	$p<0.01$
Europe	$-$1.67	$\pm$0.52	0.51	$<0.01$		0.33	0.11	0.47	$p<0.01$
USBor	$-$1.06	$\pm$0.72	0.18	$<0.01$		0.34	0.25	0.15	$p<0.05$
$>45$${}^{\circ }$ N	$-$1.24	$\pm$0.44	0.44	$<0.01$		0.41	0.13	0.49	$p<0.01$
$>60$${}^{\circ }$ N	$-$2.12	$\pm$0.75	0.45	$<0.01$		0.45	0.18	0.40	$p<0.01$
	${\text{TGS}}_{\text{END}}$					Autumn $T$
Region	Trend	unc	$r^{\mathrm{2}}$	$p$ value
ASBor	1.07	$\pm$0.79	0.16	$<0.01$		0.57	0.29	0.28	$p<0.01$
Europe	1.74	$\pm$0.66	0.42	$<0.01$		0.38	0.12	0.49	$p<0.01$
USBor	1.57	$\pm$0.69	0.35	$<0.01$		0.47	0.21	0.34	$p<0.01$
$>45$${}^{\circ }$ N	1.34	$\pm$0.47	0.45	$<0.01$		0.44	0.13	0.55	$p<0.01$
$>60$${}^{\circ }$ N	2.85	$\pm$1.04	0.43	$<0.01$		0.52	0.16	0.53	$p<0.01$
	${\text{TGS}}_{\text{LEN}}$					Annual $T$
Region	Trend	unc	$r^{\mathrm{2}}$	$p$ value
ASBor	2.46	$\pm$1.08	0.35	$<0.01$		0.45	0.16	0.45	$p<0.01$
Europe	3.41	$\pm$0.90	0.60	$<0.01$		0.35	0.10	0.61	$p<0.01$
USBor	2.63	$\pm$1.25	0.31	$<0.01$		0.43	0.16	0.41	$p<0.01$
$>45$${}^{\circ }$ N	2.57	$\pm$0.78	0.52	$<0.01$		0.40	0.10	0.67	$p<0.01$
$>60$${}^{\circ }$ N	4.97	$\pm$1.69	0.47	$<0.01$		0.43	0.11	0.63	$p<0.01$

Linear regression coefficients that describe the relationship between changes in ${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$, and temperature phase metrics at different latitude bands in the high northern latitudes (1973–2012).

${\mathrm{CO}}_{\mathrm{2}}$
	DZCP vs. ${\text{TGS}}_{\text{BEG}}$	UZCP vs. ${\text{TGS}}_{\text{END}}$	CUP vs. ${\text{TGS}}_{\text{LEN}}$
$>45$${}^{\circ }$ N	$1.02\pm 0.47$	$-0.13\pm 0.73$	$0.31\pm 0.44$
	($r^{\mathrm{2}}=0.34$, $p<0.01$)	($r^{\mathrm{2}}=0.01$, $p>0.1$)	($r^{\mathrm{2}}=0.04$, $p>0.1$)
$>60$${}^{\circ }$ N	$0.63\pm 0.27$	$-0.11\pm 0.34$	$0.14\pm 0.22$
	($r^{\mathrm{2}}=0.38$, $p<0.01$)	($r^{\mathrm{2}}=0.01$, $p>0.1$)	($r^{\mathrm{2}}=0.04$, $p>0.1$)
${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$
	DZCP vs. ${\text{TGS}}_{\text{BEG}}$	UZCP vs. ${\text{TGS}}_{\text{END}}$	CUP vs. ${\text{TGS}}_{\text{LEN}}$
$>45$${}^{\circ }$ N	$0.57\pm 0.81$	$-0.89\pm 0.51$	$-0.42\pm 0.48$
	($r^{\mathrm{2}}=0.05$, $p>0.1$)	($r^{\mathrm{2}}=0.25$, $p<0.01$)	($r^{\mathrm{2}}=0.08$, $p<0.1$)
$>60$${}^{\circ }$ N	$0.28\pm 0.48$	$-0.42\pm 0.24$	$-0.19\pm 0.24$
	(r${}^{\mathrm{2}}$=0.04, $p>0.1$)	($r^{\mathrm{2}}=0.26$, $p<0.01$)	($r^{\mathrm{2}}=0.07$, $p<0.1$)

Linear regression coefficients that describe the relationship between changes in ${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$, and ${\mathit{\delta }}^{13}\mathrm{C}$ at BRW during the overlapping time span of the data (1990–2012).

	DZCP	UZCP	CUP
${\mathrm{CO}}_{\mathrm{2}}$	$0.94\pm 0.19$	$0.59\pm 0.34$	$0.60\pm 0.37$
vs. ${\mathit{\delta }}^{13}\mathrm{C}$	($r^{\mathrm{2}}=0.84$, $p<0.01$)	($r^{\mathrm{2}}=0.38$, $p<0.01$)	($r^{\mathrm{2}}=0.34$, $p<0.01$)
${\mathrm{\Delta }}_t$${\mathrm{CO}}_{\mathrm{2}}$	$0.68\pm 0.37$	$0.95\pm 0.29$	$0.88\pm 0.42$
vs. $\mathrm{\Delta }{\mathit{\delta }}^{13}$C	($r^{\mathrm{2}}=0.41$, $p<0.01$)	($r^{\mathrm{2}}=0.70$, $p<0.01$)	($r^{\mathrm{2}}=0.48$, $p<0.01$)

Acknowledgements

We thank NOAA/ESRL for the ${\mathrm{CO}}_{\mathrm{2}}$ surface mole fraction data which is provided by NOAA/ESRL PSD, Boulder, Colorado, USA, from their website http://www.esrl.noaa.gov/psd/. We would also like to thank for making the wavelet transform code freely available at the website http://paos.colorado.edu/research/wavelets/software.html. J. M. Barlow acknowledges the centre for Earth Observation Instrumentation and the National Environmental Research Council for funding his studentship, number NE/1528818/1. P. I. Palmer thanks Donald Percival (U. Washington, Seattle) for a useful discussion and acknowledges support from his Philip Leverhulme Prize and his Royal Society Wolfson Research Merit Award. Edited by: P. Monks