Predicting resilience through the lens of competing adjustments to vegetation function.

There is a pressing need to better understand ecosystem resilience to droughts and heatwaves. Eco ‐ evolutionary optimization approaches have been proposed as means to build this understanding in land surface models and improve their predictive capability, but competing approaches are yet to be tested together. Here, we coupled approaches that optimize canopy gas exchange and leaf nitrogen investment, respectively, extending both approaches to account for hydraulic impairment. We assessed model predictions using observations from a native Eucalyptus woodland that experienced repeated droughts and heatwaves between 2013 and 2020, whilst exposed to an elevated [CO 2 ] treatment. Our combined approaches improved predictions of transpiration and enhanced the simulated magnitude of the CO 2 fertilization effect on gross primary productivity. The competing approaches also worked consistently along axes of change in soil moisture, leaf area, and [CO 2 ]. Despite predictions of a significant percentage loss of hydraulic conductivity due to embolism (PLC) in 2013, 2014, 2016, and 2017 (99th percentile PLC > 45%), simulated hydraulic legacy effects were small and short ‐ lived (2 months). Our analysis suggests that leaf shedding and/or suppressed foliage growth formed a strategy to mitigate drought risk. Accounting for foliage responses to water availability has the potential to improve model predictions of ecosystem resilience.


| INTRODUCTION
Climate projections suggest trees will soon encounter unprecedented extreme hydro-climatic conditions, which opens up important questions about their resilience.The ability of trees to maintain and recover function in response to disturbance preconditions their chances of avoiding dieback, die-off, and ultimately, local extinction.
Changes in the intensity, frequency, and timing of climate extremes (Masson-Delmotte et al., 2021) are already putting tree species at risk (Brodribb et al., 2020;McDowell et al., 2020).Droughts, for instance, have been implicated as one of the main drivers of an increasing trend in global forest mortality (Allen et al., 2015), and more widespread, frequent, and severe drought episodes are expected in the future (Ault, 2020;Masson-Delmotte et al., 2021;Ukkola et al., 2020).To gain insight into the future resilience of tree species, we need models capable of predicting functional adjustments as ecosystems undergo global change.In that context, process-based optimization approaches that shave away empirical assumptions and reduce reliance on detailed parameterization are gaining traction, as they could increase the predictive power of global vegetation models (Fisher & Koven, 2020;Harrison et al., 2021).
In land surface models (LSMs), two recently adopted optimization approaches are advancing the representation of dynamic ecosystem fluxes.First, theory around the optimal regulation of leaf nitrogen investment has been used to capture temporal changes in photosynthetic capacity, improving net/gross primary productivity predictions (Caldararu et al., 2020;Haverd et al., 2018).Second, substantial improvements in simulated canopy gas exchange have been achieved by optimally balancing photosynthetic gains and the hydraulic costs of transporting water between the roots and the leaves (Eller et al., 2020;Sabot et al., 2020).However, as much as optimality principles hold promise for future LSM developments, they are not completely exempt from shortcomings.Chiefly, the choice of fitness proxies over which to optimize function and the selection of relevant optimization timescales are critical to model outcomes.
In reality, it is unlikely plants can optimize all resource acquisition and develop strategies to cope with stress factors simultaneously (Stearns, 1989).Optimality principles that are successful in isolation may therefore not be successful when combined with competing principles that operate on different timescales, or when rendered "suboptimal" by environmental disturbances.The few studies that have considered multiple functional adjustments together have typically relied on a unique or dominant overarching optimization objective, such as to maximize long-term carbon profits (see Potkay et al., 2021;Schymanski et al., 2009Schymanski et al., , 2015)).For example, although long-term optimization objectives may be suitable for predictions of decadal carbon and/or water balance, they might neglect nonlinear interactions between processes, leading to unrealistic short-term predictions (see Kerkhoff et al., 2004).To our knowledge, there has not been an attempt to examine whether competing optimality principles can operate interactively and consistently across timescales of application, nor when system legacy is involved (i.e., lingering effects from repeated disturbance).
Multiple evidence streams-satellites (Wigneron et al., 2020), eddy covariance (Schwalm et al., 2017), tree-rings (Peltier et al., 2021;Vanoni et al., 2016), and manipulation experiments (Zweifel et al., 2020)-have shown that drought impacts on forest carbon stocks and tree growth can last several years after droughts have concluded.One of the mechanisms most often invoked to explain these legacies is hydraulic damage.Hydraulic damage is sustained from incomplete and/or lagged recovery of xylem function when cavitation is significant in trees (Klein et al., 2018).Hydraulic damage can cause leaf drop and/or hinder foliage growth, and lead to canopy dieback or even tree mortality in severe cases (McDowell et al., 2008).
Critically, hydraulic damage continues to impede water flow between the soil and the canopy after soil water is recovered.Integrating this type of system legacy with optimality approaches is a challenge we must tackle to ensure our models capture future climate risk to the vegetation appropriately.
In this study, we combine nitrogen-based and hydraulics-based plant optimality theories within a single LSM framework (Sabot et al., 2020), while also incorporating a representation of legacy effects from hydraulic damage.Doing so allows us to investigate interactions between competing optimization theories affected by disturbance, at the forefront of LSM development.
We use observations from the Eucalyptus Free Air CO 2 Enrichment (EucFACE) experiment to parameterize and assess our model.The EucFACE experiment is set in a native, mature Eucalyptus woodland that is nutrient-limited and that strongly responds to water availability.For example, Duursma et al. (2016) demonstrated that leaf flushing events occurred shortly after significant rainfall, whereas the site leaf area index (LAI) declined steadily during periods of low rainfall.The period considered here (2013-2020) encompasses frequent droughts and heatwaves, including a multi-year drought from mid-2017 onwards.
Recurring droughts and heatwaves may not only affect the hydraulic transport system, but also leaf nitrogen investments in the photosynthetic apparatus and hence photosynthetic capacity (Damour et al., 2008;Hikosaka et al., 2006;Vaz et al., 2010).Theory and experimental evidence also suggest that, besides inducing some leaf nitrogen reallocation (Ainsworth & Rogers, 2007), increasing atmospheric [CO 2 ] could alleviate plant hydraulic stress during drought (De Kauwe et al., 2021;Jiang et al., 2021).However, we do not know Our specific objectives are: (i) to infer the timescales over which vegetation function adjusts to maximize resource investment and/or resilience; (ii) to test whether contrasting optimization approaches can be applied together, including when hydraulic legacies from climate extremes affect optimization outcomes; (iii) to explore the relationship between drought-driven foliage dynamics and plant hydraulic function; and (iv) to understand whether (i), (ii), and (iii) are affected by elevated [CO 2 ].

| Overview of the model
The model described in Sabot et al. (2020) was adapted from the Community Atmosphere Biosphere Land Exchange LSM (CABLEv2.0; De Kauwe et al., 2015).In brief, our model accounts for: canopy interception; radiation scattering into sunlit and shaded components of a single-layer canopy; within-canopy micrometeorology; optimal canopy gas exchange; plant hydraulic processes; soil evaporation; and soil hydrology over six discrete layers.In this study, the canopy interception (Methods S1) and soil hydrology (Methods S2) schemes are revised to better reflect observations at EucFACE (Gimeno et al., 2018;see below).We also extend the model by embedding a scheme that optimizes leaf nitrogen allocation, thus varying photosynthetic capacity, and by incorporating hydraulic legacies alongside our existing canopy gas exchange optimization scheme (see below).
Sunlit and shaded components of the canopy are each assumed to optimally balance relative photosynthetic gains and relative hydraulic costs through their stomates: where A n (µmol m −2 s −1 ) is the net rate of carbon assimilation; A i max , (µmol m −2 s −1 ) is the instantaneous maximum A n there could be considering soil moisture limitations; the maximum root-to-leaf hydraulic conductance attenuated by soil moisture stress at a specific instant; hydraulic conductance evaluated at a given leaf water potential Ψ l (MPa); and k crit (mmol m −2 s −1 MPa −1 ) is the critically low threshold of hydraulic conductance that indicates xylem failure.
Except k crit (constant), all the elements that makeup Equation (1) vary on an instantaneous basis.As the elements in Equation ( 1) are computed from system of equations that connects photosynthetic, stomatal, and plant hydraulic processes (Methods S3 and S4), all these processes are also solved for when Equation ( 1) is solved for.

| Leaf nitrogen allocation
We draw on work that links nitrogen stoichiometry to the coordination of photosynthetic capacities (Chen et al., 1993;Field, 1983;Maire et al., 2012;Medlyn, 1996) to: (i) optimize the amount of leaf nitrogen involved in photosynthesis (N p ; mol m −2 ) and its distribution to photosynthetic compounds; whilst also (ii) coordinating key photosynthetic parameters, that is, V cmax25 (µmol m −2 s −1 ), the maximum carboxylation rate (V cmax ) at 25°C, J max25 (µmol m −2 s −1 ), the maximum rate of electron transport (J max ) at 25°C, and R d25 (µmol m −2 s −1 ), the day leaf respiration (R d ) at 25°C.N p is divided into thylakoid nitrogen and soluble protein nitrogen compounds that link to V cmax25 and J max25 (Medlyn 1996; Methods S5), so our approach allows N p and the photosynthetic parameters to co-vary with environmental conditions, even in the absence of an explicit nitrogen cycle.
Our model seeks to satisfy: (3) Both Σ ante and N p vary over a chosen time period (see Section 2.3).Importantly, the N p and photosynthetic parameters optimized through Equation ( 2) do not impose model co-limitation through time because "future" environmental conditions will differ from Σ ante .Equation ( 2) is bound by the inequalities given in Equations ( 3) and (4).Equation (3) signifies that the new optimized N p and photosynthetic parameters must a posteriori increase A n over the conditions given by Σ ante , compared to the previous time period N p (Np t−1 ).Equation (4) avoids optimization "overshoots" in the ratio of J V : Medlyn, 1996) by approximating a triose phosphate utilization (TPU) limitation on the objective function (Collatz et al., 1991).Note, however, that no TPU limitation is applied in the photosynthesis model (Methods S3).

| Hydraulic legacies
Model representations of plant hydraulics often assume perfect xylem repairs or complete xylem regrowth upon resumption of water supply following cavitation (but see Venturas et al., 2018).As such, the same plant vulnerability curve is used to describe hydraulic conductance (k; mmol m −2 s −1 MPa −1 ) decline and recovery with changes in water potential (Ψ; MPa), both before and after significant embolism: where k max (mmol m −2 s −1 MPa −1 ) is the maximum root-to-leaf hydraulic conductance; and f (Ψ) vc is the hydraulic vulnerability curve for the perfectly recoverable xylem, which captures the increase in percentage loss of conductivity (PLC) as Ψ decreases (Methods S4).
We account for delayed and imperfect xylem recovery from past embolism by substituting f (Ψ) v by a "damaged" vulnerability curve (Lu et al., 2020): where  diffuse canopy transmittance and monthly leaf litter production (Duursma et al., 2016).Figure S1 shows major model drivers.
The model was parameterized as shown in Table 1; parameter values were either retrieved from the literature or estimated from field observations (Methods S6).Other model parameters that were kept to their default values are shown in Tables S1 and S2.

| Timescales of adjustments to function
We ran 36 simulation experiments to identify the timescales over which competing adjustments to vegetation function interact most consistently.Our default model configuration assumes the optimization of canopy gas exchange only, with no change in leaf nitrogen and photosynthetic parameters over time and no hydraulic legacies.
Other configurations optimize canopy gas exchange, plus leaf nitrogen allocation and/or hydraulic legacies (Figure 1).
The timescales relevant to optimizing canopy gas exchange are not examined, as previous work indicated less than 3-hourly intervals were appropriate (Sabot et al., 2020), and half-hourly intervals are therefore used.Five model configurations consider the optimal regulation of leaf nitrogen (N opt ) every 7th, 14th, 21st, 28th, or 35th day.Another five configurations assess hydraulic legacies (H leg ) every 14th, 30th, 60th, 90th, or 180th day.The remaining 25 configurations enable joint variations in leaf nitrogen and hydraulic legacies, by combining the timescales previously listed for each process.The following matrix summarizes the model configurations, that is, the processes accounted for in addition to the optimization of canopy gas exchange and their associated timescales of application: wheredenotes the default model configuration, and N opt7 and H leg14 stand for N opt every 7th day (or N opt = 7 days) and H leg every 14th day (or N leg = 14 days), respectively.
The 36 model configurations were evaluated against observations of daily transpiration (E; see Section 2.4) from the aCO 2 rings only, hence the model behaviour is emergent at the eCO 2 rings.
Specifically, the configurations were ranked depending on: (i) their relative variability; (ii) their similarity with the observations; and (iii) their accuracy.Relative variability was determined by the ratio of modelled to observed sample standard deviation (cf. the approach of Best et al., 2015).Similarity was measured by a skill score of overlap between modelled and observed binned distributions (Perkins et al., 2007).Model accuracy was gauged via the Mean Absolute Scaled Error (Hyndman & Koehler, 2006), which establishes whether model predictions are more skilful than one-step T A B L E 1 Model parameterizations specific to the aCO 2 (R2, R3, and R6) and eCO 2 (R1, R4, and R5) rings at EucFACE.
l where a abs is the absorptance of photosynthetic photon flux density (Yang et al., 2020) and τ l is leaf transmittivity (0.05 in our model).b Values estimated from substrate-saturated carboxylation rates, Rubisco site counts, and totals of leaf nitrogen allocated to Rubisco reported by Sharwood et al. (2017) for upper canopy leaves of Eucalyptus globulus.c Values derived by Peters et al. (2021) from hydraulic vulnerability curves constructed by benchtop dehydration for E. tereticornis growing in the vicinity of EucFACE.d This corresponds to a 95% loss of hydraulic conductance, representing the maximum recoverable water stress in a number of species (Brodribb & Cochard, 2009;Sabot et al., 2020).*Parameter values shared across rings appear on the right-hand-side.**In model simulations that do not consider leaf nitrogen allocation (e.g., our default model configuration), these parameters are kept at their t 0 value through time.
forecasts of the observations.Mathematical formulations are given in Methods S7.
External contributions to change in modelled V cmax25 and J max25 from the LAI, PPFD, T a , D a , C a , and the modelled soil moisture (soil water potential, Ψ s ) were also analyzed for the best-ranked model configuration.To do so, we used Dominance Analysis (Azen & Budescu, 2003) which generalizes R 2 to all possible combinations of predictors before decomposing it for each predictor (i.e., environmental driver).Here, the method serves as an alternative to linear regression, to more meaningfully identify the relative contributions of each environmental driver to the variance in modelled V cmax25 and J max25 .

| Foliage response to drought
Finally, we carried out an experiment to infer a link between foliage response to drought (i.e., leaf abscission and/or a lack of new leaf production, without distinction between the two processes) and plant hydraulic status.We asked: how different would plant hydraulic status have been, had LAI not declined during drought?For each ring, we computed the average LAI phenology from 2013 to 2020, and we used it to force the best-ranked model configuration.The resulting predicted PLC (from the average phenology) was compared to the PLC predicted using the observed LAI, and the difference between the two simulations shed light on the hydraulic damage avoided through foliage adjustment.
We also estimated the carbon allocated to leaf growth in actual LAI versus average LAI simulations, using a data-derived, carbon budget framework specific to EucFACE (Extended Data fig.7 in Jiang et al., 2020b).We estimated net primary productivity (NPP) by applying a conversion factor accounting for autotrophic respiration to our model simulations of gross primary productivity (GPP), and the carbon allocated to leaf growth was then estimated from NPP using site-specific carbon allocation fractions (Table S3).
By probing the carbon benefits of increasing hydraulic risk by increasing LAI and, conversely, the carbon risks of drought-induced LAI drops, we gain a fuller understanding of how LAI dynamics relate to plant carbon and water status.

| Observational data sets
Observed V cmax25 , J max25 , and R d25 were used to parameterize the model (Table 1; see Methods S6), and to evaluate the temporal changes in photosynthetic capacities predicted by the optimal leaf nitrogen allocation scheme.The observed estimates were obtained from analyzing 323 repeated light-and temperature-controlled photosynthesis-CO 2 response (A-C i ) curves, using the "plantecophys" R package (Duursma 2015).A-C i curves were measured on mature (>90 days old) and newly flushed (<90 days old) leaves of the upper canopy of three to four dominant or codominant trees in each ring, from February 2013 to February 2020 (16 measurement campaigns in spring, summer, and autumn, extending previous compilations by Ellsworth et al., 2017 andWujeska-Klause et al., 2019).
Observations of foliar nitrogen concentration on either a mass or an area basis (i.e., mg g −1 or mg m −2 ) are included for a qualitative comparison with the modelled N p .Nitrogen concentrations were measured for 236 upper canopy leaf samples that were collected from the three to four (co-) dominant trees of each ring, from February 2013 to February 2018 (11 measurement campaigns; Gherlenda et al., 2016a;Wujeska-Klause et al., 2019).
Morning (9:30-11:30 h) and afternoon (13:00-15:00 h) Ψ l measured from October 2012 to January 2014 (Gimeno et al., 2016) were used to parameterize r k (Table 1; Methods S6) and to evaluate F I G U R E 1 A schematic representation of the interactions between the components of our model that simulate adjustments to vegetation function.Our default model configuration optimizes canopy gas exchange only, on a half-hourly basis.Other model configurations can either account for hydraulic legacies (H leg ), or optimize leaf nitrogen allocation (N opt ), or do both, utilizing information from the canopy gas exchange routine at a given frequency (e.g., 7 days' worth of information for N opt ), and thus feeding back into the canopy gas exchange routine over time.We used tree-level E estimates for parameterization of k max (Table 1; Methods S6).These daily E estimates were derived from tree sap velocities measured between January 2013 and September 2014, using the heat pulse compensation technique in the three to four (co-)dominant trees of each ring (Gimeno et al., 2018).The measurements were converted to tree-level E using ring-specific sapwood area per unit ground area ratios, as estimated from trees neighbouring the experimental rings.For comparison with our model simulations at the ring-level, we multiplied the tree-level E estimates by their respective ring's LAI.
Finally, root-zone soil moisture was estimated from two sets of observations described in Gimeno et al. (2018).First, daily measurements were retrieved from January 2013 to December 2019, using theta probes to a depth of 30 cm at each ring.Second, soil moisture was monitored at 25 cm intervals from 25 to 150 cm depth, and at 50 cm intervals from 150 to 450 cm depth, every 10-21 days (less frequently in 2017) between January 2013 and July 2019, using neutron probes at two locations per ring.We weighted the observed soil moisture by root fraction at depth (see Methods S2 and S6) to obtain an estimate of root water access.

| Timescales of adjustments to function
To understand whether competing optimization schemes, perturbed by hydraulic legacies, varied predictably and coherently, we examined differences between the default model configuration (optimization of canopy gas exchange only) and the 35 other configurations that combine functional adjustments.We focussed on analyzing model configuration differences on average water use efficiency (WUE; the ratio of A n to E) because WUE integrates information on both the canopy carbon and water fluxes.
Additionally, we have a theoretical understanding of how WUE ought to behave across axes of variation in LAI, soil moisture, and atmospheric [CO 2 ], and this helps assess the predictability and coherence of our combined competing optimality criteria.Within-ring WUE changes (ΔWUE config relative to the default WUE; see annotations in Figure 2) were small, ranging between −4.8% and 1.8% across rings.At most rings, the greatest change was associated with combinations of short timescale adjustments: the combined 14-day H leg and 7-day N opt .Thus, the higher the frequency of adjustments, the less water-use efficient the canopy gas exchange scheme.This result is explained by the fact that impairments from short episodes of water stress have most bearing on high frequency (short timescale) adjustments to function.For example, a 3-day reduction of C i caused by soil moisture stress ought to reflect in a weekly adjustment but not in a monthly adjustment of N p .The contrary was true for R5 (ΔWUE config > 0) where the largest gain of WUE corresponded to combinations of low frequency (long timescale) adjustments (H leg = 180 days, In all cases, N opt drove change in WUE, as H leg 's contributions to ΔWUE config were small.Indeed, the effects from H leg varied alone (e.g., H leg = 14 days, at the bottom of each plot) added to the effects from N opt varied alone (e.g., N opt = 7 days, at the bottom of each plot) did not diverge from the effects of H leg and N opt varied jointly (e.g., H leg = 14 days, N opt = 7 days, at the top of each plot), as seen from the lack of difference between the vertical line markers (added singular effects) and the horizontal bars (interactive joint effects).
The 36 model configurations' ability to simulate daily E at the aCO 2 rings is reported in Figure 3.For the accuracy metric (see Methods S7), the lower the values, the better the performance; values <1 indicate that model predictions are more skilful than onestep-ahead forecasts using the preceding observations.N opt alone and N opt combined with the longest timescale of H leg resulted in the highest accuracy (first 10 absolute ranks <0.784).For the similarity and variability metrics (see Methods S7), the closer to 1 the better.
The performance was typically higher for the higher frequency H leg and N opt (>0.795 and >0.877, respectively), which operated more functional adjustments to short episodes of water stress than lower frequency H leg and N opt .
Quantile ranks were used to summarize skill across statistical metrics, considering a configuration's performance relative to that of all other configurations within each metric.Perhaps counter intuitively, the best quantile-ranked configurations could occupy low or average positions for the accuracy metric (see the annotations in Figure 3) because there were only small differences between the 1st and 30th ranked configurations (spanning less than half the yaxis).Overall, H leg alone only marginally improved on the default configuration, and we were unable to establish clear timescales of influence.By contrast, skill distinctly increased with the frequency of N opt optimization (higher skill for N opt = 7-14 days and lower skill for

| Leaf-and stand-level dynamics
Figure 4 shows that V cmax25 and J max25 were well captured (averages within ±5% of the observations) by the best model configuration (H leg = 30 days, N opt = 7 days), for both the aCO 2 (a and c) and eCO 2 (b and d) rings.Assessing the ability of the model to reproduce the observed dynamics was complicated by tree-to-tree variability in the measured data, particularly V cmax25 (e.g., ranging between 33.3 and   116.6 μmol C m −2 s −1 at R2 on 30 January 2014).Nevertheless, the model predicted the overall observed seasonal shift from low J max25 in spring/summer (low LAI; Figure S1) to high J max25 in autumn (LAI peaks in autumn/winter).The model also predicted the increase in the ratio of J V : In the springs of 2013 and 2016, LAI was very low (c.0.6 m 2 m −2 and below 0.5 m 2 m −2 in some of the rings; Figure S1) and the simulated V cmax25 was overestimated due to underestimated J V : max cmax 25 25 .The simulated V cmax25 appeared high in the springs and summers of 2017/2018 and 2019/2020 too; but there is a paucity of observations from mid-2017 onwards.Statistical analysis revealed LAI to be the main driver of seasonality in the modelled V cmax25 and J max25 .At most rings (R2 and R3, aCO 2 ; R1 and R4, eCO 2 ), >74% of the modelled variance in V cmax25 and J max25 could be explained by variation in LAI (contributing 54%) and the environmental forcing (D a 9%, PPFD 6%, other drivers <5% each).At R6 The effect sizes associated with different combinations of simulated adjustments to vegetation function, at the aCO 2 (a) and eCO 2 (b) rings.The horizontal bars reflect magnitude differences on average water use efficiency (WUE; the ratio of net carbon assimilation to transpiration) between the default model configuration and the model configurations that consider additional functional adjustments.There are 35 such configurations, accounting for: hydraulic legacies (H leg ; singular effects), leaf nitrogen allocation (N opt ; singular effects), and combinations of H leg and N opt (joint effects) over their respective timescales of application (days; y-axis on the left side of the figure).The vertical black lines superposed on the joint effects show the corresponding singular H leg and N opt effects added together (i.e., vertical lines on both sides of the horizontal bar for the joint H leg = 14 days, N opt = 7 days correspond to the singular H leg = 14 days added to the singular N opt = 7 days).There is no x-axis or scale on the figure: the horizontal span dedicated to each ring corresponds to its magnitude effect size relative to that of all other rings.The sign of the effects at each ring is indicated by an arrow (↘ for a decrease vs. ↗ for an increase in model configuration WUE), with asterisks marking instances where the sign of the effect opposes all other effect signs within a ring (e.g., H leg = 180 days at R6).At each ring, the largest percentage change in model configuration WUE is indicated over its corresponding functional adjustments.Finally, within each panel, the rings are ordered in increasing order of soil moisture × leaf area index from left to right.[Color figure can be viewed at wileyonlinelibrary.com] (aCO 2 ) and R5 (eCO 2 ), both LAI and soil moisture were typically higher than in the other rings and the environmental drivers explained less of the variance in the modelled V cmax25 and J max25 (c.61%).The contribution of LAI was not as substantial (31%), whilst that of soil moisture (Ψ s ) increased to 17%.
At the eCO 2 rings, the down-regulation of V cmax25 averaged to 9.6% in the observations versus 11.2% in the simulations (Table 2).In the model, this down-regulation of V cmax25 partly originated from the initialization of N P using ring-specific Vcmax t 25 0 and Jmax t 25 0 (see Table 1 and Methods S6).When N P was initialized using the average Vcmax t 25 0 and Jmax t 25 0 from across the aCO 2 rings (Figures S3 and S4), our leaf nitrogen allocation scheme only predicted a 5.5% reduction in V cmax25 at the eCO 2 rings.This underestimation of the eCO 2 down- regulation of V cmax25 was explained by the model's overestimation of V cmax25 at low LAI, i.e., at R1 and R4 (see Table 2).At R5, LAI was comparatively high, and it did not drive V cmax25 as strongly; consequently, the model successfully down-regulated V cmax25 .
Notwithstanding initialization, the dynamics of N P (Figures S2 and S4) aligned with those of the observed leaf nitrogen concentration on an area basis, albeit displaying one order of magnitude less variability.Finally, similarly to the observed leaf nitrogen concentration on a mass basis, N P was lower at the eCO 2 rings than at the aCO 2 rings.−4.8% to +1.8% GPP across rings (akin to the within-ring WUE changes from Figure 2) and E virtually unchanged (−0.5% to +1.8%; r 2 ~0.7 in each case).Magnitude differences in simulated GPP between the default and the best model configurations were the F I G U R E 3 The performance ranks of the 36 different model configurations at the aCO 2 rings (R2, R3, and R6).The model configurations are evaluated against the observed daily transpiration using three statistical metrics: (i) the ratio of modelled to observed sample standard deviation (Variability metric; diamonds), (ii) a skill score of overlap between modelled and observed binned distributions (Similarity metric; circles), (iii) the Mean Absolute Scaled Error (Accuracy metric; triangles).The spans of the y-axes reflect the range of variation within a given statistical metric relative to that of other metrics, and the y-axes are arranged in either descending or ascending order for better performance to appear at the bottom of each axis.From left to right, the horizontal arrow at the bottom of the figure shows the least (default) to most (H leg = 30 days, N opt = 7 days) skilled model configurations, as given by quantile ranks summarizing all aspects of performance (i.e., each configuration's performance relative to that of all other configurations within each statistical metric).The five most skilled (or better ranked) configurations are annotated on the figure, next to their respective performance in each statistical metric.[Color figure can be viewed at wileyonlinelibrary.com]F I G U R E 4 Timeseries of the V cmax25 (a and b) and J max25 (c and d) predicted by the best model configuration (H leg = 30 days, N opt = 7 days) and compared with field observations at the aCO 2 (a and c) and eCO 2 (b and d) rings.Orange lines show the average simulations across rings (e.g., average of R2, R3, and R6 for the aCO 2 rings) and yellow shadings the range of simulations.Box and whisker plots summarize the observed distributions of V cmax25 and J max25 (line, median; box, interquartile range; whiskers, quartiles ± 1.5 times the interquartile ranges) during measurement campaigns, with each observation displayed via a dot.Fixed parameter values appear in each panel, these refer to the ring-specific parameterizations used by the default model configuration, as well as in the initialization of the leaf nitrogen allocation scheme.Note: The aCO 2 (R2, R3, and R6) and eCO 2 (R4, R1, and R5) rings appear ordered from smallest to highest average soil moisture × leaf area index, as in Figure 2. The percentage down-regulation of V cmax25 at eCO 2 is shown in italic below the average aCO 2 | eCO 2 pairing.Two sets of simulations conducted with the best model configuration (H leg = 30 days, N opt = 7 days) appear in the table: (i) simulations initialised by ring-specific parameterizations of V cmax25 t 0 (and J max25 t 0 ); (ii) additional simulations initialised by the average V cmax25 t 0 (and J max25 t 0 ) from across the aCO 2 rings applied to all rings.For context, the parameterizations of V cmax25 t 0 appear above the modelled values of V cmax25  5e,f).This lack of sustained hydraulic legacies could be assumed to have arisen from an overestimation of soil moisture, especially after 2017 (Figure S5).Nevertheless, hydraulic legacies remained small when soil moisture was prescribed from observations (Figure S6).Possible explanations for these limited hydraulic legacies and the ensuing lack of effect on the surface fluxes are: (i) that changes from legacies compensated each other (Figure S7) because diel and day-to-day fluxes were highly dynamic; and/or (ii) that legacies were realized though extra-xylary adjustments that allowed the trees to maintain function, for instance, adjustments to LAI.In each of these insets, the black line corresponds to the 1:1 of the observations and coloured lines correspond to Q-Q plots of the observed to simulated transpiration.In the absence of direct observations, the simulated PLC is assessed by proxy, by comparing distributions of observed and simulated Ψ l (morning and midday data) from the collection campaign dates marked by black dots in (e) and (f).These distributions, fitted via kernel density estimation, and their associated modes (μ) appear in the insets of (e) and (f).In all main panels (i.e., excluding the insets), the data have been smoothed with a 90-day moving window to aid visualization.Higher LAI than the average phenology (i.e., Δ rel LAI > 0) generally resulted in a net PLC increase (i.e., percentage point increase) of 1.5% on top of a background PLC of 4.2% at the aCO 2 rings.Higher LAI also led to a 2.9% net PLC increase on top of a background PLC of 3.6% at the eCO 2 rings.In mid-2014, after a multi-year peak in LAI at all rings (Figure S1), up to 45% and 78% net increases in PLC were simulated from background PLCs < 5% at the aCO 2 and eCO 2 rings, respectively.This particularly high PLC may explain the ensuing low observed LAI in 2015 (Figure S1).

T A B L E 2 A comparison of observed and simulated maximum carboxylation rates at 25°C (V
Increased PLC from the model forced with LAI higher than the average phenology was accompanied by increasing cumulative carbon storage until either 2016 (67 gC m −2 ; eCO 2 rings) or 2017 (63 gC m −2 ; aCO 2 rings).After those years, the LAI declined greatly (Figure S1), reducing carbon uptake and storage (Figure 6c,d).Still, total accumulated leaf carbon stocks were higher in 2019 than they would have been if the canopy had remained as per the average phenology (Δ y A fcum > 0 in the insets of Figure 6c,d), highlighting the beneficial nature of the observed LAI dynamics.Additionally, forcing the model with LAI lower than the average phenology deducted c.
1.8% from a background PLC of 5.8% at the aCO 2 rings, and c. 2.5% from a background PLC of 5.3% at the eCO 2 rings.In mid-2018, as much as 73.3% and 85.4% net reduction in PLC were simulated from respective background PLCs of 81.0% and 89.8% at the aCO 2 and eCO 2 rings.
On the whole, both net PLC increases from higher LAI and net PLC decreases from lower LAI were more important at the eCO 2 than at the aCO 2 rings (Figure S8), though the relationship was asymmetrical and favoured the net PLC decreases.This result implies higher eCO 2 ratios of ΔPLC:ΔLAI that were further enhanced by drought and/or low LAI.S6).
In our simulations, high canopy PLC subsided rapidly, as hydraulic legacies and their impacts on the transpiration were small.This outcome partially stemmed from our parameterization of the ratio of xylem recovery from embolism r k (Methods S6) which, at c. 0.9 (Table 1), implies a high ability to recover lost hydraulic conductance.A recent study of hydraulic legacies in two Eucalyptus hybrids (Saunders & Drew, 2021) pointed to a possible trade-off in hydraulic strategy between vulnerability to cavitation and recovery from embolism, whereby the hybrid characterized by higher hydraulic vulnerability (faster stomatal closure) recovered conductance faster than its counterpart that better resisted the formation of embolism.This mechanism could be envisioned for the E. tereticornis examined here, as their P 50 of −4.07 MPa confers less resistance to cavitation than the P 50 of other trees in the vicinity of EucFACE (e.g., E. fibrosa: −5.82 MPa, E. moluccana: −5.74 MPa, Melaleuca decora: −5.68 MPa; Peters et al., 2021).
A drawback of our approach to modelling legacy is that it is sensitive to the parameterization of r k .However, xylem recovery capacity itself may be sensitive to drought severity (Gauthey et al., 2022), and so may the value of r k .Therefore, even though we derived r k from site measurements over a period that comprised notable dry spells (October 2012-November 2013), we cannot exclude the possibility that r k may be different during a longer, more severe drought, for example, if it were derived using data from the 2017-2020 drought period.More generally, reliance on speciesspecific parameterizations of r k would increase model complexity without always impacting hydraulic capacity and/or canopy fluxes.
Meanwhile, there is mounting evidence of interconnected droughtand heat-induced legacies on nonstructural carbohydrates (NSCs) and hydraulic function (Duan et al., 2022;Poyatos et al., 2013;Ruehr et al., 2019).The key to LSMs accounting for differentiated legacies from climate extremes (i.e., depending on species and/or location) might lie in a mechanistic link between NSCs and plant hydraulics.
Representing NSCs (see e.g., Jones et al., 2020) alongside photosynthetic, stomatal, and hydraulic processes in LSMs would provide the opportunity to elucidate such a link, by way of observation-driven model experiments.

| Leaf nitrogen allocation and photosynthetic capacities
Previous tests of optimal leaf nitrogen allocation schemes in LSMs (Caldararu et al., 2020;Haverd et al., 2018) had not considered water-limited, low LAI ecosystems such as EucFACE (long-term LAI of c. 0.8 m 2 m −2 across rings), where GPP is predominantly limited by the ability to capture light (J max25 ) rather than to fix CO 2 (V cmax25 ).
Ecosystems in dry environments play a noteworthy role as drivers of interannual variability in the global carbon cycle (Ahlstrom et al., 2015;Poulter et al., 2014), and so, the dynamics of their limiting rates to photosynthesis warrant further attention.
Our optimal leaf nitrogen allocation scheme was generally able to resolve the dynamics of V cmax25 and J max25 (Figure 4), but it overestimated V cmax25 when/where LAI was very low (c. 0.6 m 2 m −2 ; see Figure S1).The model also failed to predict the eCO 2 downregulation of V cmax25 where the LAI was very low (i.e., R1 and R4 vs. R5 where it was successful; see  (Sharwood et al., 2017).The simulated decline in N P at the eCO 2 rings was also in qualitative agreement with the observations of leaf nitrogen concentration on a mass basis (cf.Crous et al., 2019;Wujeska-Klause et al., 2019), but not with the observations of leaf nitrogen on an area basis (i.e., accounting for variation in leaf mass per area) that were more often than not higher or level with the aCO 2 rings.This mismatch between the behaviour of the leaf nitrogen concentration on a mass versus an area basis, along with the sensitivity of our modelled V cmax25 to the initialization of N P , suggests drought/ LAI-driven reallocations of plant nitrogen at the eCO 2 rings consistent with the findings of Crous et al. (2019).Exploring these reallocation patterns is beyond the scope of the current study, but could be done exploiting theoretical frameworks that link nitrogen and water balance to LAI (e.g., McMurtrie et al., 2008).

| Why did combined model adjustments to function reduce water-use efficiency?
Many studies that consider optimality principles allow plants to reach optimum points without accounting for potential 'suboptimality' in behaviour (Caldararu et al., 2020).A novel aspect of our study is that we formalized our optimality criteria such that they interact whilst not serving the same purpose, and we limited optimal behaviour via a representation of hydraulic legacies.
Our two optimality criteria interacted coherently, including when paired with hydraulic legacies.Their combinations showed predictability with respect to environmental factors (soil moisture, LAI, [CO 2 ]), and behavioural coherence over their timescales of application (Figures 2 and 3).Their propagation to the stand-scale was small (i.e., <5% average decline in GPP), but of similar magnitude to the findings of other studies that considered seasonal changes in photosynthetic capacities (Bonan et al., 2011;Medvigy et al., 2013).
Here, it is also worth noting that our comparison point (i.e., the default model configuration) is a canopy gas exchange optimization approach that has been shown to significantly outperform standard LSM implementations (Sabot et al., 2020).
Contrary to leaf-level photosynthesis where Rubisco limitations on photosynthesis prevail throughout the day, our stand-scale GPP simulations were chiefly limited by RuBP regeneration, owing to relatively low ratios of J V :  et al., 2020).Due to this RuBP regeneration limitation, J max25 is a much more relevant predictor of GPP than V cmax25 at EucFACE, particularly when LAI drops during drought (at the end of 2013, from 2018 onwards; see Figures 4 and 5).Decreased GPP from the downregulation of J max25 translated into quasi-equivalent changes in WUE, because our optimization model sought to maintain E (preventing excessive reductions of C i ) to prevent the further decline of GPP during low LAI/drought periods.Unlike the transpiration simulated here, a 'standard' LSM (e.g., CABLEv2.0) would have predicted E to decline with the simulated reduction in GPP, by keeping the WUE quasi-constant.
One potentially influential physiological association that we did not examine is the link between k max and V cmax25 .This association was explored by Sabot et al. (2020) over decadal timescales and their results broadly supported a link between k max , average V cmax25 , and long-term climate.In theory, k max should vary with the length of the hydraulic path between the roots and the leaves, so it ought to adjust (sub-)seasonally as photosynthetic investments shift to different areas of the canopy (e.g., from the top of the sunlit canopy to lower shaded leaves; Peltoniemi et al., 2012).Still, whether k max adjusts in coordination with V cmax25 and/ or J max25 over timescales as short days or weeks, and whether that is to maximize GPP or to preserve WUE remain open questions.

| Implications at elevated [CO 2 ]
At eCO 2 , our best model configuration (i.e., high-frequency combinations of optimality criteria) predicted a 13.1% increase in GPP (+3.5 percentage points compared to our default model), which compared well with a model-data fusion estimate of a 12% increase in GPP over the 2013-2016 period (Jiang et al., 2020b).
Both measurements of leaf-level stomatal conductance (Gimeno et al., 2016) and our model predictions of canopy conductance agree on a eCO 2 -induced decline of the order of 20% under favourable conditions (mid-morning, no water stress).At the same time, both observations of tree-level transpiration (Gimeno et al., 2018) and our modelled stand-scale transpiration were not reduced by eCO 2 .Instead, they were slightly higher under eCO 2 .This apparent contradiction is explained by a higher average soil-to-canopy maximum hydraulic conductance under eCO 2 (k max 1.44 mmol m −2 s −1 MPa −1 ) than aCO 2 (k max 1.16 mmol m −2 s −1 MPa −1 ), accompa- nied by higher Ψ l (c.0.15-0.22MPa higher), which may indicate hydraulic acclimation and/or reflect the heterogeneity of the site, for example, differences in soil textures and water retention (Table 1).
Observational evidence of hydraulic acclimation to eCO 2 is equivocal (De Kauwe et al., 2021;Domec et al., 2017).Nevertheless, physiological adjustments of the same nature as those captured here (i.e., higher k max , higher E, higher Ψ l ) have also been observed for Eucalyptus grown in a glasshouse and exposed to dry-downs and eCO 2 (Jiang et al., 2021).Moreover, both higher k max and the higher eCO 2 ratios of ΔPLC:ΔLAI simulated in our "foliage response" experiment are consistent with the eco-hydrological equilibrium hypothesis (Sperry et al., 2019;Yang et al., 2018).

| Foliage response to drought
We wanted to understand the mechanistic causes for the adjustments in LAI that are observed with changes in rainfall at EucFACE (Duursma et al., 2016).To do so, we linked LAI, carbon storage, and hydraulic status in a model, which has seldom been attempted before (but see Trugman et al., 2018;Xu et al., 2016).This approach contrasts predominantly statistical efforts to identify legacy effects at the landscape-and ecosystem scales (see Kannenberg et al., 2020) for a review), similar to most work on ecosystem resilience to repeated droughts and heatwaves (e.g., Mitchell et al., 2014;Neumann et al., 2017).
2015 was the only year for which we could potentially link low LAI to a compromised hydraulic system (i.e., marked simulated PLC in mid-2014).Coincidentally, an outbreak of senescence-feeding psyllids (Cardiaspina fiscella) caused major defoliation at EucFACE in September 2014 (Gherlenda et al., 2016b).More generally, though, we found that foliage responses to drought (i.e., leaf abscission and/or a lack of leaf growth) prevented the loss of hydraulic function (Figure 6).Droughtinduced drops in LAI are often perceived as a pervasive consequence of water stress that negatively impacts carbon stocks (e.g., Anderegg et al., 2015).But at EucFACE, under water-limited conditions, an alternative interpretation is that this strategy could be akin to "ecological memory" benefiting the trees.Foliage growth was sufficient in the relatively wetter years to buffer carbon losses associated with lower LAI during the drier years (2013-2016/2017 vs. 2017-2020; see Figure 6).These findings concur with Eucalyptus canopy dynamics observed in the field (Battaglia et al., 1998;Pook, 1985;Whitehead & Beadle, 2004), but explicit linkage of observed hydraulic status, LAI, and carbon storage/ biomass remain rare (but see Atwell et al., 2007;Nadal-Sala et al., 2021;Poyatos et al., 2013).Future that examines this link across species ecosystems will be critical.

| CONCLUSION
In this study, we demonstrated that competing nitrogen-based and hydraulics-based optimization schemes, impacted by a representation of hydraulic damage, affected model simulations in a predictive manner along axes of variation in soil moisture, LAI, and atmospheric [CO 2 ], via interactive physiological adjustments that could also be interpreted theoretically.This result should be generalizable beyond EucFACE, opening avenues to reconcile alternative optimality principles and to help inform predictions of tree species' resilience to climate change.
We also showed that the Eucalyptus trees growing at EucFACE could suffer significant canopy PLC, but this did not lead to sustained hydraulic damage in our model simulations.Instead, legacy effects were accounted for through LAI dynamics that mostly prevented hydraulic damage.Leaves were grown during wet years to replenish carbon stores, whereas foliage declined in anticipation of water stress.Overall, capturing these dynamics during drought appears to be a key research frontier: models that ignore or do not accurately represent LAI responses to water availability risk predicting premature canopy dieback and/or plant mortality.
whether competing optimality principles can successfully predict multiple axes of functional change in response to change in [CO 2 ].Observations of LAI, water fluxes, leaf water potential, leaf photosynthetic capacity, and leaf nitrogen content, available under both ambient and elevated [CO 2 ] at EucFACE, provide us with an opportunity to probe competing combinations of optimality principles (see above) and to thoroughly interrogate physiological responses to weather and climate extremes (i.e., drought and heatwaves).
LENS OF COMPETING ADJUSTMENTS TO VEGETATION FUNCTION | 3 where A c (µmol m −2 s −1 ) and A j (µmol m −2 s −1 ) are the ribulose- 1,5-bisphosphate (RuBP) carboxylation (Rubisco) limited and the RuBP regeneration limited rates of carbon assimilation obtained from the photosynthesis model (Methods S3); and Σ ante is a set of average antecedent leaf-to-canopy scalers and canopy light environment, canopy leaf temperature (T l ; °C), [CO 2 ] in the leaf intercellular air spaces (C i ; Pa), and R d recorded daily under maximum sunlight.

2. 3
| Model simulations 2.3.1 | Environmental drivers and model parameters Our model was forced with in situ half-hourly meteorological data (photosynthetic photon flux density (PPFD), air temperature (T a ), vapour pressure deficit (D a ), precipitation, atmospheric pressure, and wind speed) and atmospheric [CO 2 ] (C a ) that were recorded, gap-filled (<0.8% of raw data), and aggregated to 30-min timesteps as per Yang et al. (2020) and Mu et al. (2021).Canopy phenology was prescribed by the observed leaf area index (LAI) at each ring.LAI was calculated following Yang et al. (2020) and Mu et al. (2021), by removing the woody contribution (0.8 m 2 m −2 ) from the plant area index based on Note that H leg and N opt can only indirectly influence one another, as they interact through the canopy gas exchange routine.[Color figure can be viewed at wileyonlinelibrary.com] the model simulations of Ψ l .The parameterization of r k was established from the three to four (co-)dominant trees in each ring, whereas the model evaluation considers 493 measurements of Ψ l made on 37 trees (4-7 trees per ring).

Figure 2
Figure 2 qualitatively illustrates the relative effects of differences in WUE between the default model configuration and all model other configurations, hereafter referred to as ΔWUE config .These effects are shown for each of the aCO 2 (panel a) and eCO 2 (panel b) rings, and the rings appear in order of increasing soil moisture × LAI from left to right within each panel.Negative ΔWUE config were simulated at each ring except R5, which indicated that the combinations of competing optimization schemes and hydraulic legacies mostly led to a reduction in WUE compared to the default model.However, overall, the effects were smaller for the eCO 2 than the aCO 2 rings (as shown by smaller relative horizontal span for the eCO 2 rings), implying that the eCO 2

N
opt = 28-35 days).The best model configuration combined H leg = 30 days with N opt = 7 days.PREDICTING RESILIENCE THROUGH LENS OF COMPETING ADJUSTMENTS TO VEGETATION FUNCTION | 7

Figure 5
Figure 5 shows comparisons between the GPP, E, and the canopy PLC simulated by the default and the best model configurations.The functional adjustments of the best model configuration had little effect on the fluxes over time, with c.

F
I G U R E 5 A comparison between the gross primary productivity (GPP; a and b), transpiration (E; c and d), and canopy percentage loss of hydraulic conductivity (PLC; e and f) simulated by the default and best (H leg = 30 days, N opt = 7 days) model configurations at the aCO 2 (a, c, and e) and eCO 2 (b, d, and f) rings.In each panel, the plain lines show the average simulations across rings (e.g., average of R2, R3, and R6 for the aCO 2 rings) and the faint lines show the range of simulations.Upscaled observations of tree-level transpiration are included for comparison with the model simulations, and the insets in (c) and (d) show the correspondence between the daily observations (x-axis) and their paired simulated estimates (y-axis).

|
Figure 6 further explores the role of LAI in response to water by showing how observed deviations in LAI from the average phenology impact the simulated canopy PLC (a and b), considering associated carbon benefits and risks (c and d).Higher LAI than the

|
Representing legacies from climate and weather extremesOur model-based approach suggests that the E. tereticornis trees at EucFACE were resilient to recurring drought (i.e., yearly soil moisture dry-down, multiyear drought after mid-2017) and heatwave events F I G U R E 6 The relative costs and benefits associated with deviations in leaf area index (LAI) from the average phenology (LAI) at the aCO 2 (a and c) and eCO 2 (b and d) rings, as simulated by the best model configuration (H leg = 30 days, N opt = 7 days).(a) and (b) show timeseries of the percentage point changes in canopy percentage loss of hydraulic conductivity (ΔPLC) associated with percentage changes in LAI relative to the average phenology of each ring (Δ rel LAI).The ΔPLC, that is, the PLC associated with the observed LAI minus the PLC associated with the average phenology (PLC LAI )are coloured to indicate the PLC there would have been in the absence of deviation from the average phenology (PLC LAI ).Overlaid purple shadings show each of the aCO 2 and eCO 2 ring's Δ rel LAI (i.e., there are individual shadings for R2, R3, and R6 in panel a).(c) and (d) show the cumulative portion of assimilated carbon that is allocated to leaf growth (A fcum ) over the hydrological year (July 1st to June 31st), depending on the LAI prescribed in the model simulations.Simulations driven by the observed LAI are coloured in the dark blue and those driven by LAI are coloured in pale blue.Inset plots show the cumulated differences in A fcum at the end of each hydrological year (Δ y A fcum ).Note, 2020 does not appear in these inset plots because the end of the 2019-2020 hydrological year is absent from our simulations.[Color figure can be viewed at wileyonlinelibrary.com] (i.e., temperatures >40°C and vapour pressure deficit >6 kPa every between 2013 and 2020.For instance, nonnegligible canopy PLC (c.9.5%-28.5%)was simulated the spring/summer dry-down of 2013, yet both LAI (Figure S1) and the transpiration flux (Figure 5) seemingly recovered to "peak" levels by autumn 2014.Further marked canopy PLC was simulated in mid-2014, mid-2016, and late 2017, although the exact magnitude PLC emerging from our model relates to its ability to represent the subsurface extraction of water by roots (cf.Figures 5 and PREDICTING RESILIENCE THROUGH LENS OF COMPETING ADJUSTMENTS TO VEGETATION FUNCTION | 13 by the low LAI at EucFACE (cf.predictions by the forest canopy model MAESPA; Yang Observational averages are different from the parameterizations because the parameter values derive from campaign averages (equal weights for each campaign; see Methods S6) whereas these averages apply equal weights to each date of collection.*Theseparametervalues are the same as the Vcmax t 25 0 presented in Table1.***Simulation averages only include predictions at dates for which observations are available.
. ******The average Vcmax t 25 0 from across the aCO 2 rings (i.e., R2, R3, and R6) is used for all rings.largest in spring/summer and fluctuated annually (e.g., −5 gC m −2 y −1 [−0.3%] vs. −86 gC m −2 y −1 [−11.7%] at R2 in 2014 [high LAI year] versus 2018 [low LAI year]).Even when GPP was markedly reduced by the best model configuration, the canopy gas exchange scheme kept E roughly level with the default model configuration, which limited further drops in GPP but reduced WUE.Departure by the best model configuration from the default was more pronounced at the aCO 2 rings (−72 gC m −2 y −1 GPP; −2.3 (+9.7% WUE and +6.8% WUE in the best and default configurations, respectively), because E was higher at the eCO 2 than at the aCO 2 rings (see Discussion).Our model captured Ψ l well (insets of Figure 5e,f).Marked canopy PLC (Figure 5e,f) was simulated in the second half of 2014, in mid-2016, and in late 2017 (p 99 > 45% vs. < 10% in 2015, 2018, Crous et al., 2021))odel's scaling from the leaf-to the canopy-level is inadequate when LAI is very low, or the leaf nitrogen allocation routine wrongly works to correct for the LAI reduction during drought, or both.In alternative work predicting seasonal variations in V cmax25 from an optimality principle,Jiang et al. (2020a)argued the need to account for site-specific coefficients of nitrogen extinction within the canopy (κ N ) correlated with LAI.With help from existing observations of leaf nitrogen, V cmax25 and J max25 depending on canopy height (e.g.,Crous et al., 2021)for EucFACE), future studies may evaluate whether varying κ Nreduced low-LAI simulation biases.Errors in the simulation of V cmax25 at low LAI depended on whether the amount of leaf nitrogen available for photosynthesis, N P , was initialised at the ring-level (see Table2).However, the simulated dynamics of N P aligned with observed adjustments in leaf nitrogen regardless(Figures S2 and S4).N P declined under elevated [CO 2 ], due to a simulated reduction in the soluble protein nitrogen invested in Rubisco concomitant to an increase in other soluble protein investments.This shift in the nitrogen invested into photosynthetic compounds is echoed by observational evidence for the upper canopy of E. globulus grown at the Hawkesbury Forest Experiment near EucFACE