a, Pixel-by-pixel linear trends in annual NDVI. b, Areas of water-limited vegetation, determined as pixels with significant (P ≤ 0.10) positive annual NDVI–precipitation correlations. Nonsignificant or negative correlations were masked out from b. Farmlands, irrigated areas and wetlands have been masked out from both panels.
Core data sets.
Normalized difference vegetation index. We obtained a time series of third-generation NDVI (NDVI3g) from the Global Inventory Modelling and Mapping Studies (GIMMS; ref. 24). This data set is gridded at 0.083° spatial resolution and was averaged from biweekly to annual time steps. The annual average for a given grid cell was determined only if >80% of biweekly values were available and was set to missing otherwise. Similarly, pixel trends were calculated only for pixels with annual time series >80% complete. Basin-specific NDVI values were obtained by averaging gridded data over basin areas.
Climatic variables. Monthly climatic fields (precipitation, minimum and maximum air temperature and shortwave radiation) were obtained from the ANUCLIM archive25. The Australia-wide data are gridded at 0.05° spatial resolution and were produced by the ANUSPLIN software package25, 26 from meteorological station data using a thin-plate smoothing spline.
An annual time series of atmospheric CO2 concentrations was obtained from National Oceanic and Atmospheric Administration Earth System Research Laboratory (NOAA ESRL; http://www.esrl.noaa.gov/gmd/ccgg/trends). The data report the mean annual CO2 concentration measured at Mauna Loa observatory in parts per million. We ignored latitudinal differences in CO2 concentration as these are small compared to the signal of interest.
Potential evapotranspiration (PET) was calculated using the Priestley–Taylor method as in ref. 27, using inputs of shortwave radiation and the mean of minimum and maximum air temperature from the ANUCLIM archive. The Priestley–Taylor method has been shown to be appropriate for estimating large-scale PET (refs 28,29) and has been adopted in other basin-scale studies14, 30, 31.
Water-balance evapotranspiration. Water-balance evapotranspiration was calculated as the difference of observed annual precipitation and streamflow integrated over the river basin area. The water-balance method remains the most firmly observationally based estimator of ET, but assumes negligible changes in soil water storage at annual to decadal timescales (see Supplementary Section 1 for further discussion). Streamflow time series were acquired from a streamflow collation for unregulated catchments across Australia32. Gaps in the water-balance ET time series (accounting for <5% of monthly records) were filled using simulations from the Australian Water Availability Project33, further detailed in Supplementary Section 1.
Study basins. The 190 study basins were chosen based on the completeness of streamflow records (>95%) and the extent of irrigated and farmed land (<5% of basin area). The basins were classified into wet, sub-humid, semi-arid and arid using the climatological aridity index A (A = PET/P, where PET = annual mean potential ET and P = annual mean precipitation) (see Supplementary Fig. 1 for basin locations and aridity classification). River basins with mean annual aridity index <1 were classified as wet, 1–2 as sub-humid, 2–5 as semi-arid and >5 as arid (adapted from UNEP (1997)34). See Supplementary Section 1 for further details on basin selection and classification criteria.
Breakpoint regression.
Five-year running mean NDVI values were binned according to their corresponding precipitation values. Following ref. 6, the 95th percentile value was determined for each 20-mm-wide precipitation bin separately for each running mean. Breakpoint regression was applied to the 95th percentile values to calculate the first regression slope marking the maximum NDVI attainable for a given precipitation and the breakpoint where the vegetation–precipitation relationship plateaus and vegetation ceases to be water-limited. We then constructed time series of the slopes and breakpoints (Fig. 3) and determined linear trends for both variables. As running means were used to construct the time series, degrees of freedom were adjusted when determining the significance of trends. Farmlands, irrigated areas and wetlands were excluded from this analysis using the Dynamic Land Cover Dataset of Australia35 (see Supplementary Section 1).
CO2 sensitivity coefficients.
Estimation of observed CO2 coefficients. Dimensionless CO2 sensitivity coefficients were calculated from NDVI and ET corrected for precipitation and PET (a function of temperature and shortwave radiation). Precipitation and PET present the main climatic constraints on plant growth36 and are the two first-order controls on ET (ref. 37). The effects of precipitation and PET were removed using linear regression: separately for each basin, annual ET (E) and NDVI were regressed against precipitation and PET and the annual corrected values were calculated as the sum of the regression residual and the 1982–2010 mean of the variable. The corrected annual variables were then log-transformed and regressed against log-transformed annual CO2 concentrations (Ca) to derive the CO2 sensitivity coefficients σET = ∂lnE/∂lnCa and σNDV I = ∂lnNDV I/∂lnCa. The sensitivity coefficients represent the fractional change in the relevant variable per unit fractional change in CO2, so that a change in ET (mm) due to CO2 is well approximated by ΔE/E ≈ σE. ΔCa/Ca for ΔE E and ΔCa Ca (as in this study). ET and NDVI sensitivities to precipitation were calculated from uncorrected data using the same principles (further detailed in Supplementary Section 2).
Prediction of theoretical ET sensitivity to CO2. The theoretical sensitivity of ET (E) to CO2 concentration (Ca) for C3 photosynthesis on a unit leaf area basis can be calculated by writing the CO2 assimilation rate (A) and E in the form of diffusion equations:
and
where gs is the stomatal conductance to CO2, χ is the ratio of internal CO2 concentration (Ci) to Ca, and D is the vapour pressure deficit. χ is a function of D and leaf temperature38, 39 and typically takes values from 0.4–0.5 in arid climates to 0.8–0.9 in wet climates. Substitution of gs from equation (1) into equation (2) yields
Differentiating with respect to Ca, holding D and χ constant, gives:
where σA is the sensitivity of A to Ca:
Equation (4) implies that the sensitivity of E to Ca approaches −1 as the CO2 fertilization effect on A saturates. However, so long as A is increasing with Ca, the sensitivity is smaller in magnitude than −1. The sensitivity of A to Ca can be calculated conservatively by invoking the coordination hypothesis (approximate equality of the carboxylation- and electron transport-limited rates of photosynthesis under field conditions: see, for example, ref. 40). With the further assumption that limitation by the maximum rate of electron transport (Jmax) is not relevant in the field (because Rubisco limitation takes over at the highest light levels), we can express the assimilation rate as
where φ0 is the intrinsic quantum efficiency of C3 photosynthesis, Iabs is the absorbed photosynthetic photon flux density and Γ∗ is the photorespiratory compensation point. Differentiating A with respect to Ca, holding χ constant, gives:
Evaluating equation (7) and then (4) at 25 °C, Ca = 370 ppm for illustration gives σE = −0.61 for χ = 0.8 and −0.38 for χ = 0.5.