Low and contrasting impacts of vegetation CO2 fertilization on global terrestrial runoff over 1982–2010: accounting for aboveground and belowground vegetation–CO2 effects

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 are still
contentious. This is partly due to eCO2-induced changes in vegetation
water use having opposing responses at the leaf scale (i.e., water-saving effect caused by
partially stomatal closure) and the canopy scale (i.e., water-consuming
induced by foliage cover increase), leading to highly debated conclusions
among existing studies. In addition, none of the existing studies explicitly
account for eCO2-induced changes to plant rooting depth that is
overwhelmingly found in experimental observations. 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 to detect the eCO2 signal on Q via vegetation feedbacks over
1982–2010. Globally, we detect a very small decrease of Q induced by
eCO2 during 1982–2010 (−1.7 %). Locally, we find a small positive
trend (p < 0.01) in the Q–eCO2 response along a resource
availability (β) gradient. Specifically, the Q–eCO2 response is
found to be negative (i.e., eCO2 reduces Q) in low-β regions
(typically dry and/or cold) and gradually changes to a small 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 1982–2010, yet we highlight that a negative Q–eCO2 response in
semiarid and arid regions may further reduce the limited water resource
there.



Introduction
Runoff (Q) is the flow of water over the Earth's surface, forming streamflow, and represents one of the most important water resources for irrigation, hydropower, and other human needs (Oki, 2006).Anthropogenic climate change is expected to alter the global hydrological cycle, with greenhouse-gas-induced climate warming intensifying the hydrological cycle (Huntington, 2006).Besides climate, terrestrial vegetation also affects the water cycle (Brown et al., 2005).It is well documented that elevated atmospheric CO 2 concentration (eCO 2 ) reduces stomatal opening, which in turn suppresses leaf-level transpiration (Field et al., 1995).If this were the only mechanism, i.e., that eCO 2 changed vege-Y.Yang et al.: Low and contrasting impacts of vegetation CO 2 fertilization tation, this would increase runoff (Q) (Gedney et al., 2006).However, eCO 2 increases vegetation foliage cover (Donohue et al., 2013;Zhu et al., 2016), leading to enhanced canopylevel transpiration and consequently reductions of Q (Piao et al., 2007).These two opposing responses of vegetation water use to eCO 2 complicate the landscape-scale net effect of eCO 2 on Q, and existing modeling results are highly debated since they focus on different aspects (i.e., physiological functioning and/or structural change) of how eCO 2 affects the plants and thus the water cycle (Fatichi et al., 2016;Gedney et al., 2006;Huntington, 2008;Piao et al., 2007;Yang et al., 2016a, b).Moreover, observational and evaluation studies of eCO 2 effects on Q remain limited, particularly at regional to global scales.
In addition to stomatal and aboveground vegetation structure responding to eCO 2 , the belowground vegetation structure (e.g., rooting depth) is also affected by eCO 2 , with eCO 2 increasing rooting depth overwhelmingly found in experimental observations (Nie et al., 2013) (Supplement Tables S1  and S2).Deeper rooting depth increases plant-available water storage capacity by allowing vegetation to access deeper 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 eCO 2 -Q modeling attempts have explicitly considered the belowground eCO 2 -induced feedback simultaneously with the two previously mentioned aboveground feedbacks; this paper fills that niche.
Here we use a parsimonious, analytical ecohydrological model based on the Budyko framework (i.e., the Budyko-Choudhury-Porporato (BCP) model; Donohue et al., 2012), in combination with an analytical rooting depth model based on ecosystem optimality theory (Guswa, 2008), an analytical CO 2 fertilization model for steady-state vegetation (Donohue et al., 2017), and observed plant stomatal response to eCO 2 (Ainsworth and Rogers, 2007), to detect the impact of eCO 2 on Q changes (dQ) via vegetation feedbacks over global vegetated lands for 1982-2010.The Budyko framework describes the steady-state (i.e., mean annual scale) hydrological partitioning as a functional balance between atmospheric water supply (i.e., precipitation, P ) and demand (i.e., potential evapotranspiration, E P ) and a model parameter that modifies the climate-hydrology relationship (Choudhury, 1999;Donohue et al., 2012).In this framework, both E P and the model parameter are affected by the response of vegetation to eCO 2 (see Methods section).The "top-down" (Sivapalan et al., 2003) developed framework allows for analytical and transparent attribution of dQ changes, which overcomes the uncertainty raised from nonlinear interactions among numerous processes when attributing dQ numerically by using "bottom-up" Earth system models (Yang et al., 2015).To examine the long-term eCO 2 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 in six 5-year means during 1982-2010, with the first period containing 4 years.Additionally, since vegetation response to eCO 2 can be greatly mediated by the availability of other resources (e.g., water, light, and nutrients) (Donohue et al., 2013(Donohue et al., , 2017;;Nemani et al., 2003;Yang et al., 2016a;Norby et al., 2010), we examine the impact of eCO 2 on Q along a resource availability gradient (Donohue et al., 2017;Friedlingstein et al., 1999) (see Methods section).Resource availability is typically low in dry (and/or cold) environments and increases as the climate becomes more humid, which enables us to detect the signal of eCO 2 on Q across a dry-wet gradient.
2 Material and methods

Methods
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).The BCP model uses the formulation by Choudhury (1999) of the Budyko curve to estimate Q (Eq. 1 below), in which the model parameter is estimated based on the relationship between Choudhury's model parameter and Porporato's model parameter (Eq. 2 below).The required rooting depth (Z r ) in estimating Porporato's parameter is calculated using the rooting depth model by Guswa (2008) (Eqs. 3-5 below).To quantify the response of Q to eCO 2 via vegetation feedbacks, the stomatal response of vegetation to eCO 2 is determined by upscaling the observed response at the site level to the biome level (Sect.2.1.4),and the leaf area index (L) response to eCO 2 is quantified based on the response of water use efficiency (WUE) to eCO 2 adjusted by the local resource availability following Donohue et al. (2017) (Sect. 2.1.5).The effects of eCO 2 on both stoma and L also affect rooting depth in the model by Guswa (2008).A flowchart of our modeling approach is summarized in Fig. 1, and detailed calculation procedures are described in Sect.2.1.1 to 2.1.5.

Runoff simulation
The BCP model adopts the formulation of the Budyko curve by Choudhury (1999), which is given as where E is the annual average actual evapotranspiration (mm yr −1 ).P is the annual average precipitation (mm yr −1 ).E P is the annual average potential evapotranspiration (mm yr −1 ) here estimated using the Shuttleworth-Wallace two-source evapotranspiration model (Shuttleworth and Wallace, 1985;see Sect. 2.1.2). n is a unitless model parameter that encodes all factors other than mean climate conditions and modifies the partitioning of P between E and Q.For assumed steady-state conditions, Q is calculated by subtracting E from P as a result of catchment water balance.
The probabilistic steady-state solution to the stochastic dynamic soil moisture model by Porporato (2004) shares a similar form with the Budyko curve (Porporato et al., 2004).Porporato's parameter ω is a dimensionless parameter, which is a function of effective rooting depth (Z r , mm), mean rainfall intensity (α, mm per event), and soil water holding capacity (WHC, mm 3 mm −3 ), and it exhibits a close relationship with Choudhury's parameter n (Yang et al., 2016b;Porporato et al., 2004).A relationship between Porporato's ω parameter and Choudhury's n parameter was built following three steps.Firstly, we obtained the numerical solution of Porporato's model for the corresponding E/P for every 0.1 increment in E P /P for six separate ω curves.Secondly, by numerically solving Choudhury's formulation of the Budyko curve, we determined the values of Choudhury's parameter (n) that correspond to the E/P values of each of the six ω curves.Thirdly and finally, we pooled all n-ω pairs together and deduced the relationship between n and ω (R 2 = 0.96, p < 0.001; Supplement Fig. S1) as Effective rooting depth (Z r ) was determined using an analytical carbon cost-benefit model based on ecosystem optimality theory proposed by Guswa (2008).The Z r model is given as where W is the ratio of the multiyear growing season mean P over potential transpiration, E P _T ; γ r is the root respiration rate (g C g −1 roots d −1 ), which is quantified using the standard Q 10 theory (Lloyd and Taylor, 1994;Ryan, 1991) with a fixed Q 10 coefficient of 2.0 (Zhao et al., 2011).The base respiration rate at 20 • C for each biome type is determined following Heinsch (2003).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 (Caldwell, 1994;Eissenstat, 1997;Fitter and Hay, 2002;Pregitzer et al., 2002).f GS is the fraction of the 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 H 2 O), which is determined for the first period (i.e., 1982-1985) from the ensemble mean from eight ecosystem models (see Data section) of annual gross primary production (GPP) and transpiration (E T ) estimates (i.e., WUE = GPP / E T ).For the following periods, WUE was estimated by considering the effects of changes in atmospheric CO 2 concentration (C a ) and vapor pressure deficit (v) on WUE (Donohue et al., 2013;Wong et al., 1979;Farquhar et al., 1993) as 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 (Donohue et al., 2017).The spatial pattern of mean annual Z r is shown in Supplement Fig. S2.

The Shuttleworth-Wallace model
The Shuttleworth-Wallace two-source evapotranspiration model (the S-W model) was used to estimate E P and its two components (potential transpiration, E P _T , and potential evaporation, E P_S ) (Shuttleworth and Wallace, 1985).The S-W model estimates evapotranspiration as https://doi.org/10.5194/hess-25-3411-2021Hydrol.Earth Syst.Sci., 25,[3411][3412][3413][3414][3415][3416][3417][3418][3419][3420][3421][3422][3423][3424][3425][3426][3427]2021 where λ is the latent heat for vaporization (MJ kg −1 ), is the gradient of the saturation vapor pressure with respect to temperature (kPa K −1 ), ρ is the air density (kg m −3 ), c p is the specific heat of air at constant pressure (MJ kg −1 K −1 ), γ is the psychrometric constant (kPa K −1 ).r a a , r c a , and r s a are the aerodynamic resistance (s m −1 ) to heat and vapor transfer between the canopy-air space and the atmosphere, between the leaf and the canopy-air space, and between the soil surface and the canopy-air space, respectively.These three aerodynamic resistance terms are estimated following Sánchez et al. (2008).r s s and r c s are soil surface resistance and stomatal resistance (the reciprocal of stomatal conductance), respectively. T estimate E P using the S-W model, r s s is set to zero, and r c s is set to its non-water-stressed value (Milly and Dunne, 2016).The non-water-stressed values of r c s for each biome type are provided in Mu et al. (2007).A is the available energy (equal to net radiation minus ground heat flux, W m −2 ), and A s is the available energy at the soil surface, which is estimated as a function of L following Beer's law (Campbell and Norman, 1998;Yang and Shang, 2013).As a result, A -A s is the available energy absorbed by the plant canopy. Te impacts of eCO 2 on E P and its two components are obtained by allowing L and r c s to vary with C a .Recently, Milly and Dunne (2016) showed that the S-W model could most satisfactorily reproduce evapotranspiration estimates under non-water-limited conditions from climate models under eCO 2 .

Attribution of runoff changes
We used the BCP model to attribute changes in Q (dQ) due to different influencing factors following Roderick and Farquhar (2011).To a first-order approximation, the change in where ∂Q/∂P , ∂Q/∂E P , and ∂Q/∂n represent the sensitivity of Q to changes in P , E P , and n, respectively, and can be expressed as The physiological (stomatal conductance, g s ) and structural (leaf area index, L, and effective rooting depth, Z r ) parameters impact both E P and n.More specifically, decreases in g s lower the transpiration rate per leaf area, whereas increases in L and Z r 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 eCO 2 on parameter n is expressed through its impact on Z r .On the one hand, increases in WUE induced by eCO 2 permit a larger vegetation carbon uptake per amount of water loss, potentially leading to more carbon allocated to roots and thus a deeper Z r .Conversely, increases in plant water demand (as quantified by potential transpiration) require vegetation to develop deeper roots to access deeper soil moisture and vice versa (Guswa, 2008).As a result, we write the functional dependencies of E P and Z r as where ∂C a dC a .The fourth and fifth terms on the right-hand side of Eq. ( 23) represent dQ caused by changes in rainfall intensity and climate-change-induced vegetation change, respectively, and we group them as one factor in the attribution of dQ.Since our primary focus was to examine how eCO 2 affects vegetation and the consequent impact on Q, as well as its relative importance to changes in P and E P _M , the other factors driving dQ were estimated as the residual of Eq. ( 23) (i.e., total dQ minus the sum of dQ induced by dP , dE P _M , and eCO 2 ).By introducing Eqs. ( 17) and (18) into Eq.( 23), the sensitivity of Q to eCO 2 (S Q_to_eCO2 , mm yr −1 ppm −1 ) is written as The sensitivities of E P and Z r to eCO 2 (i.e., ∂E P ∂C a and ∂Z r ∂C a ) are quantified by numerically running the E P model and Z r model with and without changes in C a , respectively.The difference between the two simulations under the two C a scenarios is considered the net effect of eCO 2 on Q.

Stomatal conductance response to eCO 2
The response of leaf-level stomatal conductance (g s ) response to eCO 2 was determined using 244 field experiments with artificially elevated CO 2 across a broad range of bioclimates (Ainsworth and Rogers, 2007).We linearly rescaled the reported change in g s for the magnitude of eCO 2 in each of the 244 studies to obtain the sensitivity of g s to eCO 2 , that is, the percentage change in g s per 1 % increase in C a .We then classified the 244 observations based on their biome type to construct a biome-type-based look-up table of g s sensitivity to eCO 2 .

Resource availability index and L response to eCO 2
The response of L to eCO 2 was predicted based on the response of WUE to eCO 2 adjusted by the local resource availability.We define a site resource availability index (β) based on growing season mean L following Donohue et al. (2017).This is because 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 vegetation (Donohue et al., 2017).Another advantage of this approach is that L can be readily measured directly or 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 and/or cold) to 1.0 with high resource availability (typically warm and humid) (Fig. 2).This suggests a predominant role of the climate in shaping the global vegetation pattern (Budyko, 1974;Nemani et al., 2003;Yang et al., 2015).This also implies that the resource limitations on plant growth are mainly exerted by climate, which is 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 eCO 2 was a nonlinear function of L, we estimated the relative change in L induced by eCO 2 as per Donohue et al. (2017):

Data
The BCP model is validated against observed Q in 2268 strictly selected unimpaired catchments located across the globe that cover a broad range of bioclimates (Fig. 3).Originally, daily and/or monthly Q observations were collected from more than 22 000 catchments globally (Beck et al., 2019).Three selection criteria were implemented to ensure that only catchments with continuous Q records that are negligibly affected by humans were used.First, catchments with > 5 % missing data during the entire study period  were removed.Linear interpolation was applied to fill the gaps in the remaining Q series.Second, catchments smaller than 100 km 2 were excluded.This is to ensure that at least one precipitation pixel (i.e., 0.1 • × 0.1 • or ∼ 100 km 2 ) is included for a catchment.Third, we excluded catchments where observed Q is likely to be affected by human interventions, including catchments with (i) significant forest gain or loss (> 2 % of the total catchment area) (Hansen et al., 2013), (ii) irrigated areas larger than 2 % (Siebert et al., 2005), (iii) urban areas (https://www.edenextdata.com/?q=content/esa-globcoverversion-23-2009-300m-resolution-land-cover-map-0, last access: 17 June 2019) larger than 2 %, and (iv) the presence of large dams (Lehner et al., 2011) (i.e., where the reservoir's capacity in a catchment is larger than 10 % of the catchment mean annual Q).Exactly 2268 catchments pass these selection criteria (Fig. 3).

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 2268 unimpaired catchments (Fig. 4).Spatially, the BCP model well captures the observed spatial variability in Q at the mean annual scale, with a coefficient of 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 (Fig. 4a).Temporally, trends in mean annual Q are also reasonably reproduced by the BCP model, having an R 2 of 0.71, RMSE of 0.71 mm yr −2 , and mean bias of −0.05 mm yr −2 (Fig. 4b).Additionally, we also perform a sensitivity analysis by comparing the simulated Q using the BCP model with and without considering eCO 2 .Results show that the BCP model, when considering eCO 2 , performed better in estimating Q trends than the BCP model without considering eCO 2 , as evidenced by an improvement of R 2 by 0.02, a reduction of RMSE by 0.03 mm yr −2 , and a decrease of mean bias by 0.11 mm yr −2 , averaged over all 2268 catchments (Fig. 4d).More apparent improvements of the BCP model performance with the consideration of eCO 2 are found in regions having a relatively higher resource availability index.For β of 0.4-0.6,0.6-0.8, and 0.8-1.0, the mean bias of simulated Q trends with eCO 2 is −0.02, 0.06, or −0.36 mm yr −2 but increased to 0.24, 0.20, or −0.53 mm yr −2 , respectively, when eCO 2 is not considered (Fig. 4d).These results suggest that the ana-lytical framework developed herein captures the eCO 2 signal on the observed Q changes.

Plant physiological and structural responses to eCO 2
The physiological response of plants to eCO 2 , that is, the response of g s to eCO 2 , is directly compiled from field experiments and summarized for each plant functional type in Ainsworth and Rogers (2007) (also see Supplement Fig. S3).All those field experiments report a reduction of g s in response to eCO 2 , with the largest g s reduction found in C 4 crops and lowest in shrubs for the same level of eCO 2 .On average, for a 1 % increase in C a , g s decreases by 0.47 % ± 0.12 % (mean ± 1 standard deviation), which means that g s decreases by 5.67 % ± 1.47 % under a 12.1 % increase in C a over 1982-2010(i.e., from ∼ 343.7 ppm in 1982-1985to 385.2 ppm in 2006-2010;Keeling et al., 2001).This result is consistent with a recent isotope-based study (i.e., ∼ 5 % reduction of g s during the past 3 decades; Frank et al., 2015).
For structural response, averaged across global vegetated lands, our model reveals that elevated C a has caused an increase of L by 2.12 % (0.14 % and 3.88 % for 5th and 95th percentiles) over 1982-2010 (Fig. 5a and b).Despite this relatively small fertilization effect of eCO 2 on L at the global scale, an evident gradient is found in the L-eCO 2 response: a larger eCO 2 -induced relative L increase is found in low resource availability regions (smaller β value in Fig. 2a) and vice versa (Fig. 5b).This modeled pattern of L-eCO 2 response agrees very well with observations at the Free-Air CO 2 Enrichment (FACE) observations (R 2 = 0.96, p < 0.01; Fig. 5c) and is also consistent with large-scale satellite-based observations (Donohue et al., 2013;Zhu et al., 2016;Yang et al., 2016a).
In terms of Z r , our modeling results show that elevated C a over 1982-2010 has resulted in a very minor (0.93 %, −0.12 % and 1.85 % for 5th and 95th percentiles) overall increase of Z r averaged across the globe (Fig. 5e).Since large-scale observations of Z r in response to eCO 2 are not available, we are not able to quantitatively validate the estimated response of Z r to eCO 2 .Nevertheless, the modeled result that eCO 2 increases Z r is overwhelmingly found in siteand/or plant-level experiments (Nie et al., 2013) (Tables S1  and S2).Moreover, similar to L, the response of Z r to eCO 2 also exhibits a notable difference along the resource availability gradient (Fig. 5d and e).The positive response of Z r to eCO 2 is larger in low-β regions and gradually decreases as the resource availability becomes higher.In high-β regions (e.g., tropical rainforest and Southeast Asia), Z r even shows a slight decrease in response to eCO 2 , suggesting a reduced plant water need given the range of C a over 1982-2010 in those regions.

Attribution of runoff changes over 1982-2010
Over 1982-2010, C a increased by ∼ 12.1 %.For the same period, the BCP model detected a very small reduction in Q of ∼ 1.7 % (or 2.2 mm yr −1 ) induced by eCO 2 via vegetation feedbacks across the entire global vegetated lands (Figs.6b  and 7d).This 1.7 % reduction in Q, under the context of 12.1 % increases in C a , demonstrates a muted response of Q to eCO 2 .In addition, the overall negative effect of eCO 2 on Q suggests that the structural forcing of eCO 2 on vegetation water consumption (both aboveground and belowground) outweighs the physiological effect of eCO 2 driving leaf-level water saving.Across the global vegetated lands and for the same period, the physiological response of vegetation to eCO 2 has led to an increased Q by 0.7 % (or 0.9 mm yr −1 ), with the simulated Q increases being increasingly larger as β increases (Fig. 6d).By contrast, the structural response of vegetation to eCO 2 has resulted in an overall Q reduction by 2.4 % (or 3.1 mm yr −1 ), with the decreases in Q being increasingly smaller as β increases (Fig. 6e).These two opposite responses of vegetation water use to eCO 2 along the resource availability gradient have led to a significant positive trend (p < 0.01) in the Q-eCO 2 response along the resource availability gradient, from a negative response in low-β landscapes to a positive response in high-β landscapes (Fig. 6b).Nevertheless, an exception is found in extreme arid zones (i.e., when β < 0.1; Fig. 6b).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).
We then attribute dQ to different forcing factors between 1982-1985 and 2006-2010 over the global vegetated lands (Figs. 7 and 8).Compared with the early 1980s (i.e., 1982-1985), mean observed Q over the global vegetated lands in the late 2010s (i.e., 2006-2010) increased by 29.7 mm yr −1 , and the observed pattern with comparable magnitude in dQ is well captured by the BCP model (Fig. 4b and d).Consistent with relative Q changes (in %; Fig. 6), the impacts of eCO 2 on the absolute Q change (in mm yr −1 ) also exhibit a significant upward trend as β increases (0.53 mm yr −1 per 0.1 increase in β, p < 0.01).Compared to that, increases in P led to a 43.9 mm yr −1 increase in Q, and enhanced E P _M has resulted in a decreased Q by 5.3 mm yr −1 (Fig. 7f).For the entire vegetated lands and each resource availability category, the impact of dP on Q generally dominates dQ and is often much higher than that of eCO 2 (Fig. 7).An exception is the low-β regions (β < 0.2), where the impact of eCO 2 on Q outweighs the impact of dP on Q (Fig. 8a).As for the impact E P _M on Q, it also shows a notable gradient with changes in β as detected for the eCO 2 effect, with the impact of E P _M on Q being increasingly negative as β increases (Fig. 8b-e).The combined influence of other factors including changes in rainfall intensity (Porporato et al., 2004;Westra et al., 2013) and climate-change-induced vegetation change (e.g., higher L) have, in general, exerted a negative impact on Q. https://doi.org/10.5194/hess-25-3411-2021 Hydrol.Earth Syst.Sci., 25, 3411-3427, 2021 Since changes in meteorological factors (P and E P _M ) are often considered to dominate changes in Q and have been extensively examined previously (e.g., Roderick and Farquhar, 2011;Yang et al., 2018;Zhang et al., 2018), we next examine the sensitivity of Q to eCO 2 (S Q_to_eCO2 ) and compare it with the sensitivity of Q to changes in P and E P _M .Because C a has different units from P and E P _M , we use relative units to better compare the three sensitivities (Fig. 9).Globally, an increase in C a by 1 % only leads to a decrease of Q by ∼ 0.14 % (equivalent to ∼ 1.7 % for the range of eCO 2 experienced over 1982-2010).Similar to the attribution results shown above (Fig. 6a and b), S Q_to_eCO2 is generally more negative in global arid ecosystems where β is low (Fig. 9a and b).The negative S Q_to_eCO2 diminishes quickly as β increases and becomes positive S Q_to_eCO2 in high-β regions.The overall small S Q_to_eCO2 is further manifested when comparing S Q_to_eCO2 with the sensitivities of Q to P and E P _M .Averaged across the global vegetated lands, the same relative change in P and E P would, respectively, lead to an ∼ 10 times and ∼ 4 times stronger impact on Q than eCO 2 does.This highlights the predominant role of climate in shaping the global Q regime (Figs.9c-f and S4).

Discussion and concluding remarks
Elevation in atmospheric CO 2 concentration (and other greenhouses gases) is regarded as the ultimate driver of anthropogenic climate change, with consequent impacts on Q.Although the impacts of climate change on Q have been extensively studied, the response of Q to eCO 2 through vegetation feedbacks is less understood and remains controversial (Gedney et al., 2006;Piao et al., 2007;Huntington, 2008;Cheng et al., 2014;Trancoso et al., 2017;Yang et al., 2016a;Ukkola et al., 2016a, b).Here, by developing an analytical attribution framework, we detected a very small response of global Q to eCO 2 -induced changes in vegetation structural (both aboveground and belowground) and physiological functioning (Figs.6-8), suggesting that the eCO 2 vegetation feedback only exerts a minor impact on water resources (partly due to the two opposing water effects between the structural and physiological responses to eCO 2 ) for the range of eCO 2 experienced over 1982-2010.
The overall negative impact of eCO 2 on Q detected herein suggests that increased vegetation water consumption driven by the structural response of vegetation (i.e., increases in L and Z r ) to eCO 2 outweighs the functional change of leaflevel water saving caused by the physiological effect of eCO 2 (i.e., decreases in g s ).This result is consistent with previous findings by Cheng et al. (2014), Trancoso et al. (2017), and and (e) the upper and lower box edges represent the quartile divisions, the inner horizontal line is the median, the dots indicate the mean, and the dashed lines represent the 5th and 95th percentiles.The number of grid cells in each resource availability index category is provided in Fig. 2. Ukkola et al. (2016a).In addition, we also detected a significant positive trend (p < 0.01) in the Q-eCO 2 response along the resource availability gradient (Figs.6-9).This Q-eCO 2 response pattern suggests that the structural response of vegetation (i.e., increases in L and Z r ) to eCO 2 is larger in areas with lower resource availability and gradually decreases as resources become less limiting on plant growth (Fig. 5).The positive response of Q to eCO 2 in high-β catchments (primarily located in tropical rainforests; Fig. 6a) implies a dominant effect of eCO 2 -induced partial stomatal closure over increases in L and Z r on E in these environments (Fig. 6).This is reasonable, as both theoretical predictions and in situ observations have consistently reported a negligible response of L to eCO 2 in humid and closed-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 limiting vegetation growth (Nemani et al., 2003;Yang et al., 2015), and vegetation has evolved to efficiently capture light by maximizing its aboveground structure (i.e., L).As a result, in these high-L regions, vegetation has already absorbed most of the incident light, and any extra leaves would not materially increase the light absorption (Yang et al., 2016a).By contrast, in dry regions, an eCO 2 -induced increase in vegetation water use efficiency (so less transpiration for the same amount of carbon assimilation at the leaf level) would lead to an increase in L that is directly proportional to an increase in water use efficiency which would increase canopy-level carbon fixation (Fig. 5b).This finding is consistent with satellite observations (Donohue et al., 2013) and in situ FACE experiments (Norby and Zak, 2011).https://doi.org/10.5194/hess-25-3411-2021 Hydrol.Earth Syst.Sci., 25, 3411-3427, 2021 Our findings have important implications for an 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 21st century (Lian et al., 2021;Milly and Dunne, 2016;Swann et al., 2016;Yang et al., 2018).Here we show that eCO 2 -induced vegetation feedbacks would mitigate this positive impact of climate change on Q in relatively dry regions and exacerbate the Q increase in relatively wet regions.In addition, higher C a and increased P enhance the availability of resources for vegetation growth, which increases vegetation coverage or L (Piao et al., 2020;Zhang et al., 2020a, b).As the vegetation aboveground structural responses to eCO 2 decrease with the in-crease of L, the predicted future L increases suggest that the structural response of vegetation to eCO 2 may eventually decrease, and the physiological effect of vegetation on eCO 2 may become increasingly dominant in the overall response of vegetation water use to eCO 2 , leading to an increasing watersaving effect of vegetation in response to eCO 2 under future climate change (Zhang et al., 2020b).Analyses of the stateof-the-art climate model outputs already consistently show this water-saving effect of eCO 2 globally, especially in relatively warm and humid environments where L is high (Yang et al., 2019).Nevertheless, this may partly be because only some climate models consider the physiological effect while ignoring structural responses of vegetation to eCO 2 .In addition, the impacts of eCO 2 on Q in relatively dry regions are still highly uncertain and show a great diversity between climate models (Zhang et al., 2020b).
Finally, it is worthwhile noting that there are several limitations in the developed modeling framework.First, the rooting depth model by Guswa (2008) adopted herein employs an intensive root water uptake strategy, which assumes that root water uptake occurs at a potential rate (i.e., E P _T ) until soil moisture reaches the wilting point when transpiration is completely suppressed (Guswa, 2008).This intensive root water uptake strategy differs from the root water uptake strategy employed in the stochastic soil water balance model by Porporato et al. (2004), which is a more conservative strategy under which root water uptake linearly decreases with the decrease of soil moisture (Porporato et al., 2004).Combining the two strategies into one modeling framework potentially leads to inconsistency in the theoretical aspect of the approach.In fact, a later study by Guswa (2010) incorporated the soil water balance model by Porporato et al. (2004) into Guswa's cost-benefit framework for rooting depth (re-ferred to as the Guswa-2010 approach herein).However, the Guswa-2010 approach could not provide an explicit solution for Z r , because the solution for transpiration in Porporato's model is an incomplete gamma function of Z r (Guswa, 2010;Porporato et al., 2004).As a result, to allow an analytical solution to be derived, we used Guswa (2008) for Z r in our modeling framework.According to Guswa (2010), using the conservative root water uptake strategy resulted in a slightly deeper Z r compared to when the intensive strategy was used.Despite that, the response of Z r to changes in C a under the two strategies should be similar, as the effects of eCO 2 on Z r are expressed via water use efficiency and E P _T in our parameterization, which are independent of Z r parameterizations.This means that adopting different root water uptake strategies would only lead to differences in the resultant absolute magnitude of runoff (Q) but unlikely to result in differences in the response of Q to eCO 2 , especially when the relative magnitude is used (Figs.5d, e and 6a, b, e and f).Alternatively, Rodríguez-Iturbe and Porporato (2004) incorpo- Values in parentheses represent 1 standard deviation of each response among grid cells within each resource availability index category.The number of grid cells in each resource availability index category is provided in Fig. 2. rated the intensive root water uptake strategy into a stochastic soil water balance model and obtained a steady-state solution that has a simper form than the one by Porporato et al. (2004) and also mimics the Budyko curve.This approach deserves further investigation.Second, if interpreted strictly from a theoretical perspective, the model by Porporato et al. (2004) is more suitable to estimate hydrological partitioning during growing seasons instead of over the entire year as it assumes a constant evaporative demand and precipitation regimes and does not account for snow processes.Expanding all these simplifications and acknowledging imperfect knowledge and parameterization would require further analyses to better understand how they might affect the results shown here.Nevertheless, the uncertainties caused by these simplifications in the model by Porporato et al. (2004) might be partly overcome during the empirical connection made here between Porporato's model and Choudhury's formulation of the Budyko curve, as evidenced by the overall good performance of the developed BCP model in capturing the observed Q (Fig. 4).The third limitation of the current study lies in the steady-state assumption of the modeling framework.More specifically, the steady-state assumption is made in (i) catchment water balance and (ii) vegetation functioning.For point (i), a 5-year period does not necessarily guarantee zero-storage change.Nevertheless, the imbalance in water balance calculation under a steady-state assumption at a 5-year scale is generally very small (i.e., typically less than 6 % of P in arid regions and less than 3 % of P in humid regions) (Han et al., 2020).For point (ii), both Guswa's model for Z r and Donohue's model for L (see Sect. 2.1.5)adopted herein were developed for steady-state vegetation (i.e., mature and undisturbed vegetation).Applying these two models to immature (e.g., seedlings) and/or disturbed vegetation can be problematic, because immature and/or disturbed vegetation may have very different water use and carbon allocation strategies compared to steady-state vegetation (Donohue et al., 2017;Kuczera, 1987).However, the issues of vegetation age and disturbances are extremely complex and are well beyond our scope.Moreover, global datasets of vegetation age and disturbances are currently lacking.In this light, our modeled response of Q to eCO 2 should be regarded as if all vegetation was mature and undisturbed.Further efforts are needed to better quantify the age and disturbances of vegetation and to better understand the water use and carbon allocation strategies through the entire vegetation life cycle and under various types of disturbances.Supplement.The supplement related to this article is available online at: https://doi.org/10.5194/hess-25-3411-2021-supplement.
Author contributions.YY and TRM designed the study.YY performed the calculation and drafted the article.TRM, DY, YZ, SP, SP, and HEB contributed to discussion of results and article writing.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Flowchart of using the analytical models to detect the eCO 2 impact on Q.The terminologies used are explained in the following text (Sect.2.1.1 through to 2.1.5).

Figure 3 .
Figure 3. Locations of the catchments across the globe.The grey dots show the locations of the original 21 856 catchments, and red dots are the 2268 catchments that pass the selection criteria and are used herein.

Figure 4 .
Figure 4. Validation of estimated Q at catchments.(a) Model performance in predicting mean annual Q in 2268 catchments over 1982-2010.(b) Model performance in predicting Q trend in 2268 catchments during 1982-2010.(c) Model performance in predicting mean annual Q in 2268 catchments over 1982-2010 stratified by resource availability index category.(d) Model performance in predicting Q trend in 2268 catchments over 1982-2010 stratified by resource availability index category.The number of catchments in each resource availability index category are provided at the top of each subplot.The legend from (c) applies to (d).In (c) and (d), the upper and lower box edges represent the quartile division, the inner horizontal line is the median, the dots indicate the mean, and the dashed lines represent the 5th and 95th percentiles.

Figure 5 .
Figure 5. Modeled relative changes in L and Z r caused by eCO 2 .(a) Spatial distribution of relative change in L induced by eCO 2 during 1982-2010.(b) Same as (a) but for each resource availability index category.(c) Validation of predicted L change against in situ measurement during six Free-Air CO 2 Enrichment (FACE) experiments.Note that only FACE sites with undisturbed vegetation are used (see Donohue et al., 2017).(d) Spatial distribution of relative change in Z r induced by eCO 2 during 1982-2010.(e) Same as (d) but for each resource availability index category.In panels (b) and (e) the upper and lower box edges represent the quartile divisions, the inner horizontal line is the median, the dots indicate the mean, and the dashed lines represent the 5th and 95th percentiles.The number of grid cells in each resource availability index category is provided in Fig. 2.

Figure 6 .
Figure 6.Relative Q change induced by eCO 2 during 1982-2010 across the global vegetated lands.(a) Spatial distribution of relative change in Q induced by eCO 2 .(b) Boxplot of relative change in Q induced by eCO 2 for each resource availability index category.(c) Spatial distribution of relative change in Q induced by eCO 2 when only the physiological effect is considered.(d) Boxplot of relative change in Q induced by eCO 2 when only the physiological effect is considered for each resource availability index category.(e) Spatial distribution of relative change in Q induced by eCO 2 when only the aboveground and belowground structural effects are considered.(f) Boxplot of relative change in Q induced by eCO 2 when only the aboveground and belowground structural effects are considered for each resource availability index category.In (b), (d), and (f) the upper and lower box edges represent the quartile divisions, the inner horizontal line is the median, the dots indicate the mean, and the dashed lines represent the 5th and 95th percentiles.The number of grid cells in each resource availability index category is provided in Fig. 2.

Figure 7 .
Figure 7. Attribution of changes in Q between 1982-1985 and 2006-2010 across global vegetated lands.(a) Spatial distribution of changes in Q. (b-e) Spatial distributions of changes in Q induced by (b) changes in P , (c) changes in E P _M , (d) eCO 2 , and (e) changes in other factors (mainly rainfall intensity and climate-change-induced vegetation change).(f) Attribution of changes in Q between 1982-1985 and 2006-2010 averaged over the entire global vegetated lands.Values in parentheses represent 1 standard deviation of each response among all vegetated grid cells.

Figure 8 .
Figure 8. Attribution of changes in Q between 1982-1985 and 2006-2010 at grid boxes within each resource availability index (β) category.Values in parentheses represent 1 standard deviation of each response among grid cells within each resource availability index category.The number of grid cells in each resource availability index category is provided in Fig. 2.

Figure 9 .
Figure 9. Sensitivity of Q to eCO 2 and its relative importance to P and E P _M across the globe.(a) Spatial distribution of Q sensitivity to eCO 2 (percent change in Q per 1 % change in C a ).(b) Boxplot of Q sensitivity to eCO 2 for each resource availability index category.(c) Relative importance of eCO 2 to Q compared to changes in P on Q (percent change in Q per 1 % change in C a compared to percent change in Q per 1 % change in P ).(d) Boxplot of the relative importance of eCO 2 to Q compared to changes in P on Q for each resource availability category.(e) Relative importance of eCO 2 toQ compared to changes in E P _M on Q (percent change in Q per 1 % change in C a compared to percent change in Q per 1 % change in E P ).(f) Boxplot of the relative importance of eCO 2 to Q compared to changes in E P _M on Q for each resource availability category.In (b), (d), and (f), the upper and lower box edges represent the quartile divisions, the inner horizontal line is the median, the dots indicate the mean value, and the dashed lines represent the 5th and 95th percentiles.The number of grid cells in each resource availability index category is provided in Fig. 2.
P _M is the meteorological component of E P (without considering the increases in C a ).O represents factors other than eCO 2 that affect Z r , which effectively encode the climate-change-induced vegetation change.Functions f and g describe these relationships.Changes in E P and Z r are given by The first term on the right-hand side of Eq. (23) represents dQ caused by P change, and the second term represents dQ caused by eCO 2 .The third term calculates dQ induced by changes in E P _M and is calculated as ∂Q ∂E P dE P − ∂Q ∂E P ∂E P