Observational LAI product (LAI3gV1)

The new version (V1) of the LAI data set is an update of the widely used LAI3g data set26. It was generated using an artificial neural network (ANN) and the latest version (third generation) of the Global Inventory Modeling and Mapping Studies group (GIMMS) Advanced Very High Resolution Radiometer (AVHRR) normalized difference vegetation index (NDVI) data (NDVI3g). The latter has been corrected for sensor degradation, intersensor differences, cloud cover, solar zenith angle, viewing angle effects due to satellite drift, Rayleigh scattering, and stratospheric volcanic aerosols36. The ANN model was trained with overlapping data of NDVI3g and Collection 6 Terra MODIS LAI product37,38, and then applied to the full NDVI3g time series to generate the LAI3gV1 data set. This data set provides global and year-round LAI observations at 15-day (bi-monthly) temporal resolution and 1/12 degree spatial resolution from July 1981 to December 2016. Currently, it is the only data set that spans this long period.

The quality of the previous version (V0) of the GIMMS LAI3g data set was evaluated through direct comparisons with ground-based measurements of LAI, indirectly with other estimates from similar satellite-data products, and also through statistical analysis with climatic variables, such as temperature and precipitation variability26. The LAI3gV0 data set (and related fraction vegetation-absorbed photosynthetically active radiation data set) has been widely used in various studies5,15,19,20,31,39,40,41. The new version LAI3gV1 used in our study is an update of that earlier version.

For both, observational and CMIP5 data, LAI is defined as the one-sided green leaf area per unit ground area in broadleaf canopies and as one-half the green needle surface area in needleleaf canopies. It is expressed in units of m2 green leaf area per m2 ground area. In this study, we use the annual maximum value of LAI, LAI max , to quantify the greenness level of a surface. LAI max is less influenced by cloudiness and noise; accordingly, it is most useful in investigations of long-term greening and browning trends. The drawback of LAI max is the saturation effect at high LAI values42. However, this is less of a problem in high latitudinal ecosystems which are mostly sparsely vegetated, with LAI max values typically in the range of 2–3.

The bi-monthly GIMMS LAI3gV1 data are merged to a monthly temporal resolution by averaging the two composites in the same month. Then, for model and observational data alike, the two-dimensional global fields are cropped to the northern high latitudinal band defined as 60°N to 90°N, averaged in space and temporally reduced to the annual maximum value.

Although the AVHRR data underlying the LAI data in this study have corrections for various deleterious effects36, the data may still contain residual non-vegetation-related effects. Therefore, we sought confirmation of the greening trend19, on which the current study relies, from a more reliable but shorter record from the MODIS sensors37,38. These data are well calibrated, cloud-screened, and corrected for atmospheric effects, especially tropospheric aerosols. The sensor-platforms are regularly adjusted to maintain precise orbits. All algorithms, including the LAI algorithm, are physics-based, well-tested and currently producing the sixth generation data sets. The results, not shown here for brevity, illustrate global scale greening, across all latitudinal zones and broad vegetation classes. Zhang et al.43 also reported matching greening trends between the latest (Version 6) MODIS and AVHRR (Version 3) vegetation index data sets.

We also found that the LAI max sensitivity derived with MODIS LAI data matched well with that obtained from the AVHRR LAI data (results not shown for brevity). Whether this indicates that the 17-year MODIS record from the period 2000 to 2016 captures information similar to the longer AVHRR record (1981–2016), or is simply a fortuitous occurrence, is not known, and deserves further study. In the present context, however, this adds confidence to the AVHRR LAI data used in our study.

Temperature data from ECMWF ERA-Interim

Estimates of surface air temperature at 2 m height are from the widely used global atmospheric reanalysis product ERA-Interim by ECMWF44 (for details see https://www.ecmwf.int/en/research/climate-reanalysis/era-interim). The global temperature fields were retrieved at a resolution of 0.5° × 0.5° for monthly mean estimates derived from daily means. Other reanalysis products with similar specifications (NCEP/NCAR reanalysis, University of Delaware Air Temperature & Precipitation, and GHCN/CAMS reanalysis product) were also investigated. The differences among the various products were found to be minor.

CMIP5 models used in this study

In this study, we analyze a set of the most recent climate-carbon simulations of seven ESMs participating in the fifth phase of the Coupled Model Intercomparison Project, CMIP528. The model data were obtained from the Earth System Grid Federation, ESGF (https://esgf-data.dkrz.de/projects/esgf-dkrz/). Seven ESMs provided output for the variables of interest for simulations esmHistorical, 1pctCO2, esmFixClim, and esmFdbk.

The esmHistorical simulation spans the period 1850–2005 and was driven by observed conditions such as solar forcing, emissions or concentrations of short-lived species and natural and anthropogenic aerosols or their precursors, land use, anthropogenic as well as volcanic influences on atmospheric composition. The models are forced by prescribed anthropogenic CO 2 emissions, rather than atmospheric CO 2 concentrations.

1pctCO2 is an idealized fully coupled carbon/climate simulation initialized from steady state of the pre-industrial control run and atmospheric CO 2 concentration prescribed to increase 1% yr−1 until quadrupling of the pre-industrial level. The simulations esmFixClim and esmFdbk and are set up as the 1pctCO 2 with the difference, that in esmFixClim (esmFdbk) only the radiative effect from increasing CO 2 concentration is included, while the carbon cycle sees the pre-industrial CO 2 level (vice versa)28,45.

Historical simulation with MPI-ESM higher-resolution setup

MPI-ESM-HR is the coupled high-resolution setup of the latest version of the Max-Planck-Institute Earth System Model MPI-ESM1.2, which is the baseline for the upcoming Coupled Model Intercomparison Project Phase 6 (CMIP6). Here, the atmospheric component ECHAM6.3 has 95 vertical levels and twice the horizontal resolution (~100 km) than the CMIP5 version. The ocean component MPIOM is set up on a tripolar grid at nominal 0.4° horizontal resolution (TP04) and 40 vertical levels. MPI-ESM1.2 includes the latest versions of the land and ocean carbon cycle modules, comprising the ocean biogeochemistry model HAMOCC and the land surface scheme JSBACH. The forcing components for the historical simulation are chosen from CMIP5 (Methods) as at the time the simulations were conducted CMIP6 forcing was not available46.

Atmospheric CO 2 concentration data

Monthly means of atmospheric CO 2 concentration at Point Barrow (71.3°N, 203.4°E) and Alert Nunavut (82.5°N, 297.7°E) are taken from the Global Monitoring Division measurement datasets (co2_brw_surface-insitu_1_ccgg_MonthlyData respectively co2_alt_surface-flask_1_ccgg_month) provided by the National Oceanic and Atmospheric Administration/Earth System Research Laboratory (NOAA/ESRL). Global monthly means of atmospheric CO 2 concentration are taken from the GLOBALVIEW-CO2 product (obspack_co2_1_GLOBALVIEWplus_v2.1_2016_09_02; for details see https://doi.org/10.15138/G3259Z) also available at NOAA/ESRL.

Atmospheric CO 2 inversion products

Atmospheric CO 2 inversions estimate surface–atmosphere net carbon exchange fluxes by utilizing CO 2 concentration measurements, a transport model and prior information on anthropogenic carbon emissions as well as carbon exchange between atmosphere and land respectively ocean47. We choose two products, which cover the longest time period (1980–2015) and are regularly updated, the Jena CarboScope32 (JENA, version s81_v3.8, for details see http://www.bgc-jena.mpg.de/CarboScope/s/s81_v3.8.html) and the Copernicus Atmosphere Monitoring Service33 (CAMS, version v15r2, for details see http://atmosphere.copernicus.eu/documentation-supplementary-products#greengas-fluxes) inversion systems. Both products provide monthly mean net flux estimates on a spatial resolution of 3.75° latitude and 5° longitude (JENA) and 1.875° latitude and 3.75° longitude (CAMS).

Calculation of growing degree days above 0 °C (GDD0)

The global temperature fields from the reanalysis and model data are cropped to the northern high latitudinal band and averaged in space. The resulting one-dimensional time-series is converted to GDD above 0 °C by multiplying the days of each month with the respective monthly mean estimate if it is above 0 °C. Thus, we not only capture the warming signal, but also the prolongation of the growing season.

Dimension reduction using principal component analysis

The drivers GDD0 and atmospheric CO 2 concentration vary co-linearly due to the radiative effect of increasing CO 2 concentration in the NHL. Thus, it is problematic to conduct an accurate factor separation in terms of their respective contribution to increase in LAI max . However, the co-linearity suggests that a large amount of the signal is shared. Therefore, we conduct a PCA to apply dimension reduction48.

The aim of the PCA is to find a linear combination of the driver variables that represents the one-dimensional projection with the largest possible variance. First, each driver time series x i is normalized by centering on its mean (subtracting \(\bar x_i\)) and scaling to unit variance (divide by one standard deviation σ i ). Thus,

$${\mathbf{X}} = x_i^\prime = \frac{{x_i - \bar x_i}}{{\sigma _i}}.$$ (1)

The matrix X contains the scaled time series \(x_i^\prime\) as columns. Next, we compute the covariance matrix CX,

$${\mathbf{C}}^{\mathrm{X}} = \frac{1}{n}{\mathbf{X}}^{\mathbf{T}}{\mathbf{X}}$$ (2)

where n is the length of each time series. The eigenvector u k is obtained by solving the eigenvalue problem,

$${\mathbf{C}}^{\mathrm{X}}{\mathbf{u}}_k = \lambda _k{\mathbf{u}}_k.$$ (3)

The eigenvectors u k are sorted according to the ordering of their associated eigenvalues λ k . Projecting the initial driver matrix X onto the eigenvector u 1 with the highest associated eigenvalue we arrive at the one-dimensional vector, the first principal component (PC),

$${\omega } = \left( {{\mathbf{Xu}}_1} \right)^{\mathbf{T}}.$$ (4)

Transposed to a row vector, ω denotes the time-series of the first PC, which explains the maximum variance of the two driver time series, atmospheric CO 2 concentration and GDD0.

Estimation of historical LAI max sensitivity

We derive the historical LAI max sensitivity applying a standard linear regression model (f n )

$$f_n = a + bx_n$$ (5)

where x n denotes the driver time series, a the intercept and b the gradient. We obtain the best-fit line by minimizing the squared error (s2)

$$s^2 = \frac{1}{{N - 2}}\mathop {\sum }\limits_{n = 1}^N (y_n - f_n)^2$$ (6)

where y n is the predictand time series and N is the number of data points of each time series. The resulting best-fit gradient b′ represents the sought sensitivity. The standard error of b and a are given by

$$\sigma _b = \frac{s}{{\sigma _x\sqrt N }}$$ (7)

and

$$\sigma _a = s\sqrt {\frac{1}{N} + \frac{{\bar x^2}}{{\sigma _x^2N}}}$$ (8)

respectively, for σ x being the standard deviation and \(\bar x\) being the mean value of x n .

Derivation of changes in NHL CO 2 drawdown slope

Graven et al.17 showed that NHL CO 2 drawdown mostly happens in June and July. ESMs, however, disagree on the phase, mainly due to a premature start of the growing season (Fig. 3a, b). As a consequence, the CO 2 drawdown in models peaks earlier in the season. To obtain comparability for changes in CO 2 drawdown strength, we calculate the first derivative of the CO 2 concentration time series for the observational sites and each model individually. The annual minimum of the derivative in each time series reflects the months where the increase in photosynthetic CO 2 fixation is strongest (CO 2 drawdown slope). This procedure does not require a detrending of the atmospheric CO 2 signal.

For the BRW record, the 30 years of continuous overlap with the CMIP5 historical simulations were used to calculate the drawdown slopes (1974–2005). Due to the shorter overlap in the ALT record, 30 years of data from 1985 (start of measuring campaign) to 2015 were used for comparison with models. This is legitimate, because the CO 2 concentration rate of increment for both periods are just about the same. Model time series are obtained from the near-surface CO 2 concentration using the grid box in close proximity to each site. All yearly time series are slightly smoothed with a 2-year moving window. Changes are calculated from 5-year averages at the beginning and end of the record. Here, we only present a low-end, high-end and the closest-to-observation model from the greening EC analysis, because Wenzel et al.23 already reported the behavior of the entire CMIP5 ensemble in terms of simulating the NHL CO 2 seasonal cycle.

Scaling of NPP estimates

We scale and convert the EC estimate for changes in the GPP flux ΔF GPP,EC for a doubling of the pre-industrial CO 2 level ([CO 2 ] pi ) to a NPP flux (ΔF NPP,EC ) to obtain a comparable estimate to the atmospheric CO 2 inversion datasets using

$${\mathrm{\Delta }}F_{{\mathrm{NPP}},{\mathrm{EC}}} = b\frac{{{\mathrm{\Delta }}[{\mathrm{CO}}_2]_{1980 - 2015}}}{{[{\mathrm{CO}}_2]_{{\mathrm{pi}}}}} \times {\mathrm{\Delta }}F_{{\mathrm{GPP}},{\mathrm{EC}}}$$ (9)

where Δ[CO 2 ] 1980–2015 denotes the change in atmospheric CO 2 concentration over the observational period from 1980 to 2015 and b the standard GPP to NPP conversion factor of 0.5 (assuming uncertainty of 10%)49,50.

Comparison of C fluxes between Arctic Ocean and NHL land

We require the use of a fully coupled ESM to separate between land and ocean in terms of the sign, magnitude, and seasonal cycle of the respective net carbon exchange fluxes with the atmosphere. We have access to a spatially-high resolved historical run (10 realizations) of the MPI-ESM which has the ability to reproduce seasonality in the Arctic Ocean (Methods). The terrestrial carbon pools have not been brought into equilibrium due to computational limitations in these high-resolution simulations. Thus, we use simulations from the same model but at low spatial resolution (3 realizations), the CMIP5 esmHistorical simulation, to address land carbon exchange fluxes (Methods).

The NHL land sink is approximately 2.5 times larger than the Arctic Ocean sink, on an annual basis. However, in terms of the change in carbon sink between the mid-1970s and early-2000s, the increase in CO 2 uptake by the land is about 15 times larger than the ocean. Accordingly, the Arctic Ocean can be ignored when trying to explain changes during the recent past, i.e., BRW period of CO 2 concentration measurements.

During the months from May to September (may-to-sep) when photosynthetic CO 2 drawdown is happening, the change in land sink is about 0.4 Pg C on an annual basis. Especially between May and July, the CO 2 concentration is rapidly declining, i.e., photosynthesis prevails CO 2 release processes. Thus, nearly the entire increase of 0.4 Pg C can be attributed to increasing NPP. The EC analysis shows that the MPI-ESM model is rather close to observations but generally underestimating greening sensitivity and thus also the GPP enhancement. These results are not provided as further proof of the EC estimate, although they are not contradictory, they are provided to compare the strength of NHL land and Arctic Ocean carbon sinks and why the ocean component can be neglected.

Bootstrapping for probability estimation

We apply bootstrapping to estimate the 68% confidence of the emergent linear relationship due to the small sample size of the CMIP5 ensemble. First, we randomly resample the data with replacement, where the size of the resample is equal to the size of the original sample N. Second, we compute the least-squares linear best fit for the resampled data. Third, we repeat this procedure m times (minimum m = 100) until the difference between the median best fit line of m − 1 and m computed regressions approaches zero (the actual threshold was set to a difference less than 1%). We derive the 68% confidence contours of equal probability based on the set of m random regression lines.

Calculation of probability density functions

We derive a probability density function (PDF) for the observed sensitivity \(b^\prime\) (associated standard error σ b , Methods)

$$P(b) = \frac{1}{{\sqrt {2{\mathrm{\pi }}\sigma _b^2} }}{\mathrm{exp}}\left\{ { - \frac{{(b - b^\prime )^2}}{{2\sigma _b^2}}} \right\}$$ (10)

assuming Gaussian distribution. The PDF of y for given x,

$$P(y|x) = \frac{1}{{\sqrt {2{\mathrm{\pi }}\sigma _f^2} }}\exp \left\{ { - \frac{{(y - f(x))^2}}{{2\sigma _f^2}}} \right\},$$ (11)

represent the contours of equal probability density around the best-fit linear regression, where σ f denotes the 68% confidence contours estimated by bootstrapping (Methods). As shown in Cox et al.21, for a given observation-based PDF P(b) and a model-based PDF P(y|x), the PDF of the EC on y is

$$P\left( y \right) = \mathop {\int }

olimits_{ - \infty }^\infty P\left\{ {x{\mathrm{|}}y} \right\}P\left( x \right)\,{\mathrm{d}}x.$$ (12)

The PDF of the CMIP5 unweighted multi-model mean is configured assuming Gaussian distribution.

Code availability

The code used in this study is available from the corresponding author upon request.