Precipitation and carbon-water coupling jointly control the interannual variability of global land gross primary production

Carbon uptake by terrestrial ecosystems is increasing along with the rising of atmospheric CO2 concentration. Embedded in this trend, recent studies suggested that the interannual variability (IAV) of global carbon fluxes may be dominated by semi-arid ecosystems, but the underlying mechanisms of this high variability in these specific regions are not well known. Here we derive an ensemble of gross primary production (GPP) estimates using the average of three data-driven models and eleven process-based models. These models are weighted by their spatial representativeness of the satellite-based solar-induced chlorophyll fluorescence (SIF). We then use this weighted GPP ensemble to investigate the GPP variability for different aridity regimes. We show that semi-arid regions contribute to 57% of the detrended IAV of global GPP. Moreover, in regions with higher GPP variability, GPP fluctuations are mostly controlled by precipitation and strongly coupled with evapotranspiration (ET). This higher GPP IAV in semi-arid regions is co-limited by supply (precipitation)-induced ET variability and GPP-ET coupling strength. Our results demonstrate the importance of semi-arid regions to the global terrestrial carbon cycle and posit that there will be larger GPP and ET variations in the future with changes in precipitation patterns and dryland expansion.


S1. GPP products  VPM
The Vegetation Photosynthesis Model (VPM) is a light use efficiency (LUE) model which has been tested over a variety of land cover types at CO2 eddy flux tower sites (1)(2)(3)(4)(5)(6). The GPP product we used in this paper is driven by the MODIS 8-day 500 m surface reflectance product (MOD09A1 C6), the MODIS 8-day 1 km Land surface temperature product (MOD11A2 C5), the MODIS land cover product (MCD12A2 C5), and the shortwave radiation and air temperature from the NCEP-DOE Reanalysis 2 dataset. The VPM (version 2.0) estimates GPP by using the product of absorbed photosynthetic active radiation by the chlorophyll ( , calculated by PAR) and the light use efficiency (LUE, ): Where the factor = 1.25, and 0.1 in Eq. A2 is used to adjust EVI baseline (7,8). The EVI is calculated from the MOD09A1 product and passed a rigorous quality check and gap-filling procedure (9). The LUE ( ) is reduced from a maximum light use efficiency ( ) by the temperature scalar ( ) and a water scalar ( ):

A3
Detailed information about the model parameter estimation were documented in previous publications (1,4,10).
The GPP product at 8-day temporal resolution and 500 m spatial resolution was reprojected and aggregated into 0.5 degree (latitude and longitude). The annual GPP is calculated as the sum of each 8-day GPP over a Julian year, and the interannual variation of GPP is calculated below:

∑ A4
Where and represent the annual GPP for year and average annual GPP for year 2000 to 2011, respectively.
is the total number of years.

 MPI-BGC (MTE)
The MPI-BGC estimates GPP by upscaling global in situ CO2 eddy flux observations with climate data and remote sensing fraction of absorbed photosynthetic active radiation (fAPAR) using a Model Tree Ensemble approach (11,12). This dataset has a spatial resolution of 0.5° 0.5°. This monthly dataset was first aggregated into annual sum and then the IAV over the period of 2000 to 2011 was calculated to match with other GPP products.

 MOD17
The MODIS GPP product (MOD17) employs a light use efficiency (LUE) approach to calculate GPP. The PSN model used in the MOD17 GPP product uses MODIS fPAR product as the fraction of PAR absorbed by vegetation for photosynthesis. MOD17 also uses vapor pressure deficit (VPD) as the water limitation of LUE and uses variable maximum LUE for individual biome types, which are determined by MODIS land cover map (MCD12). The MOD17 product (MOD17A2 C55) during 2000-2011 was downloaded from Numerical Terradynamic Simulation Group (http://www.ntsg.umt.edu/project/mod17#data-product) at the University of Montana. MOD17 GPP at 1-km spatial resolution and 8-day temporal resolution and was aggregated into 0.5 degree (latitude/longitude) to calculate annual GPP.

S2. Evapotranspiration and potential evapotranspiration dataset
Evapotranspiration datasets from the model tree ensemble (MTE) (Jung et al. (12)), MODIS ET (27) (MOD16A3 C5), and 10 TREDNY models are used in our study (ET from LPJ is not used due to a data problem---abnormally low ET in tropical reigons). The MTE method integrates the eddy flux tower measured ET with remote sensing satellite images and meteorological data using a machine-learning algorithm (11). The MOD16A3 product calculates ET based on the Penman-Monteith algorithm (28) and was recently improved (27). We downloaded the potential evapotranspiration (PET) and evapotranspiration (ET) from Numerical Terradynamic Simulation Group website at the University of Montana (http://www.ntsg.umt.edu/project/mod16). The PET and ET datasets were annual sums at 1 km 1 km spatial resolution and we aggregated to 0.5° 0.5°. ET from TRENDY-V4 were also interpolated into 0.5° 0.5° spatial resolution. All ET datasets from 2000 to 2011 were used.

S3. Climate dataset
The GPCC V7 global precipitation at 0.5° 0.5° is used in our analysis (29). This dataset has a monthly precipitation value from 1901 to 2013. The downward shortwave radiation was obtained from CRUNCEP V4 (http://dods.extra.cea.fr/data/p529viov/cruncep/) with a monthly sum and a spatial resolution of 0.5° 0.5°. The monthly mean temperature was from CRU TS 3.23 (30) (https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_3.23/) and also has a 0.5° 0.5° spatial resolution. All the monthly climate datasets are aggregated to annual sums. The aridity index is calculated using GPCC precipitation over PET from MOD16 for the period 2000 to 2011. The classification scheme of the aridity regions based on aridity index is from UNEP (31) and can be found in Table S1. Because CRUNCEP V4 does not provide radiation after 2010, the climate correlation is calculated using data from 2000 to 2010.

S4. Aggregation of MODIS land cover map
In order to get the annual land cover map at 0.5° 0.5° spatial resolution, we aggregate the MCD12C1 C5 product spatially. The aggregation is based on the average of area percentage of each biome type, therefore the output also includes percentage of each biome type for each grid cell. To make sure we can compare the SIF and GPP within each land cover type, we set a threshold to extract the grid cells with 'pure' land cover types. We set the threshold from 50% to 100%. With the increase of this threshold, the number of mixed pixels also increased and the available pixels decreased (Fig. S4). In this study, we used the 80% threshold, and only two biome types have less than 100 gridcells.

S5. Rationale of GPP ensemble using SIF as a reference
Recent studies have successfully retrieved SIF from satellites, and these SIF retrievals showed very good correlation with GPP from flux tower upscaled models and light use efficiency models (32). This close relationship between GPP and SIF is found to be higher within each biome type at monthly scale using data from both global simulation and in situ observations (33,34). The Soil Canopy Observation of Photochemistry and Energy fluxes (SCOPE) model is a process-based model which can also simulate energy exchange within photosynthesis process, including SIF (35). Using this SCOPE model, recent studies investigated the relationship between GPP and SIF (36,37). SIFyield (defined as SIF/APAR, absorbed photosynthetically active radiation) shows a good linear correlation with photochemical yield (defined as GPP/APAR) when APAR is higher than 400 to 600 μmol m s . However, this linear relationship varies with Vcmax. If we assume the similar Vcmax within each biome type, SIF should have a good correlation with GPP at moderate or high light conditions relevant to satellite observations.

S6. Comparing GPP ensemble using SIF as a temporal proxy and spatial proxy
Because of the close relationship between SIF and GPP within each biome type, SIF can be used as a reference to evaluate the performance of various models in GPP estimation. To reduce uncertainties of SIF from GOME-2, we proposed two strategies to reduce uncertainties: monthly SIF is averaged either temporally to yearly average, or spatially within each biome types; and we then use temporally averaged SIF as spatial representative or the spatially averaged SIF as temporal representative of photosynthetic activity. The weights for each model from these two approaches are given in Figs. S5, S6. When using SIF as a temporal reference for each biome, the difference between each individual model is rather small (Table S2), and most models have relatively high correlation. This implies most models can simulate the seasonal variation of GPP relatively well. However, when using SIF as a spatial reference within each biome, the difference between each individual model is relatively large (Table S3), and many models have relatively low correlation coefficients. In addition, the interannual difference of the weight is very small (short error bar in Fig. S5), proving the method is relatively stable.

S7. Interannual variation of SIF
SIF data from the satellite have relative high noise because of the high spectral requirement for the sensor. For the GOME-2 dataset, the reported error for monthly mean gridded data are ~ 0.1-0.4 mW m -2 nm -1 sr -1 (26). This reported error is contributed from several sources: (1) the error for each individual SIF measurement, which is affected by the radiance noise, atmospheric conditions; (2) different observing conditions (view angle), sampling locations and time; (3) cloud contaminations; (4) various systematic errors. When calculating the interannual variation, the contribution from 2-4 is limited and is ignored in our analysis.
We use the following scheme to calculate the interannual variation of SIF: 1) For each gridcell each month, not only the aggregated SIF value (average of all valid SIF retrieval) is given, the GOME-2 SIF product also gives the number of individual SIF measurements used to aggregate the gridcell and the standard deviation of all retrievals. The uncertainty of this gridcell for month (σ ) will be the standard deviation of all retrievals (regarded as the uncertainty of each individual measurement) divided by √ , where is the number of measurements used for aggregation. Therefore, we can get the uncertainty of SIF for each gridcell for each month. 2) From the SIF value and uncertainty for each month, we will further calculate the SIF mean and uncertainties for each year. Assuming for each gridcell, the uncertainty of SIF for month (ε ) follows a normal distribution (ε ∈ N 0, σ ), where σ is the uncertainty which we calculated in the previous step, the uncertainty of annual average SIF for year (ε ) will also follow a normal distribution: ε ∈ N 0, ∑ σ 12 5 3) To calculate the actual IAV of SIF, we need to eliminate the uncertainty induced variability for each year (ε ∈ N 0, σ ) from the calculated variability ( ) for each gridcell. We found a very similar pattern of ε across different years. Because of the additivity of variance, the calculated IAV of SIF ( ) can be explained by the real IAV of SIF ( ) and error of observation (σ ): σ A6 where the error of observation (σ ) is calculated as follows: The spatial patterns of and σ are shown in Fig. S10.

S8. Decompose ET variation from precipitation variation and potential evapotranspiration variation
Using the Budyko framework, similar with equation (7) in the methods, we can also get the ET sensitivity to PET (E0) change: Therefore, ET variability can be decomposed into ET variability caused by precipitation variability (ET ), and ET variability caused by PET variability (ET ).

ET A10 ET A11
We calculate the ET and ET using different n values (0.5, 1, 2) to infer the ET change from change in precipitation and change in E0. Results are shown in Fig. S15.

S9. Verification of the stability and robustness of the weighted ensemble method
To verify the stability and robustness of the weighted ensemble method to model inputs, we use 13 out of 14 models as input and test the variation of the annual GPP and GPP s.d. for both unweighted average method and the weighted ensemble method. Each time, we drop one model from the 14 models and use the rest 13 models as input. The GPP for each year is calculated as the average of the 13 model estimates for this year or using the ensemble method described in the method section. The annual GPP and GPP s.d. is then calculated from the output of each year by both methods. By repeating this process for 14 times, we can get all the possible combination of using 13 models as inputs. The variations of the annual GPP and GPP s.d. from 14 times run for both methods are regarded as the methods stability (Fig. S17).   Process-based models Data-driven models GOME-2 SIF Fig. S3. Spatial patterns of inter-annual variability (standard deviation, s.d.) of GPP from the 3 data-driven models and the 11 DGVMs from the Trendy v4 project. Maps were created using Matlab R2016a (http://www.mathworks.com/products/matlab/).             3). All the units are g C m -2 yr -1 . Maps were created using Matlab R2016a (http://www.mathworks.com/products/matlab/). Table S1. Classification scheme of aridity regions based on aridity index (AI). Aridity index is defined by FAO and the classification is inherited from UNEP Drylands are defined as AI < 0.65.