Increasing aridity will not offset CO2 fertilization in fast‐growing eucalypts with access to deep soil water

Rising atmospheric [CO2] (Ca) generally enhances tree growth if nutrients are not limiting. However, reduced water availability and elevated evaporative demand may offset such fertilization. Trees with access to deep soil water may be able to mitigate such stresses and respond more positively to Ca. Here, we sought to evaluate how increased vapor pressure deficit and reduced precipitation are likely to modify the impact of elevated Ca (eCa) on tree productivity in an Australian Eucalyptus saligna Sm. plantation with access to deep soil water. We parameterized a forest growth simulation model (GOTILWA+) using data from two field experiments on E. saligna: a 2‐year whole‐tree chamber experiment with factorial Ca (ambient =380, elevated =620 μmol mol−1) and watering treatments, and a 10‐year stand‐scale irrigation experiment. Model evaluation showed that GOTILWA+ can capture the responses of canopy C uptake to (1) rising vapor pressure deficit (D) under both Ca treatments; (2) alterations in tree water uptake from shallow and deep soil layers during soil dry‐down; and (3) the impact of irrigation on tree growth. Simulations suggest that increasing Ca up to 700 μmol mol−1 alone would result in a 33% increase in annual gross primary production (GPP) and a 62% increase in biomass over 10 years. However, a combined 48% increase in D and a 20% reduction in precipitation would halve these values. Our simulations identify high D conditions as a key limiting factor for GPP. They also suggest that rising Ca will compensate for increasing aridity limitations in E. saligna trees with access to deep soil water under non‐nutrient limiting conditions, thereby reducing the negative impacts of global warming upon this eucalypt species. Simulation models not accounting for water sources available to deep‐rooting trees are likely to overestimate aridity impacts on forest productivity and C stocks.

vapor pressure deficit and reduced precipitation are likely to modify the impact of elevated C a (eC a ) on tree productivity in an Australian Eucalyptus saligna Sm. plantation with access to deep soil water. We parameterized a forest growth simulation model (GOTILWA+) using data from two field experiments on E. saligna: a 2-year whole-tree chamber experiment with factorial C a (ambient =380, elevated =620 μmol mol −1 ) and watering treatments, and a 10-year stand-scale irrigation experiment. Model evaluation showed that GOTILWA+ can capture the responses of canopy C uptake to (1) rising vapor pressure deficit (D) under both C a treatments; (2) alterations in tree water uptake from shallow and deep soil layers during soil dry-down; and (3) the impact of irrigation on tree growth. Simulations suggest that increasing C a up to 700 μmol mol −1 alone would result in a 33% increase in annual gross primary production (GPP) and a 62% increase in biomass over 10 years. However, a combined 48% increase in D and a 20% reduction in precipitation would halve these values. Our simulations identify high D conditions as a key limiting factor for GPP. They also suggest that rising C a will compensate for increasing aridity limitations in E. saligna trees with access to deep soil water under non-nutrient limiting conditions, thereby reducing the negative impacts of global warming upon this eucalypt species. Simulation models not accounting for water sources available to deep-rooting trees are likely to overestimate aridity impacts on forest productivity and C stocks.

K E Y W O R D S
deep soil water uptake, Eucalyptus saligna, forest productivity, vapor pressure deficit, water availability

| INTRODUC TI ON
The atmospheric CO 2 concentration (C a ) has increased by about 45% since the Industrial Revolution (Ciais et al., 2014), strongly altering the climate worldwide (IPCC, 2013;Lüthi et al., 2008). Globally, forests buffer the rise in C a by removing annually ~3.6 PgC from the atmosphere, which corresponds to roughly 34% of anthropogenic greenhouse gas emissions (Canadell et al., 2007;Pan et al., 2011; see Keenan & Williams, 2018 for a review). Recent observations attributed the increase in forest carbon uptake of the last decades to rising C a (Keenan et al., , 2020Liu et al., 2019). However, forest C uptake is strongly affected by the water cycle (Denmead & Shaw, 1962;Humphrey et al., 2018;Nemani et al., 2003). Two of its components are primarily limiting: low water available in the soil (e.g., Denmead & Shaw, 1962;Grassi & Magnani, 2005; Van der Molen et al., 2011) and high atmospheric evaporative demand (Novick et al., 2016;Sanginés de Cárcer et al., 2018; see Grossiord et al., 2020 for a review). From now on, we use the term aridity as the hydric stress exerted on the vegetation by low soil water availability (land aridity) and/or by elevated evaporative demand (atmospheric aridity), following Berg et al. (2017).
Over the next decades C a is projected to keep rising, modifying Earth's climate (IPCC, 2013;Webb & Hennessy, 2015). In the absence of nutrient limitation (see Ellsworth et al., 2017), rising C a has led to an increase in forest productivity and carbon sequestration during the last few decades (e.g., Cheng et al., 2017;Haverd et al., 2020;Keenan et al., 2013;see Walker et al., 2019). More precisely, eC a enhances intrinsic water-use efficiency (WUE; De Eamus & Jarvis, 1989;Keenan et al., 2013) due to decreased stomatal conductance (e.g., Ainsworth & Long, 2005;Farquhar & Sharkey, 1982;Medlyn et al., 2001) and increased net assimilation (e.g., Barton et al., 2012;Kelly et al., 2016). Higher WUE may sustain non-limiting soil moisture for a longer period during the growing season, resulting in a further stimulation of gross primary production (GPP). In parallel, the increase in CO 2 and other greenhouse gas concentrations in the atmosphere has led to an overall augment in temperature, 2020 which is expected to keep rising in the next decades.
Following temperatures, exponential increases in vapor pressure deficit (D) are expected, thus augmenting atmospheric aridity (Berg & Sheffield, 2018;Lin et al., 2018;Scheff & Frierson, 2014). Similarly, changes in annual precipitation amounts and patterns may limit soil water availability particularly during the growing season (Doughty et al., 2015;Gazol et al., 2017;Trenberth et al., 2014;Williams et al., 2010). At declining soil water potential, trees typically close their stomata to avoid loses in hydraulic conductance (Peñuelas et al., 2017;Reich et al., 2018;Sperry & Tyree, 1990). This limits C substrate for photosynthesis, hence reducing assimilation. Additionally, trees may reduce leaf mesophyll conductance and decrease photosynthetic potential under drought stress, which further limits assimilation Grassi & Magnani, 2005;Keenan et al., 2011). Trees that are able to mobilize water from deep soil reservoirs have shown to be less affected by episodic drought events and to maintain productivity at a higher rate (e.g., McDowell et al., 2019; Ripullone et al., 2020). This ability is thought to be key for the persistence of water-demanding eucalypts in the harsh Australian climate (Christina et al., 2017;Laclau et al., 2013;Zeppel et al., 2008).
Recent findings suggest that when trees are able to access shallow water tables, even extreme sporadic drought has little impact on tree growth (Sousa et al., 2020). Nevertheless, in the long term, deep soil water availability is tightly linked to precipitation patterns over land (Berg et al., 2017). In addition, increased water consumption by vegetation due to increasing evaporative demand may reduce such deep soil water reservoirs (Mu et al., 2021). In the context of increasing aridity, this may limit productivity and survival of eucalypts relying on such deep soil water reservoirs.
Whether rising C a will continue to enhance forest productivity during the next decades, or it will be counterbalanced by aridity increases, is still under debate. On the one hand, tree mortality may increase in a hotter and drier world (Allen et al., 2015). This will limit forest C uptake and reduce forest C stock. Recent observations suggest that drought and heat stress is accelerating tree mortality in mature Amazon forests, thus counterbalancing eC a -driven biomass gains, and reducing carbon sink potential (Brienen et al., 2015). A similar trend is expected for the African tropical forests in the near future (Hubau et al., 2020). Similarly, D increases have been suggested to have the potential to offset eC a fertilization at a global level (Yuan et al. 2019). Furthermore, recent findings suggest that trees do not benefit from eC a conditions when exposed to a severe hot-drought stress (e.g., Birami et al., 2020;Duan et al., 2014;Gattmann et al., 2020). On the other hand, data from largescale decadal manipulative fertilizing experiments involving eC a (FACE experiments) of young forests growing under non-limiting nutrient availability conditions show a persistent stand biomass increase (Walker et al., 2019). A number of process-based models suggest that such eC a stimulation on forest net productivity will be sustained during the next decades, for environments ranging from semiarid Mediterranean to tropical forests (e.g., De Kauwe et al., 2016;Fatichi et al., 2016;Keenan et al., 2011;Sperry et al., 2019).
In all cases, stimulation of photosynthesis and water savings due to reduced g s at eC a were found to exceed the detrimental effects of increasing aridity (e.g., Keenan et al., 2020;Swann et al., 2016). Also, the aforementioned ability to access water-enriched deep soil layers may contribute to ameliorate stand responses, thus allowing trees to take more advantage of eC a conditions (e.g., Humphrey et al., 2018;Laclau et al., 2013;Nadal-Sala et al., 2019).
In eucalypts, early experiments with potted Eucalyptus grandis seedlings demonstrated that eC a led to an increase in net assimilation and tree biomass, if N and P were not limiting (Conroy et al., 1992). Subsequently, similar conclusions were drawn from numerous controlled experiments on different eucalypt species grown in nonnutrient limiting conditions (e.g., Duan et al., 2019;Quentin et al., 2015). However, in P-limited soil conditions, mature E. tereticornis trees grown in eC a at EucFACE did not exhibit higher aboveground growth, despite higher leaf-level photosynthesis (Ellsworth et al., 2017;Jiang et al., 2020). In growth chamber studies with saplings, eC a limited the negative effects of moderate drought stress by increasing leaf-level water-use efficiency (WUE L ) and stimulating fine root development, which delayed drought stress (e.g., Apgaua et al., 2019;Atwell et al., 2007Atwell et al., , 2009Roden & Ball, 1996). However, eucalypt seedlings exposed to a progressive, intense drought did not exhibit growth enhancement or reduced hydraulic vulnerability due to eC a (Duan et al., 2014(Duan et al., , 2018, again suggesting that the benefit of eC a was dependent on the duration and intensity of the drought stress. Therefore, whether rising C a will continue to stimulate forest productivity or will be offset by concomitant aridity increases is yet not solved, with severe implications for forest C cycling. Here, we used data from two field experiments in a young E. saligna plantation in Australia, where trees had access to deep soil water, to calibrate the GOTILWA+ process-based forest simulation model. Then, we quantified the sensitivity of E. saligna to several environmental drivers related to climate change. We particularly focused on the main and interactive effects of projected rise in D, reduced precipitation, and increasing C a , upon tree productivity and transpiration. More specifically our research questions were: (1) To what extent does increasing D and decreasing precipitation limit productivity of E. saligna? (2) Will access to deep soil water ameliorate soil and atmospheric aridity impacts on productivity? and (3) How will stand WUE respond to combined aridity and C a increases? 2 | MATERIAL AND ME THODS  Barton et al., 2010;Heroult et al., 2013). Upper soil layers are an alluvial formation of low-fertility sandy loam soils, with low organic matter content. Below 100 cm depth, there is a waterenriched clay layer that limits drainage from the upper parts of the soil and acts as a water reservoir . In the HFE, E. saligna has been observed to develop deep roots, with capacity to access water to ~4.5 m depth .

| Site description and experimental setup
In the HFE, two separate experiments were co-located and performed using E. saligna: a whole-tree chamber experiment (hereafter referred as WTC) and an irrigation experiment (hereafter referred as "Ir"). For the Ir experiment, E. saligna seedlings were planted in April 2007 at a density of 1000 trees ha −1 . Four square plots per treatment (i.e., Control and Irrigation, size = 38.5 × 41.6 m, ~160 trees per plot) were randomly determined. The plots were laid out and planting rows were ripped to a depth of 30 cm; deeper ripping and mounding was unnecessary in the sandy top soil. A 2-m wide strip was sprayed with herbicide to control weed growth prior to planting.
In each plot, a computerized irrigation system capable of supplying a uniform distribution of water across each plot at a rate of 7 mm h −1 was installed. The irrigation treatment consisted of the application of 15 mm every fourth day (delivered at night). Irrigation was withheld during periods of high rainfall to avoid leaching of nutrients or water-logging. For each plot, we performed 14 measurement campaigns between 2007 and 2018, in which we measured the diameter at breast height (DBH, in cm) at ~1.3 m for each individual tree and then we averaged the DBH per plot. under a 2 × 2 factorial combination of atmospheric CO 2 and watering treatments. The CO 2 treatment levels were ambient CO 2 (aC a , 380 μmol mol −1 ) and elevated CO 2 (eC a , 620 μmol mol −1 ); see Barton et al. (2010Barton et al. ( , 2012 for further details. The watering treatment levels were designated wet (irrigated throughout the monitoring period) and dry-down (water withheld from October 2008 to February 2009+ rain exclusion); see Duursma et al. (2011) for further details.
There were three chambers for each C a × watering treatment. The chambers were surrounded by control trees, planted at a density of 1000 trees ha −1 .

| Leaf gas exchange
Leaf-level gas exchange was measured on WTC trees during nine campaigns, from September 2007 to March 2009. Measurements were made on sun-exposed, fully expanded leaves with three replicates per WTC per measurement campaign, using two portable infrared gas analyzer systems (LI-6400, Li-Cor, Inc.) with a 6 cm 2 chamber. Measurements were taken at saturating photosynthetic photon flux density (PFD > 1800 μmol m −2 s −1 ). Assimilation in relation to [CO 2 ] in the intercellular space (A n /C i ) curves were generated by increasing C a inside the 6 cm 2 chamber in nine to 10 steps.
The first measurement was taken at ambient C a . We initially took a measurement at sub-ambient C a (ca. 37 μmol mol −1 , near CO 2 compensation point). We then increased C a progressively, to a maximum of 2000 μmol mol −1 . Leaf temperature (T leaf , in °C) and other standard variables were also measured. These measurements were taken under prevailing temperature and humidity conditions within each WTC at the time of measurement.

| Whole-tree aboveground carbon and water fluxes
For each WTC, aboveground whole-tree net carbon assimilation (A n,tot,inst = carbon assimilation minus aboveground growth and maintenance respiration, in g C m −2 leaf area), and transpiration (E inst ) per unit of leaf area (kg H 2 O m −2 leaf area), were recorded at 14min intervals from April 2008 to March 2009 (Barton et al., 2010).
We calculated average daylight D inside each WTC on a daily basis.
Daily standing leaf area (L, in m 2 ) per tree calculation is described elsewhere (Barton et al., 2012). We integrated A n,tot,inst and E inst daily, from April 2008 to March 2009, in order to obtain A n,tot (in g C m −2 leaf area day −1 ) and E (in kg H 2 O m −2 leaf area day −1 ). Periods of equipment failure, failed seals to the chamber, and people entering the chamber for measurements were excluded from the final dataset.

| Soil water content
Volumetric soil water content (SWC, in %) for each WTC was continuously measured in the top 10 cm of soil with a Theta probe (Delta T instruments) throughout the study period. Additionally, for each WTC, SWC was measured every 2 weeks to a depth of for each soil depth, we selected 0.5 m as the reference soil depth because the agreement between daily and discrete SWC measurements worsened with depth (see Figure S1). Once daily SWC 50 was obtained for each WTC, we calculated the relative SWC of the upper 0.5 m soil layer (S up ) by dividing daily SWC 50 by an average of the highest 5% daily SWC 50 values (i.e., S up = SWC 50 / SWC 50,max ). Below 0.5 m, we assumed that trees were able to access the deep soil water reservoir, with a relative SWC (S deep ) that was set constant to 0.7 based on the relative SWC below 0.5 m measured before the dry-down experiment (see Figures S2 and S3b).

| Meteorological data
Daily meteorological data for the 2008-2017 period were obtained from the Richmond RAAF base meteorological station (Station 067105, bom.gov.au; see Figure S4), a meteorological station located about 3 km east of HFE. We accessed data on daily precipitation (P, in mm day −1 ), maximum and minimum temperature (T max and T min , in °C), photosynthetic active solar irradiation (Q, in MJ m −2 day −1 ), and air relative humidity (RH, in %). One and 2-day gaps in data were linearly interpolated. Wider gaps in T max and T min data were interpolated from Windsor Bridge meteorological station (REF =067095, bom.gov.au), located about 5 km westward from the HFE, by linear regression (T max : n = 2821, t = 455.1, From RH and T min we calculated daily mean air water vapor pressure (VP, in kPa) according to the Clausius-Clapeyron relationship (Allen et al., 1998). Potential evapotranspiration (PET, in mm day −1 ) was calculated daily according to the Penman-Monteith equation (Allen et al., 1998).

| Model
GOTILWA+ is a process-based dynamic model that focuses on forest carbon and water balances (www.creaf.uab.es/gotilwa; Gracia et al., 1999;Keenan et al., 2009;Nadal-Sala et al., 2017). GOTILWA+ has been extensively tested and validated against both flux data and stand data in many previous studies (e.g., Bugmann et al., 2019;Gracia et al., 2007;Keenan et al., 2009Keenan et al., , 2010Kramer et al., 2002;Morales et al., 2005). GOTILWA+ has been applied to project climate change and forest management impacts on forest carbon and water fluxes (e.g., Reyer et al., 2017;Sabaté et al., 2002). Also, Nadal-Sala et al. (2019) found that GOTILWA+ was able to properly capture phreatic water uptake in a Mediterranean riparian forest dominated by Robinia pseudoacacia L. Briefly, GOTILWA+ simulates carbon and water fluxes in a monospecific forest stand accounting for environmental drivers (e.g., [CO 2 ], temperature, soil water, D, etc.), forest structure, and management, and assuming two leaf layers (e.g., sun and shade leaves) in a "big leaf" approach. For a detailed description of the model's process hierarchy and integration, see the references above. In the following subsections we summarize the main processes involved in our analysis.

| Leaf gas exchange in the absence of drought stress
In GOTILWA+, stomatal conductance (g s , in mol m −2 s −1 ) and leaflevel net assimilation (A n , in μmol m −2 s −1 ) are obtained iteratively at hourly time steps, while closing the leaf energy balance according to Gates (1962Gates ( , 2012. In contrast to previous versions of GOTILWA+, here g s was derived from the optimal stomatal behavior model (Medlyn et al., 2011, Equation 1): where g 0 is the residual conductance (in mol m −2 s −1 ), A n is the net assimilation, g 1 is the stomatal operating point parameter (in kPa 0.5 ), D is the vapor pressure deficit (in kPa), and C a is the atmospheric CO 2 concentration. A n is calculated according to the Farquhar-von Caemmerer-Berry model (FvCB;De Pury & Farquhar, 1997;Farquhar et al., 1980). V cmax and J max in the FvCB model are leaf temperature dependent (e.g., Farquhar et al., 1980;Harley et al., 1992). Dark respiration (R d ) is assumed to increase exponentially with leaf temperature. K c , K o , and Γ* (Rubisco oxidase and carboxylase kinetic constants, and CO 2 compensation point, respectively) are also temperature dependent.
(1) where S tot is the relative soil water content of the whole soil profile, S tot,max is the value of S tot above which there is no water limitation upon photosynthesis, and S tot,min is the critical S tot value below which β = 0. The parameter q describes the nonlinearity of the responses to decreasing relative SWC. For a given S tot , the water-limited values of where S up is the relative SWC, and S deep (fraction) is the relative SWC in the deep soil layers. Here, as noted above S deep is assumed to be constant, with a stand-specific value of 0.7 in the HFE. Finally, α is the proportion of water uptake from S deep when S up is at field capacity.
Close-to-zero α values imply that the tree accesses water from the upper soil layer, while higher α values imply increasingly greater water mobilization from the deep soil (see Figure S4). Other impacts from reductions in S tot considered in GOTILWA+ include increases in leaf temperature due to loss of leaf temperature homeostasis at declining g s (Gates, 1962;2012), which increases leaf maintenance respiration and accelerates leaf turnover by reducing mean leaf life span. Similarly, fine root carbon turnover velocity also increases at decreasing S tot , which accounts for observed higher fine root turnover rates (e.g., Bloom et al., 1985;Gaul et al., 2008;Meier & Leuschner, 2008)

| Model parameterization
For the stomatal model (Equation 1), g 0 and g 1 parameters were obtained from Herault et al. (2013). Photosynthesis activation and deactivation energies for E. saligna were obtained from Lin et al. (2013).
K c , K o , and Γ* and their activation and inactivation energies were obtained from Bernacchi et al. (2001). See Table S1 for the source of all relevant parameters.  2008). We used the "BayesianTools" package  to run the MCMCs. The first 50 k simulations for each chain were removed from the posterior as a burn-in. Convergence was assessed by comparing the three chains for each calibration, and under a conservative Gelman-Rubin's test score (Gelman & Rubin, 1992) of 1.1.

| Model validation
We validated the FvCB model against an independent A n /C i dataset of observations under no water limitation (S up ≥0.95), measured at 350 ≤ C a ≤700 μmol mol −1 (n = 77; Step 1 in Table S2). We validated the FvCB+β model with another independent dataset of A n /C i curves,  (4) reduced precipitation and reduced water uptake (R α R P , ΔP = −20% and α = 0.225). We compared simulated A n,tot corrected at D = 1 (A n,tot,corr , g C m −2 leaf day −1 kPa −1 ) with the measurements from the dry-down WTC experiment when S up was below 0.3. As a Shapiro-Wilk test confirmed non-normality of the data, significant differences in A n,tot,corr between scenarios were assessed through the Mann-Whitney-Wilcox U nonparametric test.
Modeled gas exchange fluxes at the stand scale were validated against WTC gas exchange data. We ran two simulations for the orig- Step 5). Modeled daily stand A n,tot and E were correlated with average daily daylight D, and we compared simulation results with observations from single-tree WTC under both C a treatments.
We validated GOTILWA+ growth projections against data from the Ir experiment. Two separate 10-year simulations, control without irrigation and with irrigation (3.75 mm day −1 ), were run for the 2008-2017 period in order to evaluate the ability of GOTILWA+ to capture stand-scale responses to irrigation. We compared simulated mean stand DBH against observed mean stand DBH under both irrigation treatments (Step 6). See Figure 1 for a scheme of the approach.

| Climate sensitivity scenarios
After model validation, the sensitivity of the E. saligna plantation to changes in environmental conditions was evaluated. Three environmental drivers were considered in this analysis, with upper boundaries of change based on IPCC (2013) and Webb and Hennessy (2015).
Five levels of C a were simulated (C a =350, 400, 500, 600, and 700 μmol mol −1 ). For each C a level, photosynthesis downregulation was accounted by linearly extrapolating observed changes in Notes 1). Additionally, six levels of precipitation were considered: +5%, +0%, −5%, −10%, −15%, and −20% of annual precipitation; no seasonal variations were considered. Although long-term deep soil water dynamics are strongly linked to precipitation patterns (Berg et al., 2017;Berg & Sheffield, 2018 and references therein; Seneviratne et al., 2010), there are still important uncertainties about the dynamics of such soil water reservoirs in SW Australia, mostly driven by the impact of vegetation (Mu et al., 2021). Therefore, for each precipitation reduction scenario, we assumed a linear decrease in S deep proportional to the reduction in precipitation.
To evaluate the effect of increased D, we assumed no changes in atmospheric water vapor pressure (Trenberth et al., 2007). We assumed changes in air temperature of −1, 0, 1, 2, 3, 4°C and calculated D for each air temperature scenario following the Clausius-Clapeyron relationship. Rising D due to increasing temperature is nonlinear; therefore, we report those changes as ΔT(D), in °C.
The mean annual D for the temperature change scenarios (−1°C to 4°C) varied between 1.33 and 2.21 kPa, respectively, which corresponded to a range of −11% to 48% increases in annual average D, respectively.
The initial stand density in the E. saligna plantation was 1000 trees ha −1 with initial DBH =2.6 cm. A total of 180 simulation TA B L E 1 Prior and posterior distributions for the three Bayesian calibrations performed. For all parameters, a uniform prior distribution was assumed. For each parameter, median values and 5%-95% CI of their posterior marginal distribution are noted. aV cmax,25,ref and aJ max,25,ref are V cmax and J max at 25°C under well-watered conditions for the aC a treatment (μmol m −2 s −1 ). eV cmax,25,ref and eJ max,25,ref are V cmax and J max at 25°C under well-watered conditions for the eC a treatment. S min is the minimum S tot below which A n = 0, and S max is the minimum S tot above which there is no drought limitation on photosynthesis assimilation (A n , both in %). Q is a unit-less scale coefficient according to Equation (2). α is the fraction of water uptake below 50 cm soil depth when S up is at saturation, according to Equation (3 Table S2).

| Data analysis
Data management and statistical analysis were performed using R software, version 3.6.1 (R Development Core Team, 2019). Bayesian model inversions were performed with the "BayesianTools" package . GAM analyses were performed with the "mgcv" package (Wood, 2017). For all least-squares regressions, homoscedasticity and normality of the residuals were visually checked, and data were log-transformed when necessary.

| Model parameterization and validation of leaf gas exchange
Inverse calibration of the Farquhar model captured measured observations. After that point A n,tot declined, and even turned negative at very high D (>3 kPa at aC a and >4 kPa at eC a ), when stomata were likely closed and respiration costs were higher than photosynthetic gains.
Simulations using aC a compared well with observed transpiration responses to D, although the model slightly underestimated F I G U R E 1 Scheme of model calibration validation setup. Experiment 1 is the 10-year long irrigation experiment in the Hawkesbury Forest Experiment plantation of Eucalyptus saligna (density =1000 trees ha −1 ), and Experiment 2 is the whole-tree chamber (WTC) experiment of E. saligna (see details in Barton et al., 2010Barton et al., , 2012Duursma et al., 2011). The two main references for GOTILWA+ model are provided. Colors in the boxes indicate where calibration (purple color) and validation (tan color) procedures occurred, and the origin for such data or the reference for the model calibrated in such step. aC a , ambient CO 2 atmospheric concentration (380 μmol mol −1 ); eC a , elevated CO 2 atmospheric concentration (620 μmol mol −1 ); A n , net assimilation per unit of leaf area; C i , CO 2 concentration in the sub-stomatal chamber; DBH, diameter at breast height; FcVB+β-model, Farquhar, von Caemmerer and Berry photosynthesis model (Farquhar et al., 1980)  Modeled ratio of water uptake from below 50 cm soil depth (WU deep ) with respect to the total water uptake (WU total ), as a function of upper soil water content (S up ). It was calculated with the median α coefficient ± 95% CI (black line ± gray area) from Table 1. White circles represent measured mean ± 95% CI (n = 6) fraction of WU deep WU total −1 below 50 cm during two discrete sampling campaigns,

| Model sensitivity to changes in deep soil water availability
Simulated daily A n,tot normalized to D = 1 kPa (A n,tot,corr ) was not different than WTC measurements when S up was below 0.3 (Figure 5a). However, when deep soil water uptake ( year with average precipitation) in Figure 5b. Conversely, water savings due to decreased g s at eC a reduced the number of days when S up was limiting A n ( Figure S7). Finally, average daily deep soil water uptake from S deep decreased by 44% in the R α scenario (average 0.6 ± 0.14 mm day −1 ; Figure 5c) compared to the ambient conditions (1.1 ± 0.23 mm day −1 ).

| Sensitivity of modeled canopy gas exchange to environmental drivers
The sensitivity analysis suggests that an increase in C a of 75% from 400 to 700 μmol mol −1 (leaving all other drivers unchanged) increases WUE L by 55%, driven by a 20% increase in A n and a 23% reduction in E (Figure 6). Scaling these values to the stand level by accounting for changes in leaf area, the overall change in WUE was very similar (+56%), but was composed of a larger increase in C uptake (33%) and a smaller reduction in E (−15%; Table 2; Figure S8). and GOTILWA+ outputs (red circles), and at 380 μmol mol −1 C a (aC a ) and at 620 μmol mol −1 C a (eC a ). In each panel, the colored areas represent the 95% CI interval for the response of that variable to D increases after a generalized additive model analysis with a fixed maximum of 5 df, with observed interval in blue and simulated interval in red. Note than model simulations cover a longer period than observations, and therefore the D distribution vary water uptake, and a higher reliance of E. saligna on deep water access under more arid conditions. However, the increase in WUE L to rising C a was not affected by such reduction in precipitation.
Under constant 400 μmol mol −1 C a conditions, combined increases in D and reductions in P decreased both, A n and WUE (−10% and −25%, respectively). Considering all three drivers together, above C a =500 μmol mol −1 GPP fertilization was not fully offset by aridity increases, even under the most extreme scenario (ΔD = +48% and ΔP = −20%; Figure 6).

Canopy transpiration was strongly affected by increases in D,
while decreases in precipitation only slightly affected E, due to the ability of this stand to obtain water from deep soil layers. Although F I G U R E 5 Sensitivity of GOTILWA+ model outputs to changes in deep soil water access/uptake (α) for the Eucalyptus saligna plantation in Richmond, New South Wales, Australia. (a) Comparison between GOTILWA+ projected whole-tree canopy net carbon uptake corrected at D = 1 kPa (A n,tot,corr , g C m −2 leaf day −1 kPa −1 ) under four scenarios: ambient conditions (Ambient), a 20% reduction in precipitation (R P ), a 50% reduction in α (R α ), a 20% reduction in P, and a 50% reduction in α combined (R P,α ). Both simulations and observations included correspond to the days in which observed relative soil water content in the upper 50 cm of the soil (S up ) in whole-tree-chambers under dry-down experiment was below 0.3. Common letters indicate nonsignificant differences after pairwise Mann-Whitney tests (p < 0.05). (b) Comparison of 7-day moving averages of daily whole-tree canopy GPP (g C m −2 leaf day −1 ) during the 2013 year (a year with average precipitation), for the Ambient and R α simulations. Shadowed area corresponds when S up < 0.3, in R α simulation (orange) and both in R α and Ambient simulations (red). (c) Relationship between monthly averaged daily water uptake from below 50 cm soil depth (PWU, in mm day −1 ) in relation to S up , for Ambient (green, circles) and R α (red, triangles) scenarios. Mean monthly values +95% CI are represented, as well as the linear regression for each simulation (*p < 0.05; **p < 0.01) eC a alone lowered g s by 33%, reductions in E were estimated to be about 120 mm year −1 (−15%; Table 2). Furthermore, increments in D offset eC a induced water savings, resulting in similar transpiration rates as under aC a . When considering all changes together, the annual transpiration was reduced by a mere 23 mm year −1 (3% of the reference values).

| Changes in aboveground biomass related to environmental drivers
Modeled changes in aboveground biomass stock (ABS) were strongly affected by changes in C a , with a 62% increase in ABS at 700 μmol mol −1 relative to aC a . However, increases in atmospheric evaporative demand and decreases in precipitation reduced this stimulation. For instance, in eC a , a 20% reduction in precipitation reduced ABS to a 50% increase relative to aC a . Similarly, a 48% increase in D (equivalent to an anomaly in temperature of +4°C) reduced C a fertilization to a 48% increase relative to aC a (Figure 7b; Table 2). In aC a , +4°C ΔT(D) resulted in a 14% reduction in ABS, whereas a 20% reduction in P led to a 13% reduction. The interactive effects of +4°C ΔT(D) and −20% ΔP resulted in a reduction of ABS by 27% compared to the control conditions (Figure 7c). When all three environmental drivers were assessed in combination, ABS increased by 35%, which was approximately half of the stimulation when accounting for eC a alone.

F I G U R E 6
Responses of simulated canopy leaf gas exchange-that is, net assimilation (A n ), transpiration per unit of leaf area (E), and water-use efficiency at the leaf level (WUE L, calculated as A n /E)-in a Eucalyptus saligna plantation to pairwise changes in key environmental drivers with respect to reference conditions (i.e., C a =400 μmol mol −1 , ΔT(D) =0°C and −ΔP = 0%, black dots, percent = 100%). Percentages depict extreme responses in white (increase) and red (decrease) letters

| DISCUSS ION
We investigated E. saligna responses to key environmental drivers that are expected to change in the near future. We parameterized and validated the GOTILWA+ model for an E. saligna plantation with field observations from two experiments at a site close to its natural distribution. We then used GOTILWA+ to simulate E. saligna sensitivities of key ecophysiological processes to changes in vapor pressure deficit (D), precipitation (P), and atmospheric carbon dioxide concentration (C a ), and assuming no nutrient limitation. Our results suggest that young E. saligna plantations will increase their C uptake and growth at rising C a under this assumption. Thus, the CO 2fertilizing effect will be halved, but not offset by increasing aridity.
We found that the access to the deep water-enriched soil layers will play a key role in this response, as it will ameliorate drought stress in E. saligna. Nevertheless, this highlights the sensitivity of E. saligna to the future dynamics of deep soil water reservoirs, which are highly uncertain. Furthermore, we found that a 48% increase in D will limit the eC a fertilization effect on ABS at a similar magnitude as a 20% reduction in precipitation. This emphasizes the role that D increases will play for forest functioning under global warming, and the need to properly capture tree responses to D increments in process-based models.

| Photosynthesis responses to eC a
The derived FvCB model parameters, V cmax and J max at 25°C, were similar to the ones reported for E. saligna trees growing in a nearby common garden (Heroult et al., 2013 (Table 1), which is in accordance with previous observations (e.g., Long et al., 2004;Medlyn et al., 1999).
Photosynthetic downregulation may occur due to changes in leaf morphology and biochemistry, such as reductions in specific leaf area (SLA), and leaf N and P content (e.g., Ainsworth & Rogers, 2007;Ellsworth et al., 2004;Lewis et al., 2013;Wang et al., 2018).
However, we found no differences in N dry weight per unit of leaf area (g m −2 leaf) between aC a and eC a treatments (see Notes S2 and Table S3). This is similar to observations in the nearby EucFACE experiment (Crous et al., 2019). E. saligna trees growing under eC a presented significantly lower SLA, which has been reported to reduce photosynthetic potential (e.g., De la Riva et al., 2016;Flexas et al., 2008;Sperlich et al., 2015). This might partially explain the downregulation of photosynthesis under eC a conditions. Also, in the abovementioned study by Crous et al. (2019), leaf phosphorous concentration was additionally limiting photosynthesis in trees growing TA B L E 2 Simulation outputs for different ecophysiological variables under the three environmental drivers considered (C a , ΔT(D), and ΔP), after a simulation of Eucalyptus saligna plantation stand growth over 10 years. ABS is the aboveground biomass stock (in Mg ha −1 ), GPP is the canopy gross primary productivity (in Mg C ha −1 year −1 ), NCU t is the total forest net carbon uptake (NCU t = GPP − Maintenance respiration − Growth respiration, Mg C ha −1 year −1 ), CUE is the carbon-use efficiency (GPP/NCU, in g C g −1 C), g s is the average daytime stomatal conductance during the year (in mmol m −2 s −1 ), LAI is the leaf area index (in m 2 m −2 ), E is the stand transpiration (in mm year −1 ). WUE is the canopy water-use efficiency (GPP/E) (g C kg −1 H 2 O). All variables but ABS are reported as mean annual values for the 6-to 10-year period. ABS is reported as the final ABS after the 10 simulation years at eC a conditions. Although leaf P content was not analyzed here, we acknowledge that it might have also played a role in photosynthesis downregulation, providing the low P availability in the soil of the HFE.

| Mobilizing water from deep soil layers ameliorates drought stress
Declining soil water availability reduced A n , as expected. However, we found that A n declined less in response to soil dry-down than that reported previously for E. saligna. Thus, Lewis et al. (2013) reported 80% reductions in A n on potted E. saligna seedlings after only 10 days under imposed experimental dry-down.
Contrastingly, we found A n for E. saligna in the WTCs to decline during a 4-month-long complete rain exclusion by only about 40% ( Figure 2a). This suggests that drought responses were buffered in E. saligna growing in the HFE. There, E. saligna has been suggested to be highly reliant on deep soil water uptake . Our results confirm this observation, especially during autumn, when upper soil layers become dry. The ability of fastgrowing eucalypts to obtain water from deep soil layers has already been widely reported (e.g., Calder et al., 1997;Christina et al., 2017;Drake et al., 2018;Laclau et al., 2013;White et al., 2000;Zeppel et al., 2008). In the HFE, the ability to mobilize Thus far, stand-scale simulation models tend to better describe tree water uptake than DGVMs (e.g., Mu et al., 2021;Ukkola et al., 2016).
However, this is generally done at the expense of an exhaustive soil biophysical parameterization and root vertical distribution assumptions (Vereecken et al., 2016). It is challenging for DGVMs to go beyond the stand scale because uncertainties in soil composition, in plant rooting depth and in deep soil water availability increase at larger scales (e.g., Seneviratne et al., 2010;Sitch et al., 2015). Also, unfortunately rooting depth and vertical root distribution are not routinely measured in field experiments, and actual water source partition measurements, such as the ones we used here to validate model projections, are scarcely available. Our study shows that such measurements are necessary to properly account the impacts of drought on tree productivity, and that this could potentially improve the accuracy of simulation models.

| Sensitivity of GPP and stand biomass to increasing C a
In the absence of nutrient limitation, eC a is projected to increase GPP and enhance biomass stock in the HFE plantation. If E. saligna can access water from deep soil layers, our projections suggest that an increase in aridity will only partially offset the positive eC a effect on GPP. Our simulations showed that GPP increased by 16%-33% at 700 μmol mol −1 C a depending on the degree of aridity.
Sustainment of increasing forest C stock with rising atmospheric CO 2 will depend on diverse environmental and stand-specific conditions, such as succession stage and nutrient availability (e.g.,  (Jiang et al., 2020;Riikonen et al., 2004), changes in leaf turnover rate at eC a (e.g., Duursma et al., 2016), and physiological and structural acclimation to eC a (e.g., Sperry et al., 2019), are likely to further limit biomass production stimulation at eC a .

| Sensitivity of gas exchange and stand biomass to C a × aridity
As D rises, trees reduce stomatal conductance (Creek et al., 2020;Martínez-Vilalta et al., 2014). This limits CO 2 substrate for photosynthesis, and therefore limits primary productivity. Our results suggest that in the HFE stand, this may limit forest productivity to a similar degree as projected reductions in precipitation, which supports previous observations in subhumid regions (e.g., Grossiord et al., 2020;Novick et al., 2016;Sanginés de Cárcer et al., 2018). This further may indicate an overall influence of D on reported variations in GPP responses to eC a from different FACE experiments (e.g., Medlyn et al., 2015;Walker et al., 2019).
Furthermore, properly capturing D limitations on productivity is still challenging in simulation models (e.g., Grossiord et al., 2020;. Hence, responses to rising D are still a major source of uncertainty in forest C uptake projections under global warming. Even more, D itself is affected by transpiration, as vegetation is a major driver of atmospheric humidity over land (e.g., Berg & Sheffield, 2018;Berg et al., 2017;Lemordant et al., 2018). Therefore, changes in forest transpiration due to eC a and reduced water availability could potentially increase D, hence amplifying drought limitations over forest productivity.
As expected, plants reduced stomatal conductance as C a increased. This resulted in lower transpiration rates if precipitation and evaporative demands remained unchanged. At C a =700 μmol mol −1 g s declined (−32%). However, transpiration was reduced by a mere 15%.
The slight increase in leaf area (9%) under eC a does not account for such differences. Instead, lower g s led to non-limiting SWC for a longer period, thus reducing the number of days in which E. saligna experienced severe drought stress, without severely reducing overall annual stand transpiration. This phenomenon was already observed by Fatichi et al. (2016) in water-limited ecosystems. It is also in accordance with previous reports, which did not report changes in evapotranspiration fluxes in surrounding eucalypt woodlands growing at eC a (Gimeno et al., , 2018).
Elevated C a was found to be the dominant driver of changes in GPP and biomass, while the combined effects of reduced precipitation and increased D halved eC a fertilization. Intense drought stress has been reported to potentially offset the positive effects of eC a (e.g., Andresen et al., 2016;Barbeta & Penuelas, 2017;Birami et al., 2020;Peters et al., 2018). The positive responses projected in the HFE to combined eC a and increasing aridity are related to the ability of the eucalypts to access deep soil water reservoirs. Accounting both for productivity gains and water savings, a simulated 55% increase in WUE L at C a =700 μmol mol −1 is consistent with previous reports about WUE L increases with rising CO 2 (e.g., Battipaglia et al., 2013;Eamus, 1991;Keenan et al., 2013;Mastrotheodoros et al., 2017). Interestingly, rising D has been reported to enhance intrinsic WUE (e.g., Ficklin & Novick, 2017;Zhang et al., 2019), but to reduce both WUE L and WUE (e.g., Grossiord et al., 2020;Jarvis & McNaughton, 1986). Increases in WUE L in the E. saligna saplings growing under eC a conditions in the HFE have been thoroughly discussed in Barton et al. (2012). WUE L and WUE are driven by physiological adjustments such as changes in canopy conductance, but also by biophysical drivers such as increases in transpiration per unit of stomatal conductance at increasing D. This is especially true for WUE, as both the residual transpiration due to stomatal leakiness and leaf carbon losses due to reduced photosynthesis under limiting environmental conditions scale with stand leaf area. Accordingly, our results suggest that an increase in D by 48% would reduce the stimulating effect of eC a on whole-tree WUE to 20%, less than half when considering eC a alone.

| CON CLUS ION
A simulation model was used to analyze the potential sensitivity of a E. saligna plantation to changes in C a , D, precipitation, and access to deep soil water reserves. Best-case scenario (i.e., no nutrient limitation) model projections suggest that eC a will drive E. saligna responses to climate change, resulting in an increase in stand GPP and growth. Our results also confirm the importance of accounting together, both, soil and atmospheric aridity increases when projecting climate change impacts on plantations and forests. Thus, a 48% increase in D will limit CO 2 fertilization in a similar magnitude as a 20% reduction in precipitation, under our assumptions. Those two drivers combined are projected to halve the eC a fertilization effect on GPP, and also reduce WUE gains at the stand scale. The ability of E. saligna to access deep soil water reservoirs is projected to mitigate stress from increasing aridity, potentially allowing it to take advantage of eC a conditions. Nevertheless, we expect E. saligna to be increasingly reliant on such deep soil water in the future, linking its persistence to the still uncertain deep soil water reservoir dynamics. Daniel Nadal-Sala and Nadine K. Ruehr were funded by the Emmy Noether Program (RU 1657/2-1). We thank Dr. Rüdiger Grote for his comments on an early version of the manuscript. We thank two anonymous reviewers for their insightful comments. Open access funding enabled and organized by Projekt DEAL.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N S
Daniel Nadal-Sala wrote the initial version of the manuscript, de-

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.