Thirty-eight years of CO2 fertilization has outpaced growing aridity to drive greening of Australian woody ecosystems

Climate change is projected to increase the imbalance between the supply (precipitation) and atmospheric demand for water (i.e., increased potential evapotranspiration), stressing plants in water-limited environments. Plants may be able to offset increasing aridity because rising CO2 increases water use efficiency. CO2 fertilization has also been cited as one of the drivers of the widespread “greening” phenomenon. However, attributing the size of this CO2 fertilization effect is complicated, due in part to a lack of long-term vegetation monitoring and interannualto decadalscale climate variability. In this study we asked the question of how much CO2 has contributed towards greening. We focused our analysis on a broad aridity gradient spanning eastern Australia’s woody ecosystems. Next we analyzed 38 years of satellite remote sensing estimates of vegetation greenness (normalized difference vegetation index, NDVI) to examine the role of CO2 in ameliorating climate change impacts. Multiple statistical techniques were applied to separate the CO2-attributable effects on greening from the changes in water supply and atmospheric aridity. Widespread vegetation greening occurred despite a warming climate, increases in vapor pressure deficit, and repeated record-breaking droughts and heat waves. Between 1982– 2019 we found that NDVI increased (median 11.3 %) across 90.5 % of the woody regions. After masking disturbance effects (e.g., fire), we statistically estimated an 11.7 % increase in NDVI attributable to CO2, broadly consistent with a hypothesized theoretical expectation of an 8.6 % increase in water use efficiency due to rising CO2. In contrast to reports of a weakening CO2 fertilization effect, we found no consistent temporal change in the CO2 effect. We conclude rising CO2 has mitigated the effects of increasing aridity, repeated record-breaking droughts, and record-breaking heat waves in eastern Australia. However, we were unable to determine whether trees or grasses were the primary beneficiary of the CO2-induced change in water use efficiency, which has implications for projecting future ecosystem resilience. A more complete understanding of how CO2-induced changes in water use efficiency affect trees and non-tree vegetation is needed.

Abstract. Climate change is projected to increase the imbalance between the supply (precipitation) and atmospheric demand for water (i.e., increased potential evapotranspiration), stressing plants in water-limited environments. Plants may be able to offset increasing aridity because rising CO 2 increases water use efficiency. CO 2 fertilization has also been cited as one of the drivers of the widespread "greening" phenomenon. However, attributing the size of this CO 2 fertilization effect is complicated, due in part to a lack of long-term vegetation monitoring and interannual-to decadalscale climate variability. In this study we asked the question of how much CO 2 has contributed towards greening. We focused our analysis on a broad aridity gradient spanning eastern Australia's woody ecosystems. Next we analyzed 38 years of satellite remote sensing estimates of vegetation greenness (normalized difference vegetation index, NDVI) to examine the role of CO 2 in ameliorating climate change impacts. Multiple statistical techniques were applied to separate the CO 2 -attributable effects on greening from the changes in water supply and atmospheric aridity. Widespread vegetation greening occurred despite a warming climate, increases in vapor pressure deficit, and repeated record-breaking droughts and heat waves. Between 1982-2019 we found that NDVI increased (median 11.3 %) across 90.5 % of the woody regions. After masking disturbance effects (e.g., fire), we statistically estimated an 11.7 % increase in NDVI attributable to CO 2 , broadly consistent with a hypothesized theoretical expectation of an 8.6 % increase in water use efficiency due to rising CO 2 . In contrast to reports of a weakening CO 2 fertilization effect, we found no consistent temporal change in the CO 2 effect. We conclude rising CO 2 has mitigated the effects of increasing aridity, repeated record-breaking droughts, and record-breaking heat waves in eastern Australia. However, we were unable to determine whether trees or grasses were the primary beneficiary of the CO 2 -induced change in water use efficiency, which has implications for projecting future ecosystem resilience. A more complete understanding of how CO 2 -induced changes in water use efficiency affect trees and non-tree vegetation is needed. centrated in the east, where there are large gradients of precipitation (P ) (300-2000+ mm yr −1 ) and potential evapotranspiration (PET) (800-2100 mm yr −1 ). Most eastern Australian woodlands occupy water-limited regions where annual PET far exceeds P (Fig. A1 in the Appendix), and tree species have evolved to cope with water-limited conditions (Peters et al., 2021) and high interannual rainfall variability. However, the climate is warming: 8 of Australia's 10 warmest years on record have occurred since 2005 (CSIRO and Bureau of Meteorology, 2020), and Australia's climate has warmed by ∼ 1.5 • C since records began in 1910. The warming has likely increased atmospheric demand for water (e.g., PET or vapor pressure deficit, VPD). In most woody ecosystems, the ratio of water supply (i.e., P ) to water demand (i.e., PET) has declined in recent decades (Figs. 1,2a). Eastern Australia has also been impacted by several multiyear droughts, episodic deluges of rainfall (King et al., 2020), and an increasing frequency of severe heat waves (Perkins et al., 2012) in the last few decades. Precipitation changes have been spatially variable over eastern Australia, where northern Queensland grew wetter, and southeastern Australia grew drier (Fig. 2a). In the last 2 decades, southeastern Australia experienced the two worst droughts in the observational record (2001-2009; van Dijk et al., 2013;and2017-2019;Bureau of Meteorology, 2019). Yet between these two droughts, eastern Australia experienced record-breaking rainfall in 2011 associated with a strong La Niña event. This caused marked vegetation "greening" (e.g., increased foliar cover), even in the arid interior (Bastos et al., 2013;Poulter et al., 2014;Ahlstrom et al., 2015). However, this greening contributed to record-breaking fires in the following year (Harris et al., 2018).
Theory suggests that plant physiological responses to atmospheric carbon dioxide (CO 2 ) may mitigate some of the negative effects of an aridifying climate. However, the magnitude of plant responses to increased atmospheric CO 2 has been challenging to establish in field experiments (Jiang et al., 2020b) and from observations (Frankenberg et al., 2021;Zhu et al., 2021;Walker et al., 2020) or to separate from other drivers (e.g., climate variability, disturbances, and changes in land management; Zhu et al., 2016). Studies have used data from the Advanced Very High Resolution Radiometer (AVHRR) satellites to show positive trends in the normalized difference vegetation index (NDVI) over Australia (Donohue et al., 2009). The greening trend is caused by increased leaf area, which has likely resulted from increased atmospheric CO 2 concentrations (Donohue et al., 2013;Ukkola et al., 2016). The evidence for increases in leaf area from rising CO 2 has also been supported by observations of reduced runoff in Australia's drainage basins (Trancoso et al., 2017;Ukkola et al., 2016).
Yet disentangling the CO 2 fertilization effect from other drivers of climate variability and global change has been particularly challenging for satellite-based analyses. It is challenging to attribute causes of greening because co-occurring changes in climate, land use, and disturbance are confounded by the effect of CO 2 fertilization. Furthermore, the time series of even the longest systematically collected optical vegetation index records from a single sensor is 20 years (e.g., MODIS Terra). Analysis of trends extending beyond 20 years requires merging satellite records across sensors and platforms. But this requires care to address changes in radiometric and spatial resolution of the sensor as well as drift in the solar zenith angle (Ji and Brown, 2017;Frankenberg et al., 2021) and the time of retrieval. Thus different analytical methodologies have produced disagreements over where greening has occurred (Cortés et al., 2021). One oftenused method to provide additional constraint on greening trends has been to compare remote-sensing-derived trends with modeled changes in leaf area index (LAI) from ensembles of dynamic global models (Zhu et al., 2016;Wang et al., 2020). However these model attribution approaches rely on a set of key assumptions. None of the models can accurately predict LAI changes in response to rising CO 2 (De Kauwe et al., 2014;Medlyn et al., 2016). Vegetation models have been shown to diverge in their simulation of LAI over Australia (Medlyn et al., 2016;Teckentrup et al., 2021;Zhu et al., 2016) and have bioclimatic rules for determining phenology which may not be appropriate for the highly variable Australian climate and the evergreen Eucalyptus forests (Teckentrup et al., 2021). These model simulations are typically compared with modeled LAI products derived from the red and near-infrared wavelengths of multispectral satellite sensors, of which each product carries specific algorithmic assumptions about canopy-light interception which are conditional upon estimated land cover types. In comparison, NDVI carries no ecosystem-specific assumption and is an effective proxy for leaf area in ecosystems with low to moderate canopy cover (Carlson and Ripley, 1997), a characteristic of eastern Australian woody ecosystems (Specht, 1972;Yang et al., 2018).
Here we ask how much greening trends can be explained by rising CO 2 . Using eastern Australia as a model system, we used a multi-satellite-derived NDVI record encompassing 38 years with multiple statistical techniques to isolate the influence of CO 2 from simultaneous effects of meteorological change and disturbance. Next we contrasted CO 2 effects with theoretical predictions based on water use efficiency (WUE) theory for plants and the observed rise in CO 2 . Finally, we examined whether recent NDVI greening trends have co-occurred with changes in tree or grass cover over the last 2 decades.

Climate and remote sensing datasets
We used the atmospheric CO 2 record from the deseasonalized Mauna Loa record (https://www.esrl.noaa.gov/gmd/ ccgg/trends/data.html, last access: 5 April 2020) and extracted climate data (Table A1) from the Australian Bureau of Meteorology's Australian Water Availability Project (AWAP; Jones et al., 2009). AWAP is a gridded climate product interpolated to 0.05 • from a large network of meteorological stations distributed across Australia. Vapor pressure deficit was calculated using daily estimates of maximum temperature and vapor pressure at 15:00 LT (UTC+10). PET was calculated from shortwave radiation and mean air temperature using the Priestley-Taylor method (Davis et al., 2017). The Priestley-Taylor method has been shown to be appropriate for estimating large-scale PET (Raupach, 2000) and is more suited for use in long-term analysis where CO 2 has increased than other common formulations such as the Penman-Monteith equation (Greve et al., 2019;Milly and Dunne, 2016), which explicitly imposes a fixed stomatal resistance that is incompatible with plant physiology theory (Medlyn et al., 2001). AWAP measurements of shortwave radiation only extend back to 1990, so we extended the PET record to 1982 by calibrating the ERA5-Land PET record  to the AWAP PET record (1990-2019) by linear regression for each grid cell and then gap-filled the years 1982-1989 with the calibrated ERA5 PET. PET from the Climate Research Unit record (Harris et al., 2014) was highly correlated with both the recalibrated ERA5 PET (r = 0. 91,[1982][1983][1984][1985][1986][1987][1988][1989] and the original AWAP PET (r = 0.97, 1990-2019). Next, we calculated a 30-year climatology of the meteorological variables using the period of 1982-2011 to be close to current standards (World Meteorological Organization, 2017). We used this climatology to define the mean annual P : PET (the moisture index, MI MA ) and as the reference to calculate a 12-month running anomaly of annual P : PET (MI anom ). Zonal statistics for each meteorological variable were calculated using simplified Köppen climate zones, derived from the Australian Bureau of Meteorology (Fig. 2b, Table A1).
We used surface reflectance from two satellite products to generate the NDVI record: National Oceanic and Atmospheric Administration's Climate Data Record v5 Advanced Very High Radiometric Resolution (AVHRR) Surface Reflectance (NOAA-CDR) and the National Aeronautics and Space Administration's MCD43A4 Nadir Bidirectional Reflectance Distribution Function Adjusted Reflectance (MODIS-MCD43) (Table A1; Schaaf and Wang, 2015). NDVI data were extracted from 1982-2019 at 0.05 • resolution from the NOAA-CDR AVHRR version 5 product (Vermote and NOAA CDR Program, 2018). The surface reflectance record of AVHRR extends through 2019, but the quality of the record starts to degrade in 2017 because of an increase in the solar zenith angle (Ji and Brown, 2017), causing a sensor-produced decline in NDVI during 2017-2019. For this reason we only use AVHRR surface reflectance data between 1982-2016. We composited monthly mean AVHRR NDVI (NDVI AVHRR ) estimates using only daily pixel retrievals with no detected cloud cover (quality assurance band, bit 1). Monthly NDVI AVHRR estimates aggregated from fewer than three daily retrievals were removed. They were also removed when the coefficient of variation in daily retrievals for a given month was greater than 25 %. We also removed NDVI AVHRR monthly estimates where NDVI AVHRR , solar zenith angle, or time of acquisition deviated beyond 3.5 standard deviations from the monthly mean, calculated from a climatology spanning 1982-2016.
We used the MODIS-MCD43 surface reflectance at 500 m resolution to derive NDVI for 2001-2019 (NDVI MODIS ). Monthly mean estimates of the surface reflectance were produced by compositing pixels flagged as "ideal quality" (quality assurance, bits 0-1). We also masked disturbances to have greater confidence in our attribution of the targeted drivers of NDVI MODIS change (climate and CO 2 ). The Global Forest Change product v1.7 (Hansen et al., 2013) was used to mask pixels from 2001 onwards that had experienced forest loss due to deforestation or severe stand clearing disturbance. We masked pixel locations that experienced bushfires from the year 2001 onwards. Specifically, these pixels were masked for the year of burning and the following 3 years using the 500 m resolution MODIS-MCD64 monthly burned-area product (Giglio et al., 2018). We terminated the NDVI MODIS time series in August of 2019, prior to the widespread bushfires of late 2019/early 2020. Both NDVI AVHRR and NDVI MODIS datasets were processed using Google Earth Engine (Gorelick et al., 2017) and exported at 5 km spatial resolution, which best approximated the native resolution of the NOAA-CDR AVHRR and AWAP products. Further post-processing used the "stars" (Pebesma, 2020) and "data.table" (Dowle and Srinivasan, 2019) R packages (see "Code and data availability" section).
We merged the processed 1982-2016 NDVI AVHRR with the 2001-2019 NDVI MODIS by recalibrating the NDVI AVHRR with a generalized additive model (GAM). Specifically, we used 1 million observations from the overlapping 2001-2016 494 S. W. Rifai et al.: CO 2 -driven greening portion of both records to fit a GAM using the "mgcv" R package (Wood, 2017) to model NDVI MODIS from AVHRRderived covariates as where "s" represents a penalized smoothing function using a thin plate regression spline; SZA is the solar zenith angle; NDVI AVHRR is the uncalibrated NDVI from AVHRR; TOD is time of day of retrieval; and x and y represent longitude and latitude, respectively. The fit GAM was then used to generate the recalibrated AVHRR NDVI. The merged NDVI dataset was created by joining the 1982-2000 recalibrated AVHRR NDVI with the 2001-2019 NDVI MODIS . We further reduced monthly temporal variability in NDVI by calculating a 3-month rolling mean of NDVI, which we used for subsequent statistical model fitting.

Estimating NDVI and climate trends
We estimated the relative increase in NDVI between 1982-2019 with respect to time (Eq. 2) for each grid cell with an iteratively weighted least squares robust linear model via the "rlm" function in R's MASS package (Venables and Ripley, 2002) as follows.
Here β 0 represents the estimated NDVI in 1982, the year term starts at 1982, and the sensor term is a binary covariate that accounts for residual offset differences between the recalibrated AVHRR NDVI and the NDVI MODIS . The relative temporal trends for climate variables and the MODIS vegetation continuous fractions were fit for each grid cell location using the Theil-Sen estimator, a form of robust pairwise regression, with the "zyp" R package (Bronaugh and Werner, 2019). The temporal covariate was recentered to start with the first hydrological year (where the year starts 1 month earlier in December) of the data so that the intercept term represents the mean at the start of the time series. The relative rate of change for each variable was reconstructed by calculating where β 0 and β 1 are the intercept and trend derived from Theil-Sen regression.

Estimating contribution of CO 2 and climate toward NDVI trends
We fit six forms of statistical models to the merged NDVI observations in order to quantify the impact of changes in CO 2 and meteorological variables on NDVI. Figure 1 shows that the relationship between NDVI and multi-annual mean of P : PET was nonlinear and followed a monotonic saturating sigmoidal relationship. GAMs can characterize a nonlinear response without specifying a functional form, yet the underlying spline parameters are not easily interpreted as the parameters of a fixed nonlinear function. Therefore we fit models with a set of fixed nonlinear functional forms using nonlinear least squares (nls.multstart package in R v4.01; Padfield and Matheson, 2020) and compared models fit using all pixel locations. The nonlinear functional forms included the Weibull function (Eq. 4, Fig. 5), the logistic function (Eq. 5, Fig. A5), and the Richards growth function (Eq. 6 Fig. A6). We focused on the Weibull models because they showed equivalent goodness of fit with fewer parameters than the Richards function models. Next we added a linear modifier to the Weibull function using the covariates of CO 2 (µmol µmol −1 ; C a ) and the ratio of the anomaly of P : PET (MI anom ) to the mean annual P : PET (MI MA ) as follows: Here the sensor term is a binary covariate indicating the AVHRR or MODIS sensor. Model-fitted parameters V a and V d correspond to the asymptote and the asymptote's difference from the minimum NDVI, while c ln is the logarithm of the rate constant, and q is the power to which MI MA is raised. The model was fit by individual season with 1 million observations per model fit. Corresponding goodness-of-fit metrics (R 2 and root mean square error) were calculated by season ( Fig. 5) with 1 million randomly sampled observations. We tested alternative nonlinear functional forms to characterize the effect of CO 2 upon NDVI. A logistic model was fit across space for each hydrological year as where NDVI is the hydrological year mean value of NDVI for a grid cell location, m is the midpoint, s is a scale parameter, and V A is the asymptote (plotted in Fig. A5). We also used a modified Richards growth function to characterize the CO 2 effect upon seasonal NDVI (Fig. A6) as Here the β terms act to linearly modify the core nonlinear parameters (V A , m, s, q) with the effects of CO 2 and MI f.anom . Each seasonal model component was fit across space with 1 million random samples from the total merged NDVI record (approximately 14.3 million observations).
To ensure consistent interpretation of the nonlinear response across P : PET, we also fit linear models explaining NDVI with CO 2 and MI anom by season in MI MA bin widths of 0.2 (Eq. 7, Fig. A4). Separate linear models were fit for increments of 0.15 of MI MA for each season using the merged 1982-2019 NDVI record. NDVI was modeled as where Veg. Class is the NVIS 5.1 vegetation class, and sensor is a binary variable used to account for residual differences between the recalibrated AVHRR NDVI and NDVI MODIS records. To aid the comparison of model effects, we centered and standardized the continuous model covariates before regression. The standardized CO 2 and P : PET anom effects (β) are presented in Fig. A4. Next, to produce spatially varying estimates of the CO 2 effect, we fit robust multiple linear regression models to the time series of NDVI for each of the 39 463 pixel locations. The CO 2 effect for each grid cell location was simultaneously estimated with the linear effects of the anomalies (anom) of P , PET, VPD, and MI as fractions of their mean annual (MA) values as follows: Finally we produced a GAM to estimate the CO 2 effect in addition to the nonlinear interactions between CO 2 , mean annual climate, and climate anomalies. Here instead of predetermined nonlinear functional forms, the GAM used penalized smoothing functions (s) to estimate the potentially nonlinear effect of the covariates. The GAM was fit across all pixel locations and estimated the CO 2 effect and the effects of the anomalies and mean annual values of VPD, P , and PET as well as sensor epoch as follows:

A simplified theoretical water use efficiency model
We compared the statistically attributed CO 2 amplification of NDVI with the expectation from a simple theoretical model of WUE. Following Donohue et al. (2013), WUE (W ) is defined as where A is leaf-level carbon assimilation (µmol m 2 s −1 ), E is leaf-level transpiration (µmol m 2 s −1 ), C a is atmospheric CO 2 (µmol µmol −1 ), C i is intercellular CO 2 (µmol µmol −1 ), χ is C i C a , and D is atmospheric vapor pressure deficit (mol mol −1 ). The relative rate of change in W with respect to a change in C a can be calculated as . (11) If temperature increases without a corresponding increase in humidity, D increases, which also causes transpiration to rise and thus reduces W . However, W is predicted to increase with CO 2 , which may offset increases in D. Experiments suggest that χ does not change with C a but is sensitive to D (Wong et al., 1985;Drake et al., 1997) and can be estimated as being proportional to the square root of D (Medlyn et al., 2011). By substituting 1 − χ ≈ √ (D) into Eq. (11) we can estimate the theoretical combined effect of C a and D upon W leaf as Transpiration per unit ground area is strongly controlled by water supply in warm, water-limited environments with relatively low leaf area such as eastern Australia (Specht, 1972); therefore we approximate canopy transpiration (E canopy ) as where L is leaf area. The change in E canopy can then be defined as If we assume there is no overall change in precipitation then we can assume change in E canopy is tightly coupled to the water supply; therefore we have We note that the assumption of constant precipitation is obviously unrealistic, most notably because two major droughts occurred during the observation period. However, over much longer periods (i.e., 1910-2018) precipitation in eastern Australia has been relatively stationary . NDVI is linearly related to foliar cover (F ) until LAI ≈ 3 (m 2 m −2 ) (Carlson and Ripley, 1997), which is predominantly the case when P : PET < 1. Most woody ecosystems in eastern Australia are strongly water-limited with LAI ≤ 1 (m 2 m −2 ), where NDVI is approximately proportional to the fraction of foliar cover: Then substituting Eq. (15) into Eq. (12) gives If we assume that the benefit towards W leaf from rising C a is split evenly between the relative changes in A leaf and F , we can predict the change towards NDVI to be We compared the WUE theoretical model with the robust linear models fit for each pixel location (Eq. 8) and the GAM (Eq. 9) fit across the study region. The WUE theoretical model assumes no change in P but does account for changes in VPD. Therefore in using the statistical models to compare with the WUE predictions, we generated counterfactual predictions from the statistical models with no precipitation anomaly but with the observed increases in CO 2 and VPD. One weakness with the application of this WUE theoretical model is the uncertainty regarding the assumed allocation of the W leaf benefit towards either A leaf or F (e.g., LAI; see above). Donohue et al. (2017) proposed a similar model to Eq. (18), the partitioning of equilibrium transpiration and assimilation (PETA) hypothesis, where the relative allocation to leaf area is predicted to decline with increasing resource availability (which could be inferred from growing season LAI). We calculated the expectation from the PETA hypothesis as another point of comparison with the CO 2 -attributable effect on NDVI.

Empirical attribution of the CO 2 effect
We found consistently positive NDVI responses to CO 2 across the moisture gradient of P : PET for all seasons, with the greatest increases located in regions of higher P : PET (> 0.5) (Fig. 5). The nonlinear Weibull models showed a larger CO 2 -attributable effect on NDVI in regions of higher P : PET (Fig. 5), but the effect size of the NDVI response to CO 2 was largely consistent across model forms (Figs. A4-A7). The CO 2 -attributable increase in NDVI between 1982-2019 ranged from approximately 5 % in the Arid interior regions to > 20 % in the wettest Tropical and Temperate regions (Fig. 5a). This was consistent with linear model forms when fit for individual grid cell locations (Eq. 8, Fig. 6) as well as by comparing the CO 2 effect size across 16 linear models fit for grid cell locations grouped into bins spanning 0.1 increments of P : PET (Eq. 7, Fig. A4). The GAM fit across grid cell locations also indicated a larger CO 2 effect in regions with higher P : PET (Figs. 6, 7b). Quantile regression with generalized additive models showed a pronounced response to CO 2 across the distribution of pixels with both low and high NDVI (10-97.5 percentiles) across the full aridity gradient of P : PET (Fig. A7).
A recent study found that the global CO 2 fertilization effect was halved between the 1980s and 2000s (Wang et al., 2020). In contrast to the estimates over eastern Australia from Wang et al. (2020), we found no consistent evidence of a decline in the effect of CO 2 on NDVI through time. Neither the GAM estimates nor the robust linear model estimates of the CO 2 effect showed any consistent evidence of a weakening CO 2 effect between 1982-2000 and 2001-2019 (Fig. 6). The central 25th-75th percentiles of the distribution of robust linear model effect sizes overlapped in all regions between epochs. The central 25th-75th percentiles of the GAM-estimated distributions also overlapped, with the exception of the Grassland and Arid regions, where the CO 2 effect was larger during 2001-2019. Consistent with the finding of a greater CO 2 effect in wetter regions (Fig. 5), the robust linear models and the GAM estimated the CO 2 effect to be greatest in the Equatorial and Tropical regions and lowest in the Arid region (Fig. 6).

CO 2 -driven greening and expectations from water use efficiency
The range of statistically estimated CO 2 -attributable greening responses were compared with the expectation from the theoretical CO 2 water use efficiency model. The Donohue et al. (2013) CO 2 ×WUE model (see "Methods") accounts for changes in VPD but assumes no change in water supply. Assuming an equal split in the benefits of WUE between greater carbon assimilation and increased foliage cover (see below for alternative assumptions), the model predicted an 8.7 % (10th-90th percentile range [+6.8 %, +10.2 %]) increase in NDVI (proxy for foliage cover; see "Methods"). This compared to an estimated 11.7 % ([+4.6 %, +14.6 %]) relative increase from the GAM when accounting for a simultaneous increase in VPD (which the WUE accounts for) and factoring out the effects of changing precipitation and PET (which the WUE model does not account for). We needed to assume differing levels of allocation to foliar gain (i.e., not 50 %) for the theoretical WUE model to match the statistically estimated CO 2 effect (Fig. 7e). Regions of higher P : PET (Equatorial, Tropical, Temperate, and Temperate Tasmanian) required greater allocation fractions than 50 %, whereas the allocation fraction would be between 25 %-50 % for regions with lower P : PET (Arid, Grassland, and Subtropical). In comparison, the PETA hypothesis (Donohue et al., 2017) predicted the greatest CO 2 effect on leaf area to be in regions with the lowest LAI, but this was not supported by the statistically estimated CO 2 effect on NDVI (Fig. A9). Despite having the lowest LAI, the Arid region received the smallest CO 2 effect, but it is worth noting that the Arid region also experienced the greatest increase in VPD and reduction in P over the 38-year period (Figs. 7a, A8). In contrast, the largest estimated CO 2 -attributable effects on greening were found to be in the Equatorial and Tropical areas.

Co-occurring shifts in aridity, NDVI, and vegetation cover
The shifts in P : PET and NDVI were accompanied by vegetation cover changes in some regions. Most notably, the Arid and Temperate regions experienced the strongest zonally averaged declines in P : PET (Fig. 8a) and increases in VPD (Fig. A8). Seasonal greening trends were relatively similar apart from the aforementioned exceptions in the Grassland. The MODIS vegetation continuous fraction data from 2001-2018 indicated that most regions experienced modest changes in tree vegetation cover. The largest decline in tree cover occurred in the Temperate regions (Figs. 8c, A10). Most regions experienced declines in non-vegetated (bare) cover, increases in non-tree vegetation, and modest change in tree cover (Fig. 8c); however the proportional increase in non-tree vegetation typically exceeded tree cover increases (Fig. A10).

Australian woody vegetation as model systems to quantify CO 2 fertilization
Australia is ideal to explore the CO 2 contribution towards vegetation greening because there are fewer confounding effects to drive greening. Forest trees are evergreen, and the growing season is rarely limited by temperature or radiation. The study region spans a large moisture gradient (Fig. A1a), but unlike much of the global tropics, the large majority of the study region is not so cloudy as to prevent multiple high-quality multispectral satellite retrievals per month. Australia has also not been subjected to other prominent drivers of greening such as nitrogen deposition (Ackerman et al., 2019). Nevertheless, Australia has experienced notable land use change during the study period such as high rates of deforestation in Queensland and northern New South Wales (Evans, 2016). However, we excluded affected pixel locations from the analysis.
Prior global analyses on warm arid environments have quantified the CO 2 effect on greening (Donohue et al., 2013). More expansive global-scale analyses of all terrestrial vegetation using dynamic global vegetation model attribution have found evidence to connect Australia's greening to changes in atmospheric CO 2 concentration (Winkler et al., 2021), while an earlier study did not (Zhu et al., 2016). Australian studies documented the greening trend up to 2010 using the long-term AVHRR record (Donohue et al., 2009;Ukkola et al., 2016) and have been able to partially attribute CO 2 as a driver of greening in sub-humid and semi-arid regions. Here we advanced upon prior research to separate the effects of disturbance and changes in aridity and moisture in order to quantify the CO 2 fertilization effect across the full spectrum of moisture availability experienced by Australian woody ecosystems, notably for 38 years.

Regional differences in greening and browning through time
Despite the region's high decadal-scale variability in NDVI (Fig. 4), the nearly 4-decade-long record allowed us to separate the CO 2 effect on NDVI from the anomalies caused by drought (e.g., 2003-2009) or high rainfall (e.g., La Niña 2010-2011). Although the long-term greening trends we document in the 9 years following earlier studies are generally consistent (Donohue et al., 2009;Ukkola et al., 2016), our results diverge and lead to key differences in interpretation of why NDVI has continued to increase.  (Fig. A8). These increases were largest in the most Arid and Temperate regions (12.7 % and 11 % since 1982, respectively; Fig. 6), and when coupled with seasonal reductions in precipitation (Fig. A2) these would have partially offset benefits from increased intrinsic water use efficiency (Eq. 17). In contrast, more than half of the Equatorial, Tropical, Subtropical, and Grassland regions in northern Queensland experienced increases in precipitation since 1982 ( Fig. A2; Ukkola et al., 2019), thus allowing these locations to exceed the predicted NDVI increases from the WUE model (Fig. 7d).
It should be noted that not all regions experienced consistent greening trends throughout the observation period. For example, "greening" shifted to "browning" during the austral summer (December-February) in the Arid and Grassland regions of Queensland between the 1982-2000 and 2001-2019 records (Fig. 3b). It is unclear why browning occurred during austral summertime in the Grassland region (Figs. 3b, 4, 8b). The declines in NDVI during 2001-2019 may have been meteorologically driven and related to shifts in the distribution of wet-and dry-season precipitation (Fig. A2). Alternatively, the shift could be due to changes in fire and cattle management that have been particularly prevalent across regions of Queensland in recent decades (Seabrook et al., 2006). Further, greening may have been suppressed in parts of Queensland because cattle ranching activity has intensified and has

Attributing a CO 2 fertilization contribution towards greening
Plants increase their rates of photosynthesis in response to rising atmospheric CO 2 whilst also reducing stomatal conductance, which reduces evaporative losses, and combined, the two responses lead to greater WUE (Ainsworth and Rogers, 2007;Morison, 1985). In water-limited ecosystems, it has been hypothesized that this physiological response by plants to CO 2 should result in increased leaf biomass (Donohue et al., 2013;Ukkola et al., 2016). While all of our statistical approaches indicated a year-round positive CO 2 effect, in contrast to theory the effect was consistently greater in regions with higher P : PET (Figs. 5-7, A4-A7). Our analysis also diverges with a global coarse-scale analysis that found weakening of the CO 2 fertilization effect across both southeastern and northern Australia (Wang et al., 2020). We found no meaningful difference in the CO 2 -attributable effect towards greening between the AVHRR 1982-2000 and MODIS 2001-2019 epochs to support the finding of a temporally weakening CO 2 effect (Fig. 6). The WUE model predicted similar relative rates of NDVI increase between the two epochs (4.5 % and 4.7 % for AVHRR and MODIS), yet for different reasons. CO 2 increased by 30 ppm between 1982-2000, while VPD changed minimally, and precipitation increased in Queensland and the Arid region (Fig. A8), all of which are favorable to increasing NDVI. In contrast, the larger CO 2 increase between 2001-2019 (40 ppm) was offset by a spatially ubiquitous increase in VPD (3.7 %, CI = [−0.4 %, 8.6 %]; Fig. A8) and the occurrence of two multi-year droughts (Millennium Drought, 2003, and The Big Dry, 2017-2020. Despite the widespread evidence Figure 6. Boxplot representation of the CO 2 -attributable effect upon changing NDVI. Here the 25th, 50th, and 75th percentiles of the CO 2attributable effect on NDVI are shown. The robust linear models (RLM; see "Methods", Eq. 13) were fit for each individual grid cell location, whereas the generalized additive model (GAM; see "Methods", Eq. 14) was fit using all grid cell locations. The distribution of RLMs yielded a median R 2 of 0.58 and RMSE of 0.025 over the merged period. The GAM had an overall R 2 of 0.91 and RMSE of 0.049.
of the CO 2 effect, the WUE model notably underpredicted greening in some regions (Fig. 7).

Deviations from WUE
We found that the CO 2 effect on foliar area was the smallest in the driest climate regions (Figs. 5-6, A4-A7). This was at odds with the WUE prediction (at 50 % allocation; Fig. 7e) and contrary to the expectation that the greatest WUE-derived benefit from CO 2 would be in drier climates (Donohue et al., 2017;McMurtrie et al., 2008). These deviations may have resulted from ecosystems processes beyond the scope of a simple model, such as more severe nutrient limitations, the phenology of the vegetation composition, belowground root processes, and disturbances not captured by the satellite products (e.g., small fires, grazing). Browning in the Grassland region (Figs. 3b, A3) may have been caused by distinct dry-season phenological differences between overstory woody vegetation and understory C 4 grasses (Moore et al., 2016), which are dominant there (Murphy and Bowman, 2007). C 4 grasses may have been favored over C 3 because precipitation increased during the austral summer over the course of the study (Fig. A2), and this higher concentration of rainfall during the warmest months is thought to favor C 4 grasses (Hattersley, 1983;Knapp et al., 2020;Murphy and Bowman, 2007). Finally, the linear dependency of leaf area upon VPD in the theoretical model may be ill-suited for extreme anomalously arid conditions because NDVI observations suggest a strongly nonlinear relationship with large VPD anomalies (Fig. A11). We explored how much of the higher WUE benefit would have to be allocated to match the GAM-estimated CO 2attributed changes in NDVI. Allocation rates far greater than 50 % would be required to match the WUE prediction in the Tropical and Equatorial zones (Fig. 7e), whereas allocation would need to be less than 50 % to match the GAM estimate in the Arid region. The Arid region experienced the greatest relative increase in VPD (Figs. 7a, A8), yet the theoretical model still predicted a small but positive increase in NDVI (Fig. 7c), which the GAM estimate suggests would be closer to a 10 % rather than 50 % foliar allocation level from the WUE benefit (Fig. 7d). Nevertheless, the smaller effect over the Arid region was consistent with earlier observational findings across Australia (Ukkola et al., 2016) and Figure 7. Long-term changes in vapor pressure deficit and comparison between the NDVI due to CO 2 and the expected CO 2 fertilization effect on foliar area due to gains in water use efficiency (WUE). (a) The relative increase in the annual mean of vapor pressure deficit (VPD) between 1982-2019. (b) The predicted relative increase in NDVI due to CO 2 (see "Methods", Eq. 14) with a concurrent increase in VPD. (c) The relative expected increase in NDVI following the theoretical WUE prediction where 50 % of the gain is allocated to foliar area (see "Methods", Eq. 11). (d) The difference between the predictions of relative NDVI from the WUE model with 50 % allocation and the GAM-estimated relative increase caused by CO 2 . Blue regions indicate where the GAM prediction was greater than the WUE prediction.
(e) Boxplot of the expected relative NDVI % increase from WUE gains due to CO 2 fertilization depending upon differing levels of foliar allocation from the CO 2 benefit toward WUE and the GAM-estimated NDVI% increase due to CO 2 (see "Methods", Eq. 11). The boxplot distribution of the Equatorial region appears collapsed because the Equatorial portion of the region is relatively small (a).  .

Relation to ecosystem CO 2 fertilization experiments
Notably, a 4-year-long ecosystem-scale CO 2 manipulation experiment carried out in a mature Eucalyptus woodland in Sydney (EucFACE) did not observe an increase in leaf area under elevated CO 2 (Jiang et al., 2020b). The experimental site is located upon phosphorus poor soils, typical of Australia. The lack of a leaf area growth response observed at EucFACE is not necessarily inconsistent with the greening effect observed in this study. Our observational window of 38 years is much longer than the elevated CO 2 exposure time in the experiment (4 years in Jiang et al., 2020b), covers different CO 2 increments (historical vs. future), and could imply that woodlands are eventually able to liberate belowground phosphorus to support greater biomass growth on longer timescales. Increased autotrophic soil respiration and belowground productivity were observed at EucFACE under elevated CO 2 exposure Jiang et al., 2020b), as was a brief period of enhanced nitrogen and phosphorus mineralization (Hasegawa et al., 2016). Over time, this increased investment of carbon belowground could potentially liberate sufficient phosphorus to support an expansion of leaf area. The reduced allocation to foliar area in the Arid region (Fig. 7e) may reflect that extra carbon derived from CO 2 is allocated belowground to increase water uptake or mitigate other resource limitations such as soil phosphorus (Jiang et al., 2020a).

Vegetation composition shifts
Most grid cell locations in the Temperate zone experienced simultaneous apparent declines in tree cover and increases in non-tree vegetation cover (e.g., grasses and shrubs) (Figs. 8c, A10). This is surprising because we focused this regression analysis on 2001-2018 in order to exclude the reduced tree cover due to the catastrophic megafires of 2019/20 . Some may question the ability of the MODIS Vegetation Continuous Fraction product (DiMiceli et al., 2017) to accurately distinguish Australian tree cover from non-tree vegetation (Sexton et al., 2013;Adzhar et al., 2021); however this pattern is consistent with a recent lidar-derived tree cover time series of Australia (Liao et al., 2020). Future work may seek to uncover trends in the green vegetation fraction by analyzing higher-resolution data derived from the Landsat constellation, such as Geoscience Australia's vegetation fractional cover product (Gill et al., 2017). The decline in tree cover suggests that the drought starting in 2017 was already killing trees prior to the 2019/20 megafires. Field observations of tree decline remain relatively rare, but a citizen science initiative has documented more than 300 locations of non-fire-related mass tree mortality between 2018-2020 (Atlas of Living Australia, 2020). Similarly, a study using an experimental constrained plant hydraulics model to predict the regions at risk of drought-induced tree mortality found greater risk in the same arid regions of southeastern Australian forests and woodlands . These predicted regions of mortality coincide with where we document the greening trends that fell short of the theoretical WUE expectation. These rapid shifts in vegetation underlie the need for greater continuous field vegetation monitoring to capture change imposed by climate extremes.

Conclusions
We separated the effects of disturbance and meteorological anomalies with statistical models to show that increasing CO 2 produced nearly 4 decades of widespread vegetation greening across eastern Australia. The large agreement between a theoretical model and the statistically estimated CO 2 effect indicated that greening resulted from an increase in water use efficiency. Vegetation greening occurred despite a highly variable and increasingly arid climate and on soils particularly poor in phosphorus, which have likely acted as a constraint on growth. While rising atmospheric CO 2 ameliorated what would have been a browning woody ecosystem response to declining P : PET, the CO 2 effect was insufficient to promote greening when both P and P : PET experienced long-term decline, as observed in the more arid regions in our study. Further, it is unknown whether further increases in atmospheric CO 2 will continue to enable vegetation to mitigate increases in aridity and VPD under future warming. It is also unclear if trees or grasses are the primary contributors to the recent greening trend. Future localized work is urgently needed to better understand recent changes in tree and grass competition under an increasingly arid climate, which will be essential to help forecast ecosystem resilience. Finally, our results have important implications for understanding Australia's terrestrial water availability. Greening trends signal changes in evapotranspiration and runoff and therefore need to be considered in planning for future land and water resource management on the world's driest inhabited continent.         Figure A9. (a) The fractional increase in NDVI from CO 2 is plotted following predictions from the partitioning of equilibrium transpiration and assimilation (PETA) hypothesis (green), the Generalized Additive Model (GAM; see "Methods", Eq. 6) estimates, and the pixel-level Robust Linear Model (RLM; see "Methods", Eq. 5) estimates. Note that we substitute NDVI for the leaf area formulation of the PETA hypothesis (Donohue et al., 2017) because we assume a linear relationship between relative NDVI increase and relative leaf area increase (Donohue et al., 2013;Carlson and Ripley, 1997). (b) The relative change in NDVI between 1982-2019 across the range of mean annual leaf area index calculated with MCD15A3H from 2002-2019. Figure A10. The relative shifts in land fraction by tree cover, non-tree vegetation cover (predominantly grasses), and bare ground. Annual rate (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) of change for tree cover (a), non-tree vegetation (b), and discretized change classes (c). Min change indicates that the annual rate of change in percentage cover was between −0.1 % yr −1 and 0.1 % yr −1 . Figure A11. The nonlinear dependency of NDVI upon the annual anomaly of vapor pressure deficit. Two-dimensional density plot of NDVI and the running 12-month anomaly of vapor pressure deficit. The red line shows a GAM fit. Code and data availability. All data used are publicly available from sources listed in Table A1. A Git repository with all associated data processing, analysis, and code to reproduce figures is located at https://github.com/sw-rifai/ eastern-Australia-CO2-NDVI-change (last access: 18 November 2021) and https://doi.org/10.5281/zenodo.5711964 (Rifai, 2021a). Processed data used in model fitting and figures can be accessed via the Zenodo data repository: https://doi.org/10.5281/zenodo.4340064 (Rifai, 2021b).
Author contributions. SWR, MGDK, PM, LAC, BEM, and AJP designed the study. SWR analyzed the data and produced the figures. SWR and AMU processed climate data. SWR wrote the first draft with MGDK, and all authors have contributed to writing and revising the manuscript.
Competing interests. At least one of the (co-)authors is a member of the editorial board of Biogeosciences. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.