A triple tree-ring constraint for tree growth and physiology in a global land surface model

Annually-resolved tree-ring records extending back to pre-industrial conditions have the potential to constrain the responses of global land surface models at interannual to centennial time scales. Here, we demonstrate a framework to :::::::::::: simultaneously : constrain the representation of tree growth and physiology in the ORCHIDEE global land surface model using the simulated variability of tree-ring width and carbon (∆C) and oxygen (δO) stable isotopes in six sites in boreal and temperate Europe. We exploit the ::::::: resulting : tree-ring triplet to derive integrative constraints for leaf physiology and growth 5 from well-known mechanistic relationships among the variables. The model :::::::::: ORCHIDEE : simulates ∆C (r = 0.31-0.80) and δO (r = 0.36-0.74) better than tree-ring width (r < 0.55), with an overall skill similar to that of other models. The results show a :::::::: tree-ring ::::: model :::::::::::: (MAIDENiso) :::: and :::::: another ::::::::::::: isotope-enabled :::::: global :::::::: vegetation :::::: model :::::::::: (LPX-Bern). :::: The :::::::::: comparison :::: with ::::::: tree-ring :::: data :::::: showed : that growth variability is not well represented :: in :::::::::: ORCHIDEE, and that the parameterization of leaf-level physiological responses ::::::: (stomatal ::::::: control) : to drought stress in the temperate region can be improved with ::::::::: constrained ::::: using 10 :: the :::::::::: interannual :::::::: variability ::: of tree-ring data :::: stable ::::::: isotopes. The representation of carbon storage and remobilization dynamics is critical ::::::: emerged :: as :: a :::::: critical :::::: process : to improve the realism of simulated growth variability, temporal carrying over :::::::: carryover and recovery of forest ecosystems after climate extremes. The simulated physiological response to rising CO2 over the 20th century is consistent with :::::::: Simulated :::::: forest ::::: Gross ::::::: Primary :::::::::: Productivity :::::: (GPP) :::::::: correlates ::::: with :::::::: simulated : tree-ring data in the temperate region, despite an overestimation of seasonal drought stress and stomatal control on photosynthesis. Photosynthesis 15 correlates directly with isotopic ::::: ∆C ::: and ::::: δO : variability, but the origin of correlations with :: the :::::::::: correlations :::: with :::::::: tree-ring δO is not entirely physiological. The integration of tree-ring data and land surface models as demonstrated here should ::::: guide ::::: model :::::::::::: improvements ::: and : contribute towards reducing current uncertainties in forest carbon and water cycling.


Introduction
A major challenge for the Land Surface Model (LSM) component of the Earth System models currently used to project climate change is to accurately simulate the historical and future dynamical coupling between the global biosphere and climate (Friedlingstein et al., 2014). Although LSMs are skillful at reproducing short-term (<20 years) contemporary observations of plant water and carbon cycling, their responses at longer time scales from decades to century are still highly uncertain 20 and contribute to the spread in current climate change projections (Ciais et al., 2013;Friedlingstein et al., 2014). Some of these models project that the terrestrial biosphere will continue behaving as a carbon sink of anthropogenic emissions during canopy structure in ORCHIDEE r5698 (Naudts et al., 2015;Vuichard et al., 2019), where the tree-ring functionality will also be ported (Jeong et al., 2020).

Tree-ring width
The forest management module explicitly simulates the temporal trajectory of stem growth of trees in a size-structured forest 5 stand. The size structure of the stand is prescribed by defining a starting density of trees per hectare distributed across a number of stem diameter size classes, typically following an inverted-J distribution. Initial tree height and biomass are obtained by allometric relationships with stem diameter (Bellassen et al., 2010). For each PFT, photosynthesis and NPP are computed at the stand level at half-hourly and daily steps, respectively.
At the end of the year, accumulated woody NPP increment is allocated across the different diameter classes following an 10 empirical competition rule for even-aged stands (Deleuze et al., 2004). A greater NPP share is allocated to larger dominant trees compared to smaller less vigorous trees, emulating the competition for light and resources within the stand. The absolute annual stem growth increment in the allocation rule is determined by a parameter (λ) that defines the slope between the fraction of NPP allocated to each diameter class and the mean diameter of the class. A second parameter (σ) represents a diameter threshold below which less vigorous trees receive only a fixed part of the yearly stand NPP increment. These parameters are saved 15 through the simulation and allow reconstructing the growth trajectory of each diameter class. In addition, the stem growth of the largest and smallest tree is tracked individually through the entire simulation (Bellassen et al., 2010).
Simulated tree-ring width (i.e., radial growth) for each diameter class was computed from annual increments in stem circumference as follows: (1) 20 where TRW t is the annual tree-ring width in millimetres for year t and Circ t and Circ t−1 are the stem circumferences for the current year and the previous year, respectively. The annual ring-widths simulated in this scheme account for the size-related trend in stem growth increment and PFT-specific tree allometry, providing a mean representation of radial growth for trees across a range of sizes that can be meaningfully compared with measured tree-ring width data. However, in this formulation tree-ring width variability still depends only on direct GPP allocation (i.e., it is a carbon source driven process) and does 25 not consider the dynamics of carbon storage or wood formation processes that can account for a large fraction of observed interannual tree-ring width variation (Misson, 2004).

Tree-ring carbon isotopes
To distinguish between variations in the δ 13 C of atmospheric CO 2 and the effects of plant metabolic processes, the 13 C enrichment of plant organic material is usually expressed in terms of the carbon isotope discrimination (∆ 13 C). Differences between 30 the 13 C enrichment of atmospheric CO 2 and plant material are attributed to discriminatory processes during photosynthesis.
Photosynthetic carbon isotope discrimination was estimated at half-hourly time steps using the simple formulation of Farquhar et al. (1982) for C3 plants: where a (4.4‰) is the kinetic discrimination associated with diffusion between free air and the stomatal cavity, b (27‰) is the 5 fractionation during CO 2 fixation by the Rubisco enzyme, c i is the leaf internal CO 2 concentration simulated by ORCHIDEE and c a is the atmospheric CO 2 concentration prescribed from measurements. Annual tree-ring carbon discrimination (∆ 13 C R ) was then calculated as the mean of the half-hourly discrimination values in the above-ground biomass weighted by Gross Primary Productivity (GPP) between the start (SOS) and the end (EOS) of the growing season: 10 This formulation for carbon discrimination is commonly used as a simple approximation for discrimination derived from measured δ 13 C in tree-ring cellulose (Francey and Farquhar, 1982). For simplicity, it assumes that further post-photosynthetic fractionation during photo-and dark-respiration and carbohydrate remobilization and storage is negligible. Although these processes normally have an important impact on whole-ring cellulose isotopic composition (Gessler et al., 2009;Werner et al., 2012), at least the impact of carbon remobilization is minimal in the latewood component of tree rings (Helle and Schleser, 15 2004). The c i term in Eq.
(2) integrates the gas exchange dynamics (i.e., stomatal conductance and photosynthesis) simulated by the model as a complex function of micrometeorological variability, seasonal water stress, and the long-term warming and increase in c a . Hence, despite its simplicity, simulated carbon discrimination using Eq. (2) has been shown to be valuable to constrain the integrated environmental response of land surface models using carbon discrimination derived from tree-ring δ 13 C measurements (e.g., Bodin et al., 2013;Churakova et al., 2016;Keel et al., 2016).

Tree-ring oxygen isotopes
The δ 18 O signature in tree-ring cellulose reflects primarily the isotopic composition of soil water and the evaporative enrichment of leaf water due to transpiration, but other mixing and biochemical fractionation processes during water transport along the soil-plant-atmosphere continuum also contribute to the final δ 18 O signature (McCarroll and Loader, 2004;Gessler et al., 2014). The δ 18 O isotopic composition of the source water used by plants may originate from rainfall or groundwater, while the 25 δ 18 O isotopic enrichment of leaf water depends strongly on vapour pressure deficit (i.e., the ratio of the vapour pressure in the atmosphere and the intercellular spaces of the leaves) or relative humidity (Farquhar et al., 1998;Scheidegger et al., 2000).
The δ 18 O fractionation and mixing processes in all water pools and fluxes along the soil-plant-atmosphere continuum (Risi et al., 2016) is represented following a similar formulation than in other isotope-enabled global land surface models (Aleinov and Schmidt, 2006;Haese et al., 2012). The isotopic compositions of precipitation and near-surface water vapour have to be 30 prescribed monthly when running the model stand-alone or are simulated by the LMDZ general circulation model in coupled simulations. Precipitation reaching the soil surface or intercepting the canopy fractionates during evaporation according to the 6 https://doi.org/10.5194/bg-2020-446 Preprint. Discussion started: 10 December 2020 c Author(s) 2020. CC BY 4.0 License.
Craig and Gordon equation (Craig and Gordon, 1965), which generically describes the preferential evaporation of the lighter isotope of a free water body at steady state.
The resulting isotopic composition of soil water and mixing are parameterized as a vertical profile to overcome the limitation of depth resolution in the two-layer representation of soils in the model. No isotopic fractionation is assumed to occur during 5 absorption of soil water by roots and thus the isotopic signature of xylem water is the same as that of soil water. The isotopic composition of leaf water at the evaporation sites (δ 18 O e ) is diagnosed by inverting the Craig and Gordon (1965) equation: where δ 18 O sw is the isotopic composition of soil water taken up by the roots integrating older soil water and recent precipitation, δ 18 O v is the isotopic composition of atmospheric water vapour, α is the equilibrium fractionation due to the phase change 10 from liquid water to vapour, α k is the kinetic fractionation due to diffusion of vapour into unsaturated air, and h is the relative humidity normalized to surface temperature.
Isotopic enrichment of leaf water in the mesophyll (δ 18 O lw ) results from mixing between isotopically enriched leaf water at the evaporative site and depleted xylem water (δ 18 O xw ) through the so-called Peclet effect:

15
where f = (1 − e −P )/P is a coefficient that decreases as the Peclet effect increases, with P = (E × L)/(W × D) as the Peclet parameter. E is the transpiration rate per leaf area, L is the effective diffusion length and W is the leaf water content per leaf volume. L was set to 8 mm and W was assumed to be 103 kg m −3 (Risi et al., 2016).
Half-hourly tree-ring cellulose isotopic composition (δ 18 O cell ) is calculated from the isotopic composition of leaf water (δ 18 O lw ) and xylem water (δ 18 O xw ) following the formulation of Anderson et al. (2002): Where f 0 is the fraction of leaf water exchanged with xylem water prior to cellulose synthesis, which reduces the imprint of leaf water on cellulose. For δ 18 O this exchange is estimated to be 0.42 based on the best fit relationship under controlled experiments (Roden et al., 2000). The parameter is the biochemical fractionation factor during cellulose formation associated with water carbonyl interactions and is estimated to be 27‰ (DeNiro and Epstein, 1979;Stenberg and DeNiro, 1983). 25 An estimate of growing season tree-ring isotopic composition (δ 18 O R ) is obtained by weighting δ 18 O cell by GPP as done for ∆ 13 C R in Eq. (3): A model evaluation across a network of 10 sites in Europe and North America shows that this representation is able to reproduce the main features of the seasonal and vertical variations in soil water isotope content, as well as the seasonal isotopic 7 https://doi.org/10.5194/bg-2020-446 Preprint. Discussion started: 10 December 2020 c Author(s) 2020. CC BY 4.0 License. signals in stem and leaf water (Risi et al., 2016). The isotopic variability simulated by ORCHIDEE has been used to interpret local climatic signals in boreal tree-ring δ 18 O records (Shi et al., 2011;Churakova et al., 2016) and to investigate regional and global isotopic signatures of continental recycling (Risi et al., 2013).  2004;Danis et al., 2006). In all sites, except Annecy, the available ring-width and δ 13 C and δ 18 O chronologies cover the common period 1960-2000.
The tree-ring width measurement series available for each site were standardized by dividing each series by its mean ring width (Cook et al., 1990). The resulting series of tree-ring width indices were averaged together to produce a mean site 5 chronology composed by seven to twenty-eight trees (Table 1). This simple standardization method allows the computation of average tree-ring chronologies without the average being dominated by the faster growing trees with wider rings. Depending on the site, 4 to 9 trees were selected to develop carbon ( δ 13 C R ) and oxygen (δ 18 O R ) stable isotope chronologies for each site by pooling rings from the selected trees (Table 1) using standard methods for cellulose extraction and measurement of the isotopic ratios (McCarroll and Loader, 2004;Daux et al., 2018). For oaks, earlywood and latewood sections were separated 10 but only latewood was processed. For pine, the carbon and oxygen isotopic compositions were measured for the whole ring.
Tree-ring carbon discrimination was computed by subtracting the stable carbon isotope composition of the atmosphere (δ 13 C a ) from the measurements using the expression of Farquhar et al. (1982): For Fontainebleau, the intrinsic Water Use Efficiency (iWUE) was calculated from c i and c a for the period 1901-2000 15 following the expression (Farquhar et al., 1982):

Simulations
The model was run at each tree-ring site over the period 1901-2000 using as meteorological forcing the nearest 0.5 o grid of 25 the 6-hourly CRU-NCEP dataset (Viovy, 2014). This gridded forcing dataset is a statistical merging of the monthly CRU TS station-based dataset of the Climate Research Unit (New et al., 2000) with the atmospheric reanalysis from the National Center for Environmental Prediction (NCEP). The corresponding soil type and vegetation PFT were prescribed for the sites (Table 1).
Monthly initial tree density of 1000 trees per hectare to approach current forest age and density following tree mortality over time due to self-thinning.
For a comparison with ORCHIDEE, a simulation of the MAIDENiso tree-ring model for the Fontainebleau forest was obtained from an earlier study (Danis et al., 2012). The simulation was produced using a site-specific calibration and includes

Model-data comparison
The ability of ORCHIDEE to simulate the interannual variability of the three tree-ring parameters over the course of the 10 20th century was first evaluated in the Fontainebleau forest (Table 1). Fontainebleau is a well-studied tree-ring site in France (Michelot et al., 2011(Michelot et al., , 2012Daux et al., 2018) and it has been used to evaluate the MAIDENiso tree-ring model (

15
A simulated tree-ring width chronology was derived for each site by dividing the simulated tree-ring width series of the largest model tree by its mean. Since growth allocation in the model increases almost linearly with stem size, the absolute annual ring-width vary among stem size-classes but its interannual variability remains similar across all size classes, thus the choice of size class does not affect the standardized variability. The standardization removes the effect of stem size-class and conserves the interannual and longer variability but does not remove the juvenile effect in tree-ring width. However, the juvenile 20 trend in simulated and observed ring-width does not affect the evaluation period  because at this time trees were already mature canopy individuals with their radial growth fluctuating around the mean.
Simulated δ 18 O and ∆ 13 C tree-ring chronologies were obtained by averaging the simulated half-hourly isotopic variability between May and August. Using just the summer season (Jun-August) improves results for ∆ 13 C but degrades substantially results for δ 18 O, thus May was included as a compromise in order to use a common season for the isotopes and ensure Since the focus of our evaluation is on the interannual variability and not on the absolute values, correlation and the normalised standard deviation (i.e., the standard deviation of simulated tree-ring parameters divided by the standard deviation of the observations) were used to quantitatively evaluate the skill of the models to simulate the variability in tree-ring width and 30 stable isotopes. Differences in the temporal persistence (i.e., carrying over effect) in the observed and simulated tree-ring parameters were evaluated using the first-order autocorrelation of the time series. The climate response of simulated ring width, The integrated growth-isotope responses simulated by ORCHIDEE and MAIDENiso models in Fontainebleau were qualitatively compared with observations using a bivariate response surface between tree-ring width variability and the joint variability of the two isotopes. A smoothed response surface was fitted using a data-adaptive bivariate generalized additive model (GAM) implemented with the mgcv package (Wood, 2017) in the R environment (R Core Team, 2020). This visualization device 5 extends the dual isotope conceptual model of Scheidegger et al. (2000) to illustrate the mechanistic information content of the triple tree-ring constraint for models introduced in this study. It neatly reveals the complex association of tree growth with gas exchange inferred from stable isotopes in both observations and models.
To disentangle the relative importance of source water (δ 18 O sw ) and leaf water enrichment above source water (∆ 18 O lw = δ 18 O lw − δ 18 O sw ) between May and August in determining the variability of δ 18 O R simulated by ORCHIDEE we used the 10 Lindeman-Merenda-Gold (LMG) method (Grömping et al., 2006). It allows to quantify the contribution of different correlated regressors (here δ 18 O sw and ∆ 18 O lw ) to the total r 2 of a multiple linear regression model.

Model evaluation in Fontainebleau
3.1.1 Tree-ring width and isotopic variability 15 ORCHIDEE shows a significant skill in simulating the interannual and multidecadal variability of oak tree-ring width (r=0.59, p<0.01) and latewood ∆ 13 C R (r=0.41, p<0.01) and δ 18 O R (r=0.49, p<0.01) over the 20th century in Fontainebleau (Fig.1ac). The magnitude of the interannual variability of δ 18 O R is well simulated (NSD=1.04) but that of tree-ring width is overestimated by 37% (NSD=1.37), while that of ∆ 13 C R is underestimated by about a similar magnitude (NSD=0.55). ∆ 13 C R is systematically overestimated since 1980, when the observations show a decrease by about 1‰ (Fig. 1b). Overall, the model 20 simulates 35% of observed tree-ring width variability and 17-24% of latewood isotopic variability over the 20th century.
The first-order autocorrelation or carrying over in observed tree-ring width is significant (r lag1 =0.54, p<0.001) and its magnitude indicates that ring width in the previous year explains up to 30% of current year tree-ring width variability. In contrast, simulated tree-ring width has no first-order autocorrelation (r lag1 =0.02, p>0.1) because ORCHIDEE does not account for the carrying over in radial growth associated with carbon storage and remobilization. As a result, in years with extreme summer 25 drought conditions like in 1921 and 1976 the model does not simulate any stem growth because photosynthesis is strongly suppressed (Fig. 1a). The lack of modulation of simulated radial growth by carbon storage dynamics makes the recovery after these extremes too fast compared with the observations. The first-order autocorrelation in the isotopic observations is significant only for ∆ 13 C R (r lag1 =0.37, p<0.001) while in the simulations it is marginally significant for both ∆ 13 C R (r lag1 =0.18, p<0.1) and δ 18 O R (r lag1 =0.27, p<0.01).

30
The skill of ORCHIDEE compares well with that of the specialized MAIDENiso tree-ring model (Fig. 1a-c), which was specifically calibrated for the site over the period 1953-2000. Over this period, MAIDENiso is able to simulate between 30% and 46% of the total variability of the observations of tree-ring width and isotopes, which compares to the 20% to 44% of the total observed variance simulated by ORCHIDEE over the same period with standard parameterization. MAIDENiso is considerably better than ORCHIDEE at simulating tree-ring width (r=0.68 vs r=0.51) and ∆ 13 C R (r=0.58 vs r=0.45) variability, but despite simulating well the amplitude of ∆ 13 C R and δ 18 O R (NSD=0.91-1.16) it substantially underestimates the amplitude of tree-ring width (NSD=0.68). Unlike ORCHIDEE, it is able to simulate a significant first-order autocorrelation 5 in tree-ring width (r lag1 =0.45, p<0.01) with a magnitude similar to that of the observations (r lag1 =0.50, p<0.01). This carrying over effect accounts for up to 25% of current year tree-ring width variability in the observations  and is thus an important component of the growth variability captured by its carbon remobilization dynamics.
MAIDENiso is able to simulate the observed decrease in ∆ 13 C R since 1980 better than ORCHIDEE (Fig. 1b)

Observed and simulated relationships among tree-ring variables
So far the performance statistics used above only describe the unidimensional skill of the models and do not evaluate their abil- 15 ity to simulate the joint relationships that exist among tree-ring width and ∆ 13 C R and δ 18 O R . The strength of the correlations between observed tree-ring width and isotopic variability in Fontainebleau indicates that oak radial growth has 10% and 24% of common variance with latewood δ 18 O R and ∆ 13 C R , respectively (Fig. 1d). In turn, the isotopes have 16% of common variability. The interpretation of these three-way relationships in the observations can be aided by visualizing the bivariate surface response of tree-ring width as a function of the dual latewood isotope variability (Fig. 1d). The resulting surface provides, at a 20 glance, insights on causal and non-causal (indirect) relationships between environmental variability (temperature and drought stress) and stomatal responses and growth.
The bivariate surface response shows that in this site observed ring-width is positively and linearly related to latewood ∆ 13 C R (r=0.48, p<0.01), whereas the relationship with δ 18 O R is negative (r=-0.31, p<0.05) and non-linear (Fig. 1d). This indicates that narrow rings during dry years like in 1976 are mechanistically linked to reduced latewood ∆ 13 C R because trees 25 tend to close their stomata for longer during the summer to reduce water loss at the expense of reducing their photosynthesis.
The opposite occurs during moist years, when growth, stomatal conductance and photosynthesis are high (Fig. 1d).
The apparent relationship between growth and δ 18 O R suggests interactions between the factors (air temperature, relative humidity and soil moisture) driving the growth process and δ 18 O enrichment of source and leaf water (see axes in Fig. 1d). For instance, years with high temperature tend to be also dry (e.g., 1976, 1990 and 1996) and associated with reduced atmospheric 30 humidity and soil moisture, resulting in reduced growth and increased δ 18 O enrichment of source water, vapour and leaf water due to higher evaporation rates and stomatal closure. The negative linear relationship between the two latewood isotopic variables (r=-0.40, p<0.01), apparent in the scatter of the observation points across the surface (Fig. 1d), provides evidence for a significant degree of stomatal control (g s ) of trees to avoid dehydration under warmer and drier conditions.
The two models simulate qualitatively different growth-isotope surface responses (Fig. 1e-f). ORCHIDEE simulates a surface response slightly more consistent with the observations than MAIDENiso, but none of the two simulations captures the observed non-linear relationship between tree-ring width and δ 18 O R that might arise from source water or isotopic enrichment of leaf water. Although ORCHIDEE simulates a significant relationship of tree-ring width with δ 18 O R (r=0.54, p<0.001) and 5 ∆ 13 C R (r=0.87, p<0.001), the strength of the relationships is overestimated (Fig. 1e). Subtracting the isotopic variability of to isolate the leaf signal does not change much the surface response of ORCHIDEE but increases substantially the correlation between the residual isotopic variability (∆ 18 O R ) and ∆ 13 C R and GPP-driven ring-width ( Fig. 2a-b).
MAIDENiso does not simulate any significant relationship between tree-ring width and δ 18 O R (r=0.16, p>0.1), but it cap-10 tures well the magnitude of the observed correlation (r=0.30, p<0.01) between ring-width and ∆ 13 C R (Fig. 1f). Both models largely overestimate the common variability between the isotopes, as is evident from the rather even spread of the data points across the diagonal of the surfaces (Fig. 1e-f). MAIDENiso and ORCHIDEE simulate 62% and 48% of common isotopic variability, respectively. This is three times above the observed 16%, indicating that the models overestimate the stomatal responses to summer drought stress.

Climate response
The magnitude of the seasonal climatic responses of the observed tree-ring width, ∆ 13 C R and δ 18 O R in Fontainebleau is well captured by the models, but there are some important differences in the timing and duration of the period of significant response (Fig. 3). The seasonal correlation patterns with precipitation and VPD indicate that tree-ring width and ∆ 13 C R in ORCHIDEE are too sensitive to moisture variability over the growing season compared with the observations (Fig. 3). In 20 contrast, ring-width in MAIDENiso has too little sensitivity to precipitation variability. Nevertheless, MAIDENiso captures well the observed climatic response of ∆ 13 C R to summer VPD, while the effect of summer VPD on δ 18 O R is better captured by ORCHIDEE and its isotopic forcing.

20th century change in water use efficiency
The ∆ 13 C R data show that in Fontainebleau the observed intrinsic water use efficiency (iWUE) of oak measured as the change 25 between 1901-1910 and 1990-2000 has increased by 25.5% following the gas-exchange scenario of constant c i /c a (Fig. 4).
Most of the change followed the steady increase in atmospheric CO 2 concentration since 1960, though the increasing trend in iWUE stalled around 1980 while atmospheric CO 2 continued increasing. ORCHIDEE produces a slightly lower increase in iWUE over the 20th century (21.2%), but also follows the gas exchange response of constant c i /c a . The model underestimates iWUE since around 1980 and in contrast to the observations simulates a steady increase since 1960. This model-data mismatch 30 is linked to the overestimation of ∆ 13 C R in ORCHIDEE during this recent period as described earlier, an issue that does not affect MAIDENiso (Fig. 1b).

Simulated relationship between productivity and isotopic variability
Simulated isotopic variability in Fontainebleau is significantly correlated with growing season (May-August) GPP (Fig. 4b).
The magnitude of the correlations indicate that simulated δ 18 O R (r=-0.62, p<0.001) and ∆ 13 C R (r=0.84, p<0.001) explain 38% and 71% of GPP variability over the 20th century, respectively. Since stem growth in ORCHIDEE depends directly on the 5 allocation of GPP, simulated tree-ring width variability is by definition strongly correlated with growing season GPP (r=0.90, p<0.001) as it can be seen in Fig. 4b.

Model performance across sites
The performance of ORCHIDEE for tree-ring width, ∆ 13 C R and δ 18 O R varies substantially across sites (Fig. 5), with no clear pattern along the climate gradient from Finland to France or between species for any parameter. However, it is clear that the 10 isotopic variability is better simulated than tree-ring width and also that δ 18 O R is the tree-ring variable best simulated by the model. Tree-ring width is well simulated (25-30% of the observed variability) at only two out of six sites. It is not clear whether the fact that the best simulations are for the two southernmost sites (Fontainebleau and Annnecy) is a coincidence or suggest that the model processes and/or parameters are biased in favour of deciduous temperate forests. In the remaining sites the simulations account for less than 10% (r<0.32) of the observed variability. The lower ability of ORCHIDEE to simulate 15 tree-ring width is to a large extent due to the present inability of the model to simulate the significant carrying over effect of growth (autocorrelation) evident in the observations (Fig. 7a).
Although ORCHIDEE is able to simulate about 10-64% (r=0.31-0.80) of the observed ∆ 13 C R variability, it tends to underestimate its amplitude by 30-60%, particularly in the northernmost sites of pine (Kessi and Sivakkovaara) where the isotopic ratios were measured over the whole ring (Fig. 5). With the exception of the northernmost site, the amplitude of δ 18 O R variabil-20 ity is simulated within ±20% and the simulations account for 13-55% (r=0.36-0.74) of the variance of the observations. The primary driver of simulated δ 18 O R variability in five out of the six sites is the isotopic composition of source water (δ 18 O sw ).

25
ORCHIDEE overestimates the correlation between δ 18 O R and ∆ 13 C R in the temperate sites in France, but it simulates very well the magnitude of the isotopic coupling observed in the boreal oak and pine sites in Finland (Fig. 6b). This means that the simulated stomatal control and responses to atmospheric humidity are overestimated in the temperate deciduous PFT, as is apparent in a stronger correlation between simulated isotopic variability and VPD in these sites compared with the observations (Fig. 6c-d). The significant relationship found between simulated δ 18 O R and GPP in Fontainebleau is also observed across all the other sites (Fig. 8a). Its strength (r=-0.47 to -0.73) appears to increase from the boreal to the temperate region, accounting for up 5 to 22-53% of GPP variability. This empirical relationship is driven by a synergistic effect of source water and leaf water enrichment above source water because both variables correlate negatively with GPP (Fig. 8b). Leaf water enrichment tends to correlate with GPP more strongly than source water, though in the northernmost pine site only source water is correlated with GPP. ∆ 13 C R is also significantly correlated with GPP in most of the sites, but correlations are insignificant or change sign in the colder northernmost sites (Fig. 8a). Tree-ring width is by definition highly correlated with GPP (Fig. 8a).

Integrating tree-ring width and carbon and oxygen isotopes for insights on growth and gas exchange
The enhanced information content obtained by combining multiple and complementary tree-ring variables has long been recognized in multi-proxy dendroclimatology (McCarroll et al., 2003;Loader et al., 2008;Hilasvuori et al., 2009;Daux et al., 2011;Schollaen et al., 2013;Loader et al., 2015) and dendroecology (Guerrieri et al., 2009;Savard, 2010;Leonelli et al., 2012;15 Shestakova and Martínez-Sancho, 2020), but its potential remained largely untapped in ecological modelling. Tree-ring widths provide a historical record of annual aboveground biomass increment (Clark et al., 2001;Bouriaud et al., 2005;Babst et al., 2018;Cernusak and English, 2015;Foster et al., 2016;Dye et al., 2016;Evans et al., 2017;Shestakova et al., 2019). The stable carbon isotope ratio of plant material, usually reported as δ 13 C, is related to the ratio of intercellular (c i ) and atmospheric CO 2 concentration (Farquhar et al., 1982). As this ratio varies according to changes in stomatal aperture (g s ; supply of CO 2 ) and 20 assimilation rate (A; CO 2 demand), it integrates the physiological response of plants to environmental changes such as drought stress and increasing atmospheric CO 2 concentration (McCarroll and Loader, 2004;Gessler et al., 2014).
Changes in c i can result from changes in either A or g s , thus the interpretation of changes in carbon isotope alone is difficult. The imprint of leaf evaporative enrichment on δ 18 O of tree-ring cellulose is not affected by A and like the carbon isotope depends on g s (McCarroll and Loader, 2004;Gessler et al., 2014). Hence, the leaf enrichment signal of tree-ring δ 18 O 25 (∆ 18 O lw ) is an independent proxy for g s and can be used to disentangle the physiological drivers (A or g s ) of variations in carbon isotope fractionation and growth in the same tree ring (Saurer et al., 1997;Scheidegger et al., 2000;Barnard et al., 2012).
The concurrent variability of tree-ring width, ∆ 13 C R and δ 18 O R is driven by coordinated ecophysiological responses to external environmental factors that, depending on species and location, can imprint a varying degree of covariability among 30 them as seen in the tree-ring triplet (Fig. 1d). In a temperate location such as Fontainebleau, it is expected that dry summer conditions associate with narrow tree rings, decreased ∆ 13 C R and increased oxygen isotope enrichment in the leaf water (thus higher δ 18 O R ) because trees close their stomata to reduce g s and avoid dehydration during drought stress at the expense of photosynthesis (Saurer et al., 1997;Scheidegger et al., 2000;Barnard et al., 2012). Consistent with this expectation, our bivariate response surface shows that years with narrow rings like 1976, 1990 and 1996 are associated with low ∆ 13 C R and high δ 18 O R in latewood (Fig. 1d). The opposite occurs during moist years with wider rings, when conditions for stomatal aperture and photosynthesis are optimal. This demonstrates that projecting tree-ring width variability into the isotopic space allows unravelling the underlying mechanistic relationships between tree growth and gas exchange to obtain an integrated 5 picture of how trees respond to drought stress, an eventually also to temperature anomalies if δ 18 O R can be interpreted as a temperature proxy (Fig. 1d-f).
The often dominant effect of δ 18 O sw on δ 18 O R variability (Fig. 6a) dilutes the leaf enrichment signal and reduces the correlation between δ 18 O R and ∆ 13 C R (Fig. 2a). Studies using the dual isotope conceptual model to interpret physiological signals typically remove or minimize its effect (Scheidegger et al., 2000;Barnard et al., 2012;Roden and Siegwolf, 2012). In mod-10 elling studies, this signal can be quantified and used to evaluate the quality of the isotopic forcings (precipitation and vapour) of modelled δ 18 O R and attribute data-model mismatches. In addition, δ 18 O sw carries a well-known paleo-temperature signal (McCarroll and Loader, 2004) that can be used to disentangle the climatic forcing of past growth anomalies and physiological responses inferred from ∆ 13 C R well beyond the duration of the meteorological records.
The three-way correlations among tree-ring width and ∆ 13 C R and δ 18 O R (or ∆ 18 O R ) variability provide a measure of the 15 interannual coupling between growth and gas exchange, revealing direct and indirect associations driven by common or correlated environmental forcings (Fig. 1d). The strength of the direct causal relationships expected from mechanistic understanding, like the correlation between isotopes (Scheidegger et al., 2000) and between photosynthesis or biomass increment and carbon isotopes (Belmecheri et al., 2014;Battipaglia et al., 2013;Lévesque et al., 2014;Fernandes et al., 2016), is a simple benchmark to constrain the physiological parameterization of models as we discuss below.

Critical processes for the concurrent simulation of multiple tree-ring variables in global land surface models
Our evaluation in Fontainebleau and other five sites across a boreal-to-temperate climate gradient in Europe showed that ORCHIDEE simulates better the interannual variability of tree-ring stable isotopes than ring-width (30-64% and <30% of the observed variability, respectively), with a general performance for the stable isotopes similar to MAIDENiso and LPX-Bern models ( Figs.1 and 5). The lower skill for tree-ring width results from the inability of ORCHIDEE to simulate the significant 25 temporal autocorrelation observed in ring-width variability (Fig. 7a). This temporal carrying over effect is common in tree-ring width (Cailleret et al., 2018;Breitenmoser et al., 2014) and results from carbon remobilization from previous years (Kagawa et al., 2006) and to some extent from cambial dynamics (Vaganov et al., 2011). It varies considerably with species and location but typically accounts for 20-25% of current year ring-width variability (Breitenmoser et al., 2014).
The ORCHIDEE version used in this study (r898) lacks the representation of the influence of carbon storage on simulated 30 tree-ring width and instead growth reflects only current year GPP variability (see Figs. 3b and 6a). This simplified representation of growth is currently the major limitation of the model to simulate the tree-ring triplet (Fig. 5). Nevertheless, in Fontainebleau ORCHIDEE still simulates 26% of the observed interannual variability but overestimates its amplitude and the speed of recovery from drought extremes like 1976 (Fig. 1a). Legacy effects of reduced growth following drought events can persist for 1 to 4 years in drought-sensitive ecosystems (Anderegg et al., 2015;Cailleret et al., 2018). Such prolonged legacy effects are typically not simulated by global carbon cycle models because most of them, like ORCHIDEE, lack the representation of the significant dependency of interannual tree growth on carbon remobilization from storage in their carbon allocation schemes (Anderegg et al., 2015).
Unlike global models, MAIDENiso explicitly represents the autocorrelation in tree-ring width in its carbon allocation scheme 5 (Misson, 2004;Danis et al., 2012) and as a result it was able to capture the enduring effect of the extreme drought of 1976 in Fontainebleau even when tree growth is represented through GPP allocation as in ORCHIDEE (Fig. 1a). This result demonstrates that a simple approach to represent the dependency of tree growth on carbon remobilization might produce a substantial improvement in ORCHIDEE and allow a direct tree-ring constraint for the simulation of ecosystem recovery from climate extremes. Implementing wood formation dynamics would be a fully process-based approach to represent tree-ring growth (Fritts 10 et al., 1999;Friend et al., 2019), but in a global model it implies an important tradeoff with computing time.
ORCHIDEE was able to simulate 30-64% of the observed ∆ 13 C R variability along the climate gradient but the amplitude of the interannual variations was underestimated by 30-60% (Fig. 5) and the observed decrease in ∆ 13 C R since 1980 in Fontainebleau was not captured (Fig. 1b). An earlier study with ORCHIDEE also found a similar underestimation of the interannual variability of ∆ 13 C R for larch in northeastern Yakutia, where the model simulated 26% (r = 0.51) of the observed 15 variability (Churakova et al., 2016). The low amplitude of the simulated variability can be related to missing fractionation and mixing processes and the parameterization of soil and processes that affect c i such as photosynthesis and moisture responses.
A simple parameter sensitivity test in Fontainebleau (not shown) indicated that the amplitude of simulated ∆ 13 C R is very sensitive to soil depth and photosynthetic capacity (Vcmax). This suggests that simulated ∆ 13 C R and the associated drought and stomatal responses can be better parameterized using tree-ring data. Nevertheless, besides an improved paramaterization, 20 using a more complete formulation for carbon discrimination combined with a scheme of carbohydrate mixing (Hemming et al., 2001;Ogée et al., 2009;Danis et al., 2012) should contribute to improve the simulation of ∆ 13 C R and stomatal responses in ORCHIDEE.
The long-term evaluation in Fontainebleau shows that ORCHIDEE simulates no overall change in 20th century ∆ 13 C R (Fig   1b), implying that the simulated c i /c a ratio remained roughly constant as atmospheric CO 2 concentrations increased by about 25 25% (Saurer et al., 2004). Under this type of gas-exchange response, Eq. (9) shows that iWUE should increase proportionally to the relative increase in atmospheric CO 2 (Saurer et al., 2004). Indeed, the simulated centennial increase in iWUE is 21.2% with respect to the 1901-1910 period (Fig. 4a). This is still comparable with the observed increase in iWUE of 25.5%, indicating that the decrease of about 1‰ in ∆ 13 C R (decrease in c i ) since 1980 (Fig 1b) was not sufficient to shift the constant c i /c a type of response of oak to rising CO 2 in Fontainebleau (Fig. 4a). Constant c i /c a ratio over the 20th century is the most common 30 physiological response reported for trees in Europe Frank et al., 2015) and has also been correctly simulated by the LPX-Bern vegetation model Keller et al., 2017), which uses the same formulation for ∆ 13 C R (Eq. 2) than ORCHIDEE. The implication of this stomatal response for land-atmosphere interactions is a centennial-scale reduction in transpiration as trees under present environmental conditions use less water for the production of the same amount of biomass compared with earlier decades of the 20th century.
Consistent with earlier palaeoclimatic and modelling studies (Roden et al., 2000;Robertson et al., 2001;Treydte et al., 2014;Danis et al., 2006Danis et al., , 2012Keel et al., 2016), most of the simulated variability of δ 18 O R (46-85%) in ORCHIDEE is driven by δ 18 O sw (Fig. 6a). Leaf water enrichment (∆ 18 O lw ) accounted for 14-53% of δ 18 O R variability, without showing any latitudinal 5 pattern in its contribution (Fig. 6a). Nevertheless, the coolest northernmost pine site had the lowest contribution of ∆ 18 O lw to δ 18 O R variability as it would be expected from lower transpiration rates in a cool environment .
The general pattern of relative contributions of source and leaf water enrichment highlights the strong dependence and sensitivity of δ 18 O R to the choice of the precipitation isotopic forcing in modelling studies. Differences in the isotopic drivers might account for most of the observed differences among models (Fig. 5). The isotopic forcings for LPX-Bern (Haese et al., 10 2013;Keel et al., 2016) and ORCHIDEE (Risi et al., 2010) were produced by isotope-enabled global circulation models driven or nudged by reanalysis products, whereas for MAIDENiso simple regressions with temperature and precipitation were used to produce daily precipitation and water vapour δ 18 O forcings (Danis et al., 2012).
The LPX-Bern simulations of δ 18 O R consistently compared better with the observations than the simulations of ORCHIDEE and MAIDENiso, both in terms of correlations and amplitude of the variability (Fig. 5). Notably, LPX-Bern was able to 15 simulate 74% of the observed δ 18 O R variability in the grid-box corresponding to Fontainebleau, which is considerably higher than the variance accounted for by ORCHIDEE (44%) and MAIDENiso (29%) in the site. This pattern suggests that the higher performance of LPX-Bern might be related to a better isotopic forcing.
Overall, a scheme of carbon storage and remobilization dynamics (e.g., Misson, 2004) or wood formation (Fritts et al., 1999;Friend et al., 2019) should be implemented in ORCHIDEE as part of the ongoing developments to produce novel observational 20 benchmarks from tree-ring width data (Jeong et al., 2020). Such development will allow capturing the significant autocorrelation in tree-ring width (Fig. 7a) and improving the simulation of ∆ 13 C R variability, resulting in a better representation of the impacts of climate extremes on forest ecosystems. Since the choice of isotopic forcing has a large impact on the simulation of δ 18 O R , more than one forcing should be used to evaluate the magnitude of the uncertainty of source water signals in simulated

Constraining model processes with the growth-isotope tree-ring triplet
The novel simulation of the growth-isotope tree-ring triplet (ring width, ∆ 13 C R and δ 18 O R ) in a global land surface model enabled us to use known mechanistic relationships between isotopes (Scheidegger et al., 2000) and between growth and carbon isotopes (Francey and Farquhar, 1982;Cernusak and English, 2015;Shestakova et al., 2019) to benchmark the physiological responses of the model beyond the traditional use of the interannual variability or trends (Panek and Waring, 1997;Danis et al., 30 2012;Bodin et al., 2013;Churakova et al., 2016;Keel et al., 2016;Keller et al., 2017;Ulrich et al., 2019). A much stronger negative correlation between simulated ∆ 13 C R and δ 18 O R than in the observations suggests that ORCHIDEE might overestimate the limitation of photosynthetic carbon assimilation by stomatal control in the temperate region (Rennes, Fontainebleau and Annecy; Fig. 6b).
Removing the δ 18 O sw effect from simulated δ 18 O R to highlight the effect of leaf water enrichment in Fontainebleau further increases the strength of the isotopic coupling (Fig. 2). It implies that growth and gas exchange for the temperate Broadleaf Summergreen PFT of the model are likely too sensitive to atmospheric evaporative demand and drought stress. This is consistent with correlations between VPD and simulated isotopic variability being stronger than in the observations in these sites ( Fig.   5 6c-d), and with an overestimated sensitivity of simulated tree-ring width to precipitation and VPD in Fontainebleau (Fig. 3). In the boreal region (Kessi, Sivakkovvara and Bromarv), ORCHIDEE simulated a level of correlation between isotopes (Fig. 6b) and sensitivity to VPD (Fig. 6c) in good agreement with the observations, indicating a better gas-exchange parameterization for the boreal Broadleaf Summergreen and Needleleaf Evergreen PFTs. These results illustrate how tree-ring data can be used to evaluate and identify critical processes in the parameterization of a global land surface model. The parameterization of the 10 temperate PFT can be improved by assimilating tree-ring data together with other ecosystem observations at shorter time scales using a data assimilation technique (Peylin et al., 2016;Thum et al., 2017).
We did not remove the effect of δ 18 O sw on δ 18 O R in the observations because the actual isotopic composition of precipitation is unknown. Thus, we caveat that the interpretation of the direct correlations between the two isotopes ( Fig. 1) might violate the basic assumption of a dominant leaf-water signal in the dual-isotope approach (Roden and Siegwolf, 2012). Nevertheless, 15 removing simulated δ 18 O sw from δ 18 O R reduces the correlation between the isotopes in the observations for the temperate sites and further increases the correlation in the simulations (not shown), reinforcing the finding that ORCHIDEE overestimates stomatal control in the temperate PFT. The strengthening of the isotopic coupling after removing the effect of δ 18 O sw from simulated δ 18 O R can be seen in Fontainebleau (Fig. 2).
Observed tree-ring width was significantly and positively correlated with ∆ 13 C R in Fontainebleau (Fig. 1d). This correlation 20 is common in pine and oak forests in the mid-latitudes of Europe and reflects that tree growth and leaf physiology are both similarly driven by water stress (Shestakova et al., 2019). ORCHIDEE (r898) largely overestimated this relationship because of an excessive drought sensitivity and lack of the representation of carbon remobilization, which would buffer the effect of the instantaneous leaf response on tree growth (Fig. 1e). Because of the latter process MAIDENiso captured the right magnitude of the expected coupling between leaf physiology and radial growth in the site (Fig. 1f).

Simulated relationship between isotopic variability and photosynthesis
Our results showed that the simulated interannual variability in δ 18 O R was significantly and negatively correlated with GPP across all sites in the climate gradient (r = -0.47 to -0.73, p < 0.01) because both source water (δ 18 O sw ) and leaf water enrichment (∆ 18 O lw ) related negatively with GPP (Fig. 8). Although ∆ 18 O lw dominates the correlation in most sites, it is the synergistic effect of δ 18 O sw variations that contributes to the apparent spatial consistency of the δ 18 O R -GPP correlation. The On the other hand, higher temperatures also lead to increased VPD, which increases leaf water enrichment and reduces GPP because of stomatal closure (Fig. 1d).
An earlier study attributed the correlation between δ 18 O R variability and satellite-based estimates of Net Primary Productivity (NPP) at large spatial scales to the common but opposite effect of VPD on leaf water enrichment and NPP, without considering the possible effect of source water (Levesque et al., 2019). Our finding suggests that δ 18 O R -productivity correlations should be interpreted with caution because they are not entirely driven by leaf physiology, and source water is often a 5 dominant driver of δ 18 O R (Fig. 6a).
The simulated interannual variability in ∆ 13 C R also correlated with GPP in most sites, but the relationship was less spatially consistent since its significance and sign varied across the boreal-to-temperate climate gradient (Fig. 8a). A positive correlation between simulated ∆ 13 C R and GPP in all temperate and boreal oak sites (r = 0.45 to 0.85, p < 0.01) is coherent with a stomatal limitation of carbon assimilation (Shestakova et al., 2019). In the two boreal pine sites, correlations were weaker and only the 10 northernmost site had a significant negative correlation (r = -0.40, p < 0.01). In such cool and moist environment, internal leaf CO 2 concentration is controlled by photosynthetic rate, which is limited by temperature and sunshine (McCarroll and Loader, 2004;Hilasvuori et al., 2009).
This contrast in ∆ 13 C R -GPP correlations between the temperate and boreal region mirrors a shift from a positive to a negative relationship between ∆ 13 C R and tree-ring width due to a change from water-limited to temperature-limited leaf

Conclusion
We demonstrated the potential of tree-ring width and carbon and oxygen stable isotopes to constrain the representation of tree 20 growth and physiology in the global land surface model ORCHIDEE (r898), bridging the long-standing gap between the treering and land surface modelling communities. ORCHIDEE had an overall performance to simulate tree-ring isotopic variability similar to that of the specialized MAIDENiso tree-ring model and the global vegetation model LPX-Bern. However, the lack of representation for the dependency of tree growth and carbon isotopes on carbon remobilization from storage in the model is currently a major limitation for the concurrent simulation of the ring width-isotope triplet (ring width, ∆ 13 C R and δ 18 O R ). 25 The large contribution of the isotopic signature of source water on δ 18 O R makes its simulation sensitive to the choice of the prescribed isotopic drivers. Future modelling studies would benefit from quantifying uncertainties on δ 18 O R variability from the isotopic drivers.
The simulated long-term physiological response of constant leaf c i /c a ratio under rising CO 2 during the 20th century is consistent with the observations in the temperate region, despite an overestimation of seasonal drought stress and the limitation isotopic version of the land-surface model ORCHIDEE: Implementation, evaluation, sensitivity to hydrological parameters, Hydrology