Thirty-eight years of CO2 fertilization have 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 interannual to decadal-scale climate variability. In this study we 5 asked the question, how much has CO2 contributed towards greening? We focused our analysis on a broad aridity gradient spanning eastern Australia’s woody ecosystems. Next we analysed 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 10 pressure deficit, and repeated record-breaking droughts and heatwaves. 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 15 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. 1 https://doi.org/10.5194/bg-2021-218 Preprint. Discussion started: 23 August 2021 c © Author(s) 2021. CC BY 4.0 License.

1 Introduction 20 Australia is the world's driest inhabited continent. Predicting how climate change will affect ecosystem resilience and alter Australia's terrestrial hydrological cycle is of paramount importance. Australia's woody ecosystems are mostly concentrated 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), and tree species have evolved to cope with water-limited conditions (Peters et al.) and high interannual rainfall variability. 25 However, the climate is warming: eight of Australia's ten 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 multi-year droughts, episodic deluges of rainfall (King et al., 2020), and an increasing frequency of 30 severe heatwaves (Perkins et al., 2012) in the last few decades. Precipitation changes have been spatially variable over eastern Australia, where northern Queensland grew wetter and southeast Australia grew drier (Fig. A2). In the last two decades, southeast Australia experienced the two worst droughts in the observational record (2001van Dijk et al. (2013van Dijk et al. ( ) and 2017van Dijk et al. ( -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), 35 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), from observations (Zhu et al., 2016;Walker et al., 2020), or to separate 40 from other drivers (e.g. climate variability, disturbances, and changes in land management). 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 resulted from increased atmospheric CO 2 concentrations (Donohue et al., 2009;Ukkola et al., 2016). The evidence for increases in leaf area from rising CO 2 have also been supported by observations of reduced runoff in Australia's drainage 45 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 of co-occurring changes in climate, land-use, and disturbance are confounded with 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 50 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 often-used 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 55 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 60 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 of 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 can greening trends be explained by rising CO 2 ? Using eastern Australia as a model system, we used a multi satellite derived NDVI record encompassing 38 years 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-useefficiency (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 two decades.

Study area
The study region encompasses the dominant woody ecosystems of eastern Australia (Fig. A1b). We used the National Vegetation Information System 5.1 land cover dataset (Dep ; Table A1) to select locations designated as "Acacia Forests and Woodlands", "Acacia Open Woodlands", "Callitris Forests and Woodlands", "Casuarina Forests and Woodlands", "Eucalypt 75 Low Open Forests", "Eucalypt Open Forests", "Eucalypt Open Woodlands", "Eucalypt Tall Open Forests", "Eucalypt Woodlands", "Low Closed Forests and Tall Closed Shrublands", "Mallee Open Woodlands and Sparse Mallee Shrublands", "Mallee Woodlands and Shrublands", "Melaleuca Forests and Woodlands", "Other Forests and Woodlands", "Other Open Woodlands", "Rainforests and Vine Thickets", and "Tropical Eucalypt Woodlands/Grasslands". 80 We used the atmospheric CO 2 record from the deseasonalized Mauna Loa record (https://www.esrl.noaa.gov/gmd/ccgg/trends/ data.html), 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 hours. PET was calculated from shortwave radiation and mean air temperature using 85 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 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 90 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-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 (Organization, 2017). We used this climatology to define the mean annual P:PET (MI MA ), and as the reference to calculate a 95 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).

Climate and remote sensing datasets
We used surface reflectance from two satellite products to generate the NDVI record: National Oceanic and Atmospheric Ad-  (Table A1; Schaaf & Wang, 2015). NDVI data was 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 105 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 less than three daily retrievals were removed. They were also removed when the coefficient of variation of 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. 110 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 & 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 115 that experienced bushfires from the year 2001 onwards. Specifically, these pixels were masked for the year of burning and the following three 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/2020. Both NDVI AVHRR and NDVI MODIS datasets were processed using Google Earth Engine (Gorelick et al., 2017), and exported at 5 km spatial processing used the 'stars' (Pebesma, 2020) and 'data.table' (Dowle and Srinivasan, 2019) R packages (see code 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 one million observations from the overlapping 2001-2016 portion of both records to fit a GAM using the 'mgcv' R package (Wood, 2017) to model NDVI MODIS from AVHRR derived covariates where 's' represents a penalized smoothing function using a thin-plate regression spline, SZA is the solar zenith angle, N DV I AV HRR is the uncalibrated NDVI from AVHRR, T OD 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.

Estimating NDVI and climate trends
We estimated the relative increase in NDVI between 1982-2019 with respect to time (equation 2) for each grid cell with an 135 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 140 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 Consortium, 2019). The temporal covariate was recentered to start with the first hydrological year (where the year starts one 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 used the merged NDVI observations to fit multiple statistical models to quantify the impact of changes in CO 2 and meteorological variables on NDVI. The relationship between NDVI and the running 12-month mean of P:PET was strongly nonlinear and followed a monotonic saturating sigmoidal relationship as indicated by GAM fits (methods equation 6, see below). 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 used nonlinear least squares (nls.multstart package (Padfield and Matheson, 2020) Fig. A6). The focus 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 (ppm) 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 platform. 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 one million observations per model fit. Corresponding goodness-of fit metrics (R 2 and root mean square error) were calculated by season ( Fig. 5) with one million randomly sampled observations. Alternative nonlinear functional forms were also fit to characterize the effect of CO 2 165 upon NDVI. A logistic model was fit across space for each hydrological year as where N DV I 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 M I f.anom .
Each seasonal model component was fit across space with one million random samples from the total merged NDVI record (approx 14.3 million observations).

175
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 (equation 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 MI anom is the annual anomaly of P:PET, 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 before regression. The standardized CO 2 and P:PET anom .
Next, 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, 185 VPD, and MI as fractions of their mean annual values (MA) as follows.
Finally we estimated the CO 2 effect across the study region using a GAM with a penalized smoothing function (s) characterizing the effect of the anomalies and mean annual values of VPD, P, and PET and 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 assimilation (umol m 2 s −1 ), E is leaf level transpiration (mol m 2 s −1 ), C a is atmospheric CO 2 (umol umol −1 ),

195
C i is intercellular CO 2 (umol umol −1 ), χ is Ci Ca , 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: 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 200 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 into equation (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: The change in E canopy can then be defined as: If we assume there is no long-term overall change in precipitation then we can assume change in E canopy is tightly coupled to the water supply, therefore we have: NDVI is linearly related to foliar cover (F ) until LAI ≈ 3 (m 2 m −2 ) (Carlson and Ripley, 1997), which is the predominantly the case when P:PET < 1. Most woody ecosystems of eastern Australia are strongly water limited with LAI ≤ 1 (m 2 m −2 ), 215 where NDVI is approximately proportional with the fraction of foliar cover: Then substituting equation (15) into equation (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 220 can predict the change towards NDVI to be We compared the WUE theoretical model with the robust linear models fit for each pixel location (equation 8), and the GAM (equation 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 225 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 230 the PETA hypothesis as another point of comparison with the CO 2 attributable effect on NDVI. and NDVI for the period of 1982-1986 are shown as histograms above and to the right of the main panel. Note that the majority of arrows shift towards higher aridity (lower P:PET) and higher 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 250 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-7). 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 (equation 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 (equation 7; Fig. A4). The GAM 255 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 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 wangRecentGlobalDecline2020 over, we found no consistent evidence of a

CO 2 driven greening and expectations from water use efficiency
The range of statistically estimated CO 2 -attributable greening responses was compared with the expectation from the theo-270 retical CO 2 water use efficiency model. The Donohue et al. (2013) CO 2 x 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% (10-90% 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 simultaneous increase in VPD (which the WUE ac-275 counts 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. 7d). 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 280 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 the Arid region also experienced the greatest increase in VPD and reduction in P over the 38 year period (Fig. 7a, A8).
In contrast, the largest estimated CO 2 -attributable effect on greening were found to be in the Equatorial and Tropical areas. The largest decline in tree cover occurred in the Temperate regions (Fig. 8c,A10. Most regions experienced declines in non-290 vegetated (bare) cover, increases in non-tree vegetation, and modest change in tree-cover (Fig. 8c), however the proportional increase of non-tree vegetation typically exceeded tree cover increases (Fig. A10). The relative expected increase in NDVI following the theoretical WUE prediction where 50% of the gain is allocated to foliar area (methods -eq 11).

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 295 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 300 New South Wales . 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), yet a more expansive global-scale analysis of all terrestrial vegetation using dynamic global vegetation model attribution did not connect Australia's greening to changes in atmospheric CO 2 concentration (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 305 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 of NDVI (Fig. 4), the nearly four decade long record allowed us to separate 310 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 nine 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. First, while the relative increase in NDVI between the 1982-2000 AVHRR epoch ( 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 320 NDVI increases from the WUE model (Fig. 7d).
It should be noted not all regions experienced consistent greening trends throughout the observation period. For example, 'greening' shifted to 'browning' during the austral summer (Dec -Feb) 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 summer time 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.). Further, greening may have been suppressed in parts of Queensland because cattle ranching activity has intensified and has driven forest conversion to managed pasture in the region (McAlpine et al., 2009).

Attributing a CO 2 fertilization contribution towards greening 330
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 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. 335 5-7;A4-7). 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 temporally weakening CO 2 effect (Fig. 6). The WUE model predicted similar relative rates of NDVI increase between the two epochs

Deviations from WUE
We found the CO 2 effect on foliar area was the smallest in the driest climate regions . This was at odds with the WUE prediction (at 50% allocation; Fig. 7e), and contrary to 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 the phenology of the vegetation composition, and disturbances not 350 captured by the satellite products (e.g. small fires, grazing). Browning in the Grassland region (Fig. 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, 355 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 2 -attributed 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.

360
The Arid region experienced the greatest relative increase in VPD (Fig. 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 experimentation from the Nevada Desert FACE experiment (Smith et al., 2014).

Relation to ecosystem CO 2 fertilization experiments
Notably, a four 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 370 than the elevated CO 2 exposure time in the experiment (four 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 (Drake et al., 2016;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 375 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 380 non-tree vegetation cover (e.g. grasses/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 veracity of the MODIS Vegetation Continuous Fraction product (DiMiceli et al., 2017) to accurately distinguish Australian tree cover from non-tree vegetation, however this pattern is consistent with a recent LiDAR derived treecover time series of Australia (Liao et al., 2020). This suggests the drought starting in 2017 was already killing trees prior to 385 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 (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 southeast Australian forests and woodlands (De Kauwe et al., 2020). 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 increasing CO 2 produced nearly four decades of widespread vegetation greening across eastern Australia. The large agreement between a theoretical 395 model and the statistically estimated CO 2 effect indicated that greening resulted through 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 400 further increases of atmospheric CO 2 will continue to enable vegetation to mitigate increases of aridity and VPD under future warming. It is also unclear if trees or grasses are the primary contributors to the recent greening trend. Future localised 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 405 considered in planning for future land and water resource management on the world's driest inhabited continent.
Code and data availability. All data used are publicly available from sources listed in Table A1. Processed data used in model fitting and figures can be accessed via Zenodo data repository: [future repository link]. 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.  N DV I3mo = s(CO2) + s( P P ETMA ) + s(CO2, P P ETMA ) + s( P P ETanom. )

Appendix A: Figures and tables in appendices
Here 's' represents a thinplate spline smoothing function, P P ET M A is the 30-year mean annual P:PET, and P P ETanom is the annual anomaly of P P ET .