Low and contrasting impacts of vegetation CO2 fertilization on terrestrial runoff over the past three decades: Accounting for above- and below-ground vegetation-CO2 effects

1 State Key Laboratory of Hydroscience and Engineering, Department of Hydraulic Engineering, Tsinghua University, Beijing, China 2 CSIRO Land and Water, Black Mountain, Canberra, ACT 2601, Australia 10 3 Australian Research Council Centre of Excellence for Climate Extremes, The Australian National University, Canberra, Australia 4 Key Laboratory of Water Cycle and Related Land Surface Processes, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing, China 5 Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, Peking University, Beijing 15 100871, China.

Abstract. Elevation in atmospheric carbon dioxide concentration (eCO2) affects vegetation water use, with consequent impacts on terrestrial runoff (Q). However, the sign and magnitude of the eCO2 effect on Q is still contentious. This is partly due to the poor understanding of the opposing eCO2-induced water effects at different scales, being water-saving caused by partial stomatal closure at the leaf-level contrasting with increased water-consumption due to increase foliage cover at the canopy level, leading 25 to highly debated findings among existing studies. None of the existing studies implicitly account for eCO2-induced changes to below-ground vegetation functioning. Here we develop an analytical ecohydrological framework that includes the effects of eCO2 on plant leaf, canopy density, and rooting characteristics to attribute changes in Q and detect the eCO2 signal on Q over the past three decades.
Globally, we detect a very small decrease of Q induced by eCO2 during 1982-2010 (-1.69%). When 30 assessed locally, along the resource availability (α) gradient, a positive trend (p<0.01) in the Q-eCO2 response is found ranging from a negative response (i.e., eCO2 reduces Q) in low α regions (typically dry) to a positive response (i.e., eCO2 increases Q) in high α areas (typically warm and humid). Our findings suggest a minor role of eCO2 on changes in global Q over the past three decades, yet highlights the negative Q-eCO2 response in semi-arid and arid regions which may further reduce the limited water 35 resource there.

Introduction
Runoff (Q) is the flow of water over the Earth's surface, forming streamflow, becoming one of the most important water resources for irrigation, hydropower and other human needs (Oki and Kanae, 2006).
Anthropogenic climate change is expected to alter the global hydrological cycle, with greenhouse gas-40 induced climate warming intensifying the hydrological cycle (Huntington, 2006). Besides climate, terrestrial vegetation also affects the water cycle. It is well-documented that elevated atmospheric CO2 concentration (eCO2) reduces stomatal opening, which in turn suppresses leaf-level transpiration (Field et al., 1995;Donohue et al., 2013). If this were the only mechanism that eCO2 changed vegetation this would increase runoff (Q) (Gedney et al., 2006). However, eCO2 increases vegetation foliage cover 45 (Donohue et al., 2013;Zhu et al., 2016), leading to enhanced canopy-level transpiration and consequently reductions of Q (Piao et al., 2007). The two opposite responses of vegetation water use to eCO2 complicate the net effect of eCO2 on Q, and existing modeling results are highly debated since they focus on different aspects of how eCO2 affects the plants and thus the water cycle (Gedney et al., 2006;Piao et al., 2007;Huntington, 2008;Cheng et al., 2014;Yang et al., 2016a;Ukkola et al., 2016). 50 Moreover, those previous modelling results have not been thoroughly validated against observations. In addition to stomatal and above-ground vegetation structure responding to eCO2, the below-ground vegetation structure (i.e., rooting depth) is also affected by eCO2, with eCO2 increases rooting depth overwhelmingly found in observations (Nie et al., 2013) (Supplementary Tables S1 and S2). Deeper rooting depth means larger plant-available water storage capacity by allowing plants to access deeper 55 soil moisture, which potentially increases transpiration water loss and reduces Q, especially during dry spells (Trancoso et al., 2017;Yang et al., 2016b). To date, no previous eCO2-Q modeling attempts have explicitly considered the below-ground eCO2-induced feedback (with the above-ground feedbacks): this paper fills that niche.
Here we use a parsimonious, analytical eco-hydrological model based on the Budyko framework (i.e., 60 the Budyko-Choudury-Porporto, BCP model; Donohue et al., 2012), in combination with an analytical rooting depth model based on ecosystem optimality theory (Guswa, 2008), an analytical CO2 fertilization model for steady-state vegetation (Donohue et al., 2017) and observed plant stomatal response to eCO2 (Ainsworth and Rogers, 2007), to detect the impact of eCO2 on Q changes (dQ) over the past three decades (i.e., 1982-2010). The Budyko framework describes the steady-state (i.e., mean 65 annual scale) hydrological partitioning as a functional balance between atmospheric water supply (i.e., precipitation, P) and demand (i.e., potential evapotranspiration, EP) and a model parameter that modifies the climate-hydrology relationship (Choudhury, 1999;Donohue et al., 2012). In this framework, both EP and the land surface parameter are affected by the response of vegetation to eCO2 (see Methods). The developed framework allows analytical and transparent attribution of dQ changes, which overcomes the 70 uncertainty raised from non-linear interactions among numerous processes when attributing dQ numerically by using bottom-up earth system models (Yang et al., 2015). To examine the long-term eCO2 impact and to minimize year-to-year "transient" effects (i.e., water storage changes), we performed our analyses using sequential 5-year periods (Yang et al., 2016a;Han et al., 2020), resulting https://doi.org/10.5194/hess-2020-548 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License. 4 in six 5-year-means during 1982-2010, with the first period containing 4 years. Additionally, since 75 vegetation response to eCO2 can be greatly mediated by the availability of other resources (e.g., water, light and nutrients) (Donohue et al., 2013;Donohue et al., 2017;Nenami et al., 2003;Yang et al., 2016a;Norby et al., 2010), we examine the impact of eCO2 on Q along a resource availability gradient (Donohue et al., 2017;Friedkubgstein et al., 1999) (see Methods). The resource availability is typically low in dry environments and increases as the climate becomes more humid, which enable us to detect 80 the signal of eCO2 on Q across a drywet gradient.

Runoff simulation
The Budyko-Choudhury-Porporato (BCP) model was adopted here to simulate Q and to attribute changes in Q (Yang et al., 2016b;Donohue et al., 2012). Choudhury's formulation of the Budyko curve 85 is (Choudhury, 1999): where E is the actual evapotranspiration (mm yr -1 ). P is the precipitation depth (mm yr -1 ). EP is the potential evapotranspiration (mm yr -1 ) and is estimated here using the Shuttleworth-Wallace two-source evapotranspiration model (Shuttleworth and Wallace, 1985) with the assumption of full soil moisture 90 supply (by taking soil surface resistance equal to zero and stomatal resistance equal to its non-waterstressed value) while allowing leaf area and leaf-level conductance to vary with atmospheric CO2 concentration (Ca) (Milly and Dunne, 2016). A recent study by Milly and Dunne (2016)  dimensionless parameter, which is a function of effective rooting depth (Ze, mm), mean rainfall 100 intensity (β, mm per event) and soil water holding capacity (θ, mm 3 mm -3 ) and exhibits a close relationship with the Choudhury's parameter n (Yang et al., 2016b;Porporato et al., 2004). Taking data from Porporato et al. (2004), we deduced the relationship between n and ω as (R 2 =0.96, p<0.001; Supplementary Figure S1): Effective rooting depth (Ze) was determined using an analytical carbon cost-benefit model based on ecosystem optimality theory proposed by Guswa (2008). The Ze model is given as (Guswa, 2008): where W is the ratio of mean annual P over potential transpiration, EP_T. γr is the root respiration rate (g C g -1 roots day -1 ), which is quantified using the Q10 equation (Lloyd and Taylor, 1994). RLD is the root length density (cm roots cm -3 soil) and SRL is the specific root length (cm roots g -1 roots). We fixed RLD to be 0.1 cm roots cm -3 soil and SRL to be 1500 cm roots g -1 , representing the median value of these two parameters reported in the literature, respectively (Caldwell, 1994;Eissenstat, 1997;Fitter 115 and Hay, 2002;Pregitzer et al., 2002). fGS is the fraction of growing season within a year, with the growing season length quantified according to Zhu et al. (2016). WUE is the photosynthetic water use efficiency (g C cm -3 H2O), which is determined for the first period (i.e., 1982-1985) from the ensemble means from eight Earth system models (described later) of annual gross primary production (GPP) and ET estimates (i.e., WUE=GPP/ET). For the following periods, WUE was estimated by considering the 120 https://doi.org/10.5194/hess-2020-548 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License.
6 effects of changes in Ca and vapor pressure deficit (v) on WUE (Donohue et al., 2013;Wong et al., 1979;Farquhar et al., 1993): where t is time in year. Note that the above equation implicitly assumes the same upscaling factor when converting the leaf-level assimilation and transpiration to the canopy level for a given location 125 (Donohue et al., 2017). The spatial pattern of mean annual Ze is shown in Supplementary Figure S2.

Attribution runoff changes
We used the BCP model to attribute changes in Q (dQ) due to different influencing factors following the method developed by Roderick and Farquhar (2011). To first order, change in Q (dQ) is: where ⁄ , P ⁄ and ⁄ represent the sensitivity of Q to changes in P, Ep and n, respectively, and can be expressed as: The physiological (stomatal conductance, Cs) and structural (Leaf area index, L, and effective rooting depth, Ze) impact both EP and n. More specifically, decreases in Cs lower the transpiration rate per leaf area, whereas increases in L and Ze enhance the canopy level transpiration rate. Additionally, increases in L also reduce soil evaporation by shading the soil surface (Shuttleworth and Wallace, 1985). The impact of eCO2 on parameter n is expressed through its impact on Ze. On one hand, increases in WUE 140 induced by eCO2 permit a larger vegetation carbon uptake per amount of water loss, potentially leading https://doi.org/10.5194/hess-2020-548 Preprint. Discussion started: 14 December 2020 c Author(s) 2020. CC BY 4.0 License. 7 to more carbon allocated to roots and thus a deeper Ze. Conversely, increases in plant water demand (as quantified by potential transpiration) would require plants to develop a deeper root to access to soil moisture at deeper depths, and vice versa (Guswa, 2008). As a result, we write EP and Ze as: and changes in EP and Ze are given by: where EP_M is the meteorological component of EP (without considering Ca). O represents factors other 150 than eCO2 that affects Ze, which effectively encodes the climate change-induced vegetation change.
Combining Eqs. (2), (7), (13) and (14) The first term on the right hand of Eq. (15) represents dQ caused by P change, the second term represents dQ caused by eCO2 and the third term calculates dQ induced by changes in EP_M. To 155 maintain simplicity, we calculate EP_M-induced dQ by subtracting the effect of eCO2-caused changes in EP on Q from dQ caused by changes in EP (i.e., P P − P P a a ). The fourth and fifth terms on the right hand of Eq. (15) represent dQ caused by changes in rainfall intensity and climate changeinduced vegetation change, respectively, and we group them as one factor in the attribution of dQ (i.e., other factors in Fig. 3). Since our primary focus was to examine how eCO2 affect Q and its relative The sensitivities of EP and Ze to eCO2 (i.e., P a and e a ) are quantified by numerically running the EP 165 model and Ze model with and without changes in Ca, respectively. The difference between the two simulations under the two Ca scenarios is considered the net effect of eCO2.

Stomatal conductance response to eCO2
The response of leaf-level stomatal conductance (Cs) response to eCO2 was determined using 244 field observations across a broad range of bioclimates (Ainsworth and Rogers, 2007). We linearly rescaled 170 the reported change in Cs for the magnitude of eCO2 in each study to obtain the sensitivity of Cs to eCO2: that is, the percentage change in Cs per 1% increase in Ca. We then classified the 244 observations based on their biome type to construct a biome type-based look-up table of Cs sensitivity to eCO2.

Resource availability index and leaf area index response to eCO2
175 The response of L to eCO2 was predicted based on the response of WUE to eCO2 adjusted by the local resource availability. We define the site resource availability index (α) based on growing season mean L following Donohue et al. (2017). This is because that observed L at a site is the net response to the local growing conditions and provides an effective proxy of the growing conditions experienced by plants (Donohue et al., 2017). Another advantage of this approach is that L can be readily measured directly or 180 remotely. We calculated α as, where τ is an exponential extinction coefficient, which typically varies from 0.3 to 1.2 (Campbell and Norman, 1998) and is set to be 0.7 herein. Broadly across the globe, α also corresponds well with climate aridity. The calculated α increases from 0.0 with low resource availability (typically dry) to 1.0 185 with high resource availability (typically warm and humid) (Figure 1). This suggests a predominant role of the climate in shaping the global vegetation pattern (Nemani et al., 2003). This also implies that the resource limitations on plant growth are mainly exerted by climate, consistent with the framework of climate limitation on vegetation proposed in previous studies (Nemani et al., 2003;Budyko, 1974;Yang et al., 2015). Then following Norby and Zak (2011), who showed that the observed response of L to 190 eCO2 was a non-linear function of L, we estimated the relative change in L induced by eCO2 per Donohue et al., (2017):

Data
To focus on the impacts of eCO2 on Q via feedbacks through vegetation and to eliminate potential 195 human impacts on Q, we limit our analyses to 2,268 strictly selected unimpaired catchments located across the globe ( Figure 2 than 10% of the catchment mean annual Q). Exactly 2,268 catchments pass the above selection criteria based on AVHRR GIMMS-3g NDVI data (Pinzon and Tucker, 2014). Land cover classification was acquired from the Moderate Resolution Imaging Spectroradiometer (MODIS) land use map (MOD12Q1) available from the NASA data center (Friedl et al., 2010). The global C4 vegetation fraction was obtained from the NASA data center 220 (http://webmap.ornl.gov/ogcdown/dataset.jsp?ds_id=932). Soil texture data at 30'' spatial resolution was acquired from the Harmonized World Soil Database (HWSD) (Nachtergaele, 2009), which was used to determine the value of θ according to the US Department of Agriculture (USDA) soil classification (Saxton and Rawls, 2006). These gridded data were further aggregated for individual catchments at a mean annual scale (i.e., five years).

225
3 Results and discussion

Validation of the BCP model in runoff estimation
The validity of the BCP model is tested by comparing the estimated Q with observed Q, in terms of both spatial and temporal variability, at the 2,268 unimpaired catchments (Figure 3). Spatially, the BCP model well captures the observed spatial variability in Q at the mean annual scale, with a coefficient of 230 determination (R 2 ) of 0.93, root-mean-squared error (RMSE) of 87.9 mm yr -1 and mean bias (estimated Q minus observed Q) of -11.4 mm yr -1 (Figure 3a) Temporally, trends in mean annual Q are also reasonably reproduced by the BCP model, which produces an R 2 of 0.71, RMSE of 0.71 mm yr -2 and mean bias of -0.05 mm yr -2 ( Figure 3b) In addition, we also perform a sensitivity analysis by comparing the simulated Q using the BCP model with and without considering eCO2. Results show that the BCP 235 model, when considering eCO2, performed differentially better in estimating Q trends than the BCP model without considering eCO2, suggesting that the developed analytical framework herein can well capture the eCO2 signal on the observed Q changes (Figure 3d).

Plant physiological and structural responses to eCO2
The physiological response of plant to eCO2, that is, the response of stomatal conductance (Cs) to eCO2 240 is directly compiled from field experiments and summarized for each plant functional type in Ainsworth and Rogers (2007) (also see Supplementary Figure S3). All those field experiments report a reduction of Cs in response to eCO2, with the largest (lowest) Cs reduction found in C4 crops (shrubs) for the same level of eCO2. On average, for a 1% increase in atmospheric CO2 concentration (Ca), Cs decreases by 0.47% ± 0.12% (mean ± one standard deviation), which means that Cs decreases by 5.67% ± 1. In terms of Ze, elevated Ca over the past three decades has resulted in a very minor (~1%) overall increase of Ze averaged across the globe (Figure 4e). Since large-scale observations of Ze in response to eCO2 are not available, we are not able to quantitatively validate the estimated response of Ze to eCO2. Nevertheless, the modelled result that eCO2 increases Ze is overwhelmingly found in site-and/or plant-260 level observations (Nie et al., 2013) (Supplementary Tables S1 and S2). Moreover, similar with L, the response of Ze to eCO2 also exhibits a notable difference along the resource availability gradient (Figure   4d and 4e). The positive response of Ze to eCO2 is larger in low α regions and gradually decreases as the resource availability becomes higher. In high α regions (e.g., tropical rainforest and southeast Asia), Ze even shows a slightly decreases in response to eCO2, suggesting a reduced plant water need in a high Ca 265 atmosphere in those regions.

Attribution of runoff changes over 1982-2010
During the last three decades, atmospheric CO2 concentration (Ca) increased by ~12.1%. For the same period, the BCP model detected a very small reduction in Q of ~0.73% (or 2.8 mm yr -1 ) induced by eCO2 across the 2,268 studied catchments ( Figure 5). The overall negative effect of eCO2 on Q suggests 270 that the structural forcing of eCO2 on vegetation water consumption (both above-and below-ground) outweighs the physiological effect of eCO2 driving leaf-level water saving. Despite the overall small effect averaged across all catchments, a significant positive trend (p<0.01) in the Q-eCO2 response is found along the resource availability gradient, from a negative response in low α catchments to a positive response in high α catchments ( Figure 5). 275 We then attribute dQ to different factors between 1982-1985 and 2006-2010 for the study catchments ( Figure 6). Compared with the early 1980s, mean observed Q over the 2,268 catchments in the late 2010s decreased by ~5.8 mm yr -1 , and the observed pattern with comparable magnitude in dQ is well captured by the BCP model ( Figure 6). The impact of eCO2 on dQ is estimated to be -2.3 mm yr -1 averaged over all 2,268 catchments. Consistent with relative Q changes (in %; Figure 5), the impacts 280 eCO2 on the absolute Q change (in mm yr -1 ) also exhibit significant upward trend as α increases (~0.97 mm yr -1 per 0.1 increase in α, p<0.01). Compared to that, decreases in P led to a 2.7 mm yr -1 decreases in Q, and enhanced EP has resulted in a decreased Q by 1.6 mm yr -1 (Figure 6a). The comparable magnitudes of dQ induced by dP and eCO2 only exist when averaged across all 2,268 catchments, while for each resource availability category, the impact of P on Q generally dominates dQ and is often much 285 higher than that of eCO2 ( Figure 7). As for the impact EP on Q, it also shows a notable gradient with changes in α as detected for the eCO2 effect, with the impact of EP on Q being increasingly negative as α raises (Figure 6b-f). Other factors include changes in rainfall intensity (Porporato et al., 2004) and climate change-induced vegetation change (e.g., higher L) have, in general, exerted a small negative impact on Q. 290 The same conclusions that the impacts of eCO2 on vegetation have limited yet contrasting (between warm-humid, high α regions and dry, low α regions) feedbacks on Q retain beyond the 2,268 catchments (Figure 7a and b Figures 7a and b). This is because in extremely dry areas, the availability of water defines the outcome and the sensitivity of Q to any changes in land surface properties is very small (Donohue et al., 2013;Roderick et al., 2014). The negative SQ_to_eCO2 diminishes quickly as α increases and turns to be a positive SQ_to_eCO2 in high α 300 regions. The overall small SQ_to_eCO2 is further manifested when comparing SQ_to_eCO2 with the sensitivities of Q to climate variables (i.e., P and EP). Averaged over the globe, a same relative change in P and EP would respectively lead to a ~10-times and ~4-times stronger impact on Q than eCO2 does, highlighting a predominant role of climate in shaping the global Q regime (Figure 7c-f and Supplementary Figure S4).

Discussion and concluding remarks
Elevation in atmospheric CO2 concentration is regarded as the ultimate driver of anthropogenic climate change, with consequent impacts on terrestrial runoff. Although the impacts of climate change on Q has been extensively documented in previous studies, the response of Q to eCO2 through vegetation feedbacks is less understood and remains controversial in existing studies (Gedney et al., 2006;Piao et 310 al., 2007;Huntington, 2008;Cheng et al., 2014;Yang et al., 2016a;Ukkola et al., 2016). Here, by developing an analytical attribution framework, we detected a very small response of global Q to eCO2induced changes in vegetation functioning ( Figure 5-7), suggesting that the eCO2 vegetation feedback only exert a minor impact on water resources for the range of eCO2 that we have experienced over the past three decades. 315 We also detected a significant positive trend (p<0.01) in the Q-eCO2 response along the resource availability gradient (Figure 5-7), which is consistent with field experiments (Norby and Zak, 2011;De Kauwe et al., 2013;Körner and Arnone, 1992) (Figure 4c), satellite observations (Donohue et al., 2013;Zhu et al., 2016;Yang et al., 2016a), and model attributions (Cheng et al., 2014;Lian et al., 2018). This Q-eCO2 response mechanism suggests that the structural response of vegetation to eCO2 (i.e., increases 320 in L and Ze) is larger in areas with lower resource availability, and gradually decreases as resources become less limiting on plant growth (Figure 4). The positive response of Q to eCO2 in high α catchments (primarily located in tropical rainforests) implies a dominant effect of eCO2-induced partial stomatal closure over increases in L and Ze on E in these environments. This is reasonable, as both theoretical predictions and in-situ observations have consistently reported a negligible response of L to 325 eCO2 in humid and close-canopy environments (Donohue et al., 2017;Yang et al., 2016a;Norby and Zak, 2011;Körner and Arnone, 1992). In such environments, water is generally abundant with light and/or nutrient availability being the most limiting resources for vegetation growth (Nemani et al., 2003;Yang et al., 2015), and plants have evolved to efficiently capture light by maximizing their above-ground structure (i.e., L). As a result, in these tropical rainforests plants have already absorbed 330 most of the incident light and any extra leaves would not materially increase the light absorption (Yang et al., 2016a).
Our findings have important implications for improved understanding of the global hydrological cycle and managing the world's water resources in a changing climate. Climate models have predicted an increased Q that is primarily driven by an increased P for the 21 st century (Milly and Dunne, 2016;335 Swann et al., 2016;Yang et al., 2018). Here we show that eCO2 would mitigate this positive impact of climate change on Q in relatively dry regions but exaggerate the Q increase in relatively wet regions via its impacts on vegetation water use. In addition, higher Ca and increased P enhance the availability of resource for vegetation growth, which increases vegetation coverage or L (Piao et al., 2020;Zhang et al., 2020a;Zhang et al., 2020b). This suggests that the structural response of vegetation to eCO2 may eventually decrease and the physiological effect of vegetation to eCO2 may become increasingly dominant in the overall response of vegetation water use to eCO2, leading to an increasing water-saving effect of plant in response to eCO2 under future climate change (Zhang et al., 2020b). In fact, analyses of the state-of-the-art climate model outputs have already consistently shown this water-saving effect of eCO2 at the global scale and especially in relatively warm and humid environments where L is high 345 (Yang et al., 2019). Yet, the impacts of eCO2 on Q in relatively dry regions are still highly uncertain and show a great diversity between climate models (Zhang et al., 2020b). In this light, our findings based on the well-validated analytical framework provide insightful guidance for climate model development that improves the models' capability in representing the vegetation and hydrological responses to eCO2.

350
All data for this paper are properly cited and referred to in the reference list.

Author contribution
YY and TRM designed the study. YY performed the calculation and drafted the manuscript. TRM, DY, YZ, SP, SP, and HEB contributed to results discussion and manuscript writing.

355
The authors declare that they have no conflict of interest. States

Figure 7
Sensitivity of Q to eCO2 and its relative importance to P and EP across the globe. a, Spatial distribution of Q sensitivity to eCO2 (% change in Q per 1% change in Ca). b, Boxplot of Q sensitivity to eCO2 for each resource availability category. c, Relative importance of eCO2 on Q compared to changes in P on Q (% change in Q per 1% change in Ca compared to % change in Q per 1% change in P). d, Boxplot of the relative importance of eCO2 on Q compared to changes in P on Q for each resource availability category. e, Relative importance of eCO2 on Q compared to changes in EP on Q (% 545 change in Q per 1% change in Ca compared to % change in Q per 1% change in EP). f, Boxplot of the relative importance of eCO2 on Q compared to changes in EP on Q for each resource availability category. In b, d and f the upper / lower box edges represent the quantile divisions, the inner horizontal line is the median, the dots indicate the mean value, and the dashed line represent the 1% and 99% percentiles.   Figure   590 3(d).