Modelling tree ring cellulose δ 18 O variations in two temperature-sensitive tree species from North and South America

Oxygen isotopes in tree rings (δOTR) are widely used to reconstruct past climates. However, the complexity of climatic and biological processes controlling isotopic fractionation is not yet fully understood. Here, we use the MAIDENiso model to decipher the variability in δOTR of two temperature-sensitive species of relevant palaeoclimatological interest (Picea mariana and Nothofagus pumilio) and growing at cold high latitudes in North and South America. In this first modelling study on δOTR values in both northeastern Canada (53.86 N) and western Argentina (41.10 S), we specifically aim at (1) evaluating the predictive skill of MAIDENiso to simulate δOTR values, (2) identifying the physical processes controlling δOTR by mechanistic modelling and (3) defining the origin of the temperature signal recorded in the two species. Although the linear regression models used here to predict daily δ18O of precipitation (δOP) may need to be improved in the future, the resulting daily δOP values adequately reproduce observed (from weather stations) and simulated (by global circulation model) δOP series. The δOTR values of the two species are correctly simulated using the δOP estimation as MAIDENiso input, although some offset in mean δOTR levels is observed for the South American site. For both species, the variability in δOTR series is primarily linked to the effect of temperature on isotopic enrichment of the leaf water. We show that MAIDENiso is a powerful tool for investigating isotopic fractionation processes but that the lack of a denser isotope-enabled monitoring network recording oxygen fractionation in the soil–vegetation–atmosphere compartments limits our capacity to decipher the processes at play. This study proves that the eco-physiological modelling of δOTR values is necessary to interpret the recorded climate signal more reliably.


Introduction
Oxygen isotopes in tree rings (δ 18 O TR ) are increasingly used as indicators of past climatic changes in temperate areas (Cernusak and English, 2015;Hartl-Meier et al., 2014;Saurer et al., 2008).They have been widely used to reconstruct past atmospheric conditions such as air temperature (Naulier et al., 2015), drought (Labuhn et al., 2016), precipitation amount (Rinne et al., 2013), isotopic composition of precipitation (Danis et al., 2006), relative air humidity (Wernicke et al., 2015), cloud cover (Shi et al., 2012), and even atmospheric circulation patterns (Brienen et al., 2012).This diversity of climatic targets possibly reconstructed based on oxygen isotopes hints at the challenge of understanding the complexity of the climatic and biological processes that control isotopic fractionation of oxygen in trees (Treydte et al., 2014).Uncertainties arise because different poorly measured factors influence δ 18 O TR values.Isotopic signals in tree rings cellulose are strongly influenced by isotopic signature of soil water taken up by the roots and by evaporative and physiological processes occurring at the leaf level and during downstream metabolism (Barbour et al., 2005;Gessler et al., 2014).Thus, a comprehensive approach that embraces existing mechanistic understanding of the fractionation processes involved is required.
A few isotopic process-based models have been developed to investigate the mechanistic rules governing the δ 18 O TR variations (Guiot et al., 2014): the Péclet-modified Craig-Gordon model (Kahmen et al., 2011) and the Roden's model (Roden et al., 2000) are able to estimate, at a daily time step, the δ 18 O values of soil and xylem waters as well as the isotopic fractionation occurring in the leaves due to evapotranspiration.Versions of these models are integrated in more complete forest ecophysiological models simulating the ensemble of forest water and carbon fluxes: (1) MAIDEN (Modeling and Analysis In DENdroecology) (Gea-Izquierdo et al., 2015;Misson, 2004), which contains the isotopic module MAIDENiso (Danis et al., 2012), and (2) MUSICA (Ogée et al., 2003(Ogée et al., , 2009)).Both account for important postphotosynthetic factors and are able to link photosynthesis and carbohydrate allocation to stem growth.
In this paper, we use the MAIDENiso model to decipher the δ 18 O TR variability in American temperature-sensitive species (Picea mariana in northeastern Canada and Nothofagus pumilio in western Argentina).The selected sites are of special interest for palaeoclimatology given that their δ 18 O TR chronologies carry strong temperature signals.A summer temperature reconstruction was already developed at the North American site (Gennaretti et al., 2017b;Naulier et al., 2015) and a calibration study conducted at the South American one highlighted the strong potential of δ 18 O TR values to reflect variations in summer-autumn temperatures over a large region south of 38 • S (Lavergne et al., 2016).However, up to now, the climate-δ 18 O TR relationships were analysed using a black box approach based on linear models.Here, we specifically aim at (1) evaluating the predictive skill of MAIDENiso to simulate δ 18 O TR values, (2) identifying the physical processes controlling δ 18 O TR by mechanistic modelling and (3) defining the origin of the temperature signal recorded in the two species.

Sampling sites and tree ring data
Two high-latitude American native species were studied here: (1) Picea mariana (Mill.B.S.P.; black spruce), which is a conifer widely distributed over the American boreal forest (Viereck and Johnston, 1990), and (2) Nothofagus pumilio (Poepp. et Endl. Krasser; lenga), which is an angiosperm deciduous species dominating the high-elevation forests along the Patagonian Andes from 35 to 55 • S (Donoso, 1981;Schlatter, 1994).We selected two sites of P. mariana in the centre of the Québec-Labrador Peninsula in northeastern Canada (L01 and L20; from 53 • 51 N-72 • 24 W to 54 • 33 N-71 • 14 W; ∼ 480 m elevation; see Gennaretti et al., 2014, andNaulier et al., 2014, for details) and three sites of N. pumilio in northern Patagonia, western Argentina (NUB, ALM and CHA; from 41 • 09 S-71 • 48 W to 41 • 15 S-71 • 17 W; 1270-1610 m elevation; see Lavergne et al., 2016Lavergne et al., , 2017, for details), for details).Climate in northeastern Canada is mostly continental and subarctic with short, mild and wet summers and long, cold and dry winter.Total annual precipitation averages 825 mm, with up to 46 % falling during the growing season in summer (June to September) (Naulier et al., 2014).In western Argentina, precipitation is largely concentrated from late autumn to early spring (May-November) followed by a drier and mild period during summer and early autumn (December-April) (López Bernal et al., 2012).
Four trees per site were collected for both species.The selection of the samples and analytical procedure for δ 18 O TR measurements were described in Lavergne et al. (2016) and Naulier et al. (2014).The developed δ 18 O TR chronologies covered the 1950-2005 and 1952-2011 periods at the northeastern Canadian and western Argentinian sites, respectively.The chronologies that were built for each species were significantly correlated between stands (Fig. 1).This supported the construction of a combined isotope chronology for both the northeastern Canada and western Argentina sites.

Modelling oxygen isotopes in tree ring cellulose with MAIDENiso
MAIDENiso is a process-based model that can simulate in parallel phenological and meteorological controls on photosynthetic activity and carbon allocation (Danis et al., 2012).It explicitly allocates carbohydrates to different carbon pools (leaves, stem, storage and roots) on a daily basis using phenological stage-dependent rules (see Gennaretti et al., 2017a, for details on the construction of the main MAIDEN model).
It also simulates the fractionation of carbon and oxygen isotopes during growth processes.In this study, the calculation of the daily δ 18 O TR in tree ring cellulose (‰) is based on the Danis et al. (2012) formulation of the Craig-Gordon model (Craig and Gordon, 1965): This equation summarizes how δ 18 O TR is determined by the following: i.The δ 18 O of the source (xylem) water (δ 18 O XW ), which is computed by averaging the δ 18 O SW values of the different soil layers weighted by the volume of water taken up by the roots in each layer.The isotopic effects of water mixing and soil evaporation on the δ 18 O SW values of the different soil layers are computed by a mass and isotopic balance (Danis et al., 2012).It is worth noting that fractionation occurs neither during water uptake by roots (Wershaw et al., 1966) nor during the transport of water from the roots to the leaves.
ii.The 18 O enrichment of the leaf water due to transpiration is described by (ε Craig and Gordon (1965), where a. ε * is the equilibrium fractionation due to the change of phase from liquid water to vapour at the leaf temperature (fixed at 21.4 iii.The biochemical fractionations (ε 0 ) due to oxygen exchange between carbonyl groups (C = O) in the organic molecules and water (DeNiro and Epstein, 1979;Farquhar et al., 1998).
iv.The dampening factor f o reflecting the exchange of the oxygen atoms between sucrose and xylem water during cellulose synthesis in the xylem cells of tree rings.
As previously evoked (i), δ 18 O XW of Eq. (1) depends on δ 18 O SW and thus on δ 18 O P values.However, long continuous time series of δ 18 O P are not available in the studied area.Here, we tested the impact of using two different methods for deriving δ 18 O P time series.
First, a linear model was used to estimate the daily values of δ 18 O P and subsequently δ 18 O V based on the primary drivers of their temporal variability (Dansgaard, 1964;Horita and Wesolowski, 1994), which are air temperature (T air ; • C) and precipitation at the corresponding site (P ; mm): with ε * T air the fractionation due to the change of phase from liquid water to vapour at the mean air temperature.The coefficients a and b were allowed to vary over a plausible range (or prior range) in the calibration process together with other MAIDENiso parameters, while coefficient c was fixed to a likely value (see Table 1 and Sect . 2.4).This estimated set of data is referred to in the following as the estimated δ 18 O P dataset.
Second, we run the model with the series of the daily δ 18 O P derived from two general circulation models (GCMs) with different spatial resolutions and enough available data at our site locations: (1) the Melbourne University model (MUGCM; Noone and Simmonds, 2002) forced by varying sea surface temperature (SST) from the HadISST dataset for the 1950-2003 period (2 • × 2 • resolution; hereafter referred to as the MUGCM δ 18 O P dataset), and (2) the Laboratoire de Météorologie Dynamique Zoom model (LMDZ5A; Hourdin et al., 2013;Risi et al., 2010) with the horizontal winds guided by those of the National Centers for Environmental Protection 20th Century Reanalysis (NCEP20) for the 1950-2008 period (Compo et al., 2011) (2.5 • × 3.75 • resolution; hereafter referred to as the LMDZ-NCEP20 δ 18 O P dataset).
The final δ 18 O TR time series are the annual average of the δ 18 O TR daily values (Eq. 1) weighted by the daily simulated stand gross primary production (GPP), assuming a proportional allocation of carbon to the trunk.For the northeastern Canadian sites, the GPP simulated by MAIDENiso was optimized using observations from an eddy covariance station (see Gennaretti et al., 2017a).Unfortunately, such observations were not available for N. pumilio, and therefore the parameterization obtained for the GPP of P. mariana was also used for the western Argentinian sites but constraining the simulations with phenological observations extracted from the literature.For example, to respect the annual cycle of the leaf area index (LAI) for N. pumilio (Magnin et al., 2014;Rusch, 1993), in MAIDENiso we used a seasonal LAI annual cycle with a development of leaves (LAI increase) between October and November, a maximum LAI (set at 5 leaf area/ground area) from November to April, a decreasing LAI (leaf fall) between April and May, and finally a leafless period (null LAI) from June to September (Magnin et al., 2014;Rusch, 1993).Furthermore, based on the finding that δ 18 O TR annual time series were more correlated with climate variables of specific months of the growing season (Lavergne et al., 2016), we also computed δ 18 O TR annual values by weighting the δ 18 O TR daily values (Eq. 1) with synthetic GPP time series maximizing the correspondence between observations and simulations.

Meteorological and atmospheric CO 2 data
At the western Argentinian sites, we did not have long daily records of observed climate data.Therefore, daily minimummaximum temperature and precipitation data were derived from the 20th Century Reanalysis V2c (2 • × 2 • resolution; Compo et al., 2011), which is one of the few reanalysis products entirely covering the 20th century.The temperature daily time series of the reanalysis were corrected in order to respect the monthly mean values detected at Bariloche, the nearest meteorological station from our sampling sites (∼ 48 km from the sites; 41  Hutchinson et al., 2009).In addition to these data we also used modelled daily data from the GCMs described above for both the western Argentinian and northeastern Canadian sites (see Table 2 with the input data used for each tested configuration).

Estimation of parameters influencing δ 18 O TR
We used a Bayesian method for the simultaneous calibration of the various MAIDENiso parameters specific to the study species and site.A set of 50 plausible blocks of parameters (posterior values) was selected according to the method described in Gennaretti et al. (2017a) using Markov chain Monte Carlo (MCMC) sampling (Table 1).The following prior plausible ranges were considered:  1).
3. The range for the kinetic fractionation ε k , which has been set to 26.5 ‰ in Farquhar et al. (1989) but that can vary over larger ranges (Buhay et al., 1996), was taken between 10 and 30 ‰ here.
4. The range for the dampening factor f o was allowed to vary between 0.3 and 0.5 following Saurer et al. (1997).
We tested the sensitivity of the MAIDENiso model to the calibrated parameters by modifying them within their respective prior calibration range.To control the robustness of the calibrated parameters, we performed the calibration of these parameters over two equal length intervals (1950-1977 and 1978-2005 for P. mariana; 1952-1981 and 1982-2011 for N. pumilio) keeping the second half for independent validation of the parameter estimates.Once the model was calibrated for the two species, the MAIDENiso model's performance to simulate P. mariana and N. pumilio δ 18 O TR interannual data was evaluated using the correlation coefficients (r) and the root mean square errors (RMSEs) between observed and simulated values.This is a standard approach to evaluate how well a mechanistic model is simulating δ 18 O TR variations (e.g., Danis et al., 2012;Lorrey et al., 2016).

Disentangling leaf-level fractionation processes and
source water influences on δ 18 O TR signature To define the relative contributions to the δ 18 O TR signature of the isotopic signal of the source water (xylem water) and of the fractionation processes due to transpiration taking place in the leaves, we designed two experimental simulations with MAIDENiso based on Eq. ( 1): 1. To quantify the influence of the variability in the isotopic composition of the xylem water on δ 18 O TR , we compared the reference simulations to those where the relative humidity (h air ) and the isotopic composition of atmospheric vapour (δ 18 O V ) were assumed to be constant.The constant values for h air and δ 18 O V were defined as the averages of the respective MAIDENiso outputs (h air = 0.62 and 0.9, and δ 18 O V = −26.28 and −17.34 ‰, respectively, for northeastern Canada and western Argentina; the XW source experiment simulation hereafter), 2. To quantify the influence of the isotopic enrichment of the leaf water due to transpiration on δ 18 O TR , we compared the reference simulations to those where the δ 18 O XW series were assumed to be constant.The constant value for δ 18 O XW was estimated as the average of the δ 18 O XW MAIDENiso outputs (δ 18 O XW = −13.81and −7.03 ‰, respectively, for northeastern Canada and western Argentina; the leaf water enrichment-driven experiment simulation hereafter).
Comparison between the experimental and reference simulations (i.e. using the optimal values of the parameters) was achieved through the calculation of the coefficient of determination (R 2 ).

Estimated versus modelled and observed δ 18 O P values
The modelled δ 18 O P series from the GCMs are similar to the GNIP datasets, with mean values ranging from −12 to −8 ‰ over June-September in northeastern Canada (Fig. S2a) and from −7 to −3 ‰ over December-April at the western Argentinian sites (Fig. S2b).In general, δ 18 O P series from LMDZ-NCEP20 model in western Argentina are slightly displaced toward higher values (+1 ‰) in comparison with the GNIP and MUGCM data.The estimated δ 18 O P values based on plausible values of coefficients a and b agree well with those of the models and observations in northeastern Canada.For the western Argentinian sites, they are 2-3 ‰ lower from April to October, i.e. late spring-early autumn (Fig. S2).1).The horizontal dot lines in (b) are the mean values of the observed δ 18 O TR .

Sensitivity of the model to the calibrated parameters
Most of the calibrated parameters have an influence on the correlations between observed and simulated δ 18 O TR series and/or on the mean levels of the simulated series (Fig. 2).The temperature and precipitation dependences of δ 18 O P values (respectively a and b coefficients) have the strongest influence on correlations.Increasing a and b values increase the mean δ 18 O TR levels, more strongly in western Argentina than in northeastern Canada (Fig. 2).Changes in the dampening factor (f o ) and in the biochemical fractionation (ε 0 ) have almost no effect on correlation, but their increase induces significant decrease of the mean levels of δ 18 O TR series.Finally, increasing the kinetic fractionation (ε k ) leads to lower correlations and to higher mean levels of δ 18 O TR (Fig. 2).

MAIDENiso performance in reproducing observed δ 18 O TR series
Split-period verifications of the calibrated relationships for P. mariana and N. pumilio when using estimated δ 18 O P series from Eq. ( 2) indicate that the calibration over either the firsthalf or the second-half periods provides similar posterior densities of the calibrated parameters than the ones obtained when calibrating over the whole periods (Fig. S3).One exception is observed in the calibration of coefficient a in northeastern Canada over the two half periods, where the posterior densities of a are different from the one obtained by calibrating over the entire period.Over the entire periods, observed and simulated δ 18 O TR series are significantly correlated in northeastern Canada (r = 0.56, p< 0.01 and RMSE = 0.67; Fig. 3a) and in western Argentina (r = 0.48, p< 0.01 and RMSE = 0.63; Fig. 3c).The correlations between observed and simulated δ 18 O TR series are slightly improved when we used synthetic daily GPP (r = 0.62 and r = 0.52, p< 0.01, respectively, for northeastern Canada and western Argentina; Fig. 3b and d).It is worth noting that the mean levels of the simulated δ 18 O TR series for the Argentinian sites are lower than those of the observations (offset of around −2.5 ‰; Fig. S4).The series were therefore corrected to respect the mean values detected in the observations (Fig. 3c and d).In contrast, the correlations between observation and simulation considerably decrease when we used modelled δ 18 O P from MUGCM models or LMDZ-NCEP20 reanalysis data.They only reach r = 0.13 (p> 0.05) to 0.23 (p< 0.05) in northeastern Canada and r = 0.23 to 0.26 (p< 0.05) in western Argentina, respectively (Fig. 4).
3.4 Influence of source water and leaf water isotopic enrichment to the δ 18 O TR signature The relative contributions to the δ 18 O TR signature of the isotopic signal of the source (xylem) water and of the 18 O enrichment of the leaf water due to transpiration were investigated.In both regions, the leaf water enrichment experimental simulations and the reference simulations have a higher correlation coefficient (R 2 centred on 0.9 and 0.95, respectively, for northeastern Canada and western Argentina; Fig. 5) than are the XW source simulations with the reference simulations (R 2 centred on 0.65 and 0.8, respectively, for northeastern Canada and western Argentina).This suggests that, with the model, the variability in δ 18 O XW has a weaker influence on δ 18 O TR variations than the changes of the leaf water isotopic enrichment do.Notably, P. mariana in northeastern Canada appears to be more sensitive to both influences than N. pumilio in western Argentina (Fig. 5).

Precipitation δ 18 O P variations and estimation
Although the regression models used to predict daily δ 18 O P values are likely too simplistic, the resultant monthly averaged values adequately reproduce the distribution of the observed (from GNIP stations) and modelled (by GCMs) (2) leaf water enrichment-driven experiment simulation (δ 18 O XW set as constant, green) in Québec (bold line) and Argentina (dashed line).
monthly δ 18 O P series in northeastern Canada.In western Argentina, the distribution of monthly δ 18 O P values is also well reproduced but the amplitude of variation in the predicted values is too high, leading to simulated values lower than the measured ones during the colder months.The temporal δ 18 O P variations are positively related to air temperature given the positive coefficient a.In agreement with the simple Rayleigh distillation model (Dansgaard, 1964), as air temperature decreases, the specific humidity at saturation decreases, and water vapour condenses.H 18 2 O condenses preferentially, and the residual water vapour becomes more and more depleted as condensation proceeds.Consequently, in the tropics, the 18 O / 16 O ratio in the meteoric water has been observed to decrease with increasing amount of precipitation and/or relative humidity (Rozanski et al., 1993).In extratropical regions, δ 18 O P may also correlate with precipitation amount (negative coefficient b), since both variables depend on the meteorological conditions.
The results of the linear regressions show comparatively lower influence of precipitation on δ 18 O P in western Argentina than in northeastern Canada (Table 1).This suggests that the imprint of the precipitation amount on δ 18 O P in western Argentina is low and that δ 18 O P variations are mainly controlled by seasonal changes in temperature, which is in agreement with previous work (Rozanski et al., 1995).However, due to the strong west-to-east precipitation gradient in this region (orographic rain shadow), large δ 18 O P variations occur over short distances (Rozanski et al., 1995;Smith and Evans, 2007;Stern and Blisniuk, 2002).Therefore, the daily precipitation dataset extracted from the gridded reanalysis data, which has a low spatial resolution (> 200 km), may not represent the daily variations in precipitation at a local scale faithfully.Therefore, the model may underestimate the contribution of precipitation on δ 18 O P variability in this particular area.
In contrast, in northeast Canada, both temperature and precipitation amount equally control the δ 18 O P variations.The high amount of precipitation falling in summer (∼ 46 %) should have a strong effect and decrease the δ 18 O P values in the condensed water, while high temperatures counteract this effect by increasing this ratio.Before reaching northeastern Canada, the air masses pushed by the dominant westerly winds discharge most of their humidity over the land, leading to a depleted δ 18 O P signal at our sites (for the same reason, δ 18 O TR values at L20, which is located 110 km northeast of L01, are ∼ 1 ‰ lower).Moreover, the δ 18 O P signal in the Canadian sites is comparatively more depleted than in the Argentinian sites, because of their higher latitude.It is worth noting that the resolution of the gridded meteorological dataset used for the Canadian sites is relatively high (∼ 10 km), which means that the local processes are likely well represented.

Relative performance in modelling δ 18 O TR values
The simulated δ 18 O TR series based on daily δ 18 O P estimation from the regression models reproduce the observations better than the ones based on δ 18 O P values derived from GCMs (Fig. 4).This is in part due to the greater number of parameters to optimize, as the calibration process can more easily find a solution that fits the observations better.This may however reflect error compensations especially in western Argentina where the estimated annual variability in δ 18 O P is too large.Conversely, in northeastern Canada, the annual variations in δ 18 O P that are estimated, simulated by GCMs and observed are in good agreement (Fig. S2).Although isotopeenabled atmospheric global models can reproduce the mean annual precipitation isotopic values and seasonality for many areas (Risi et al., 2010), results at specific sites, especially in mountainous regions such as at our western Argentinian site, can be less accurate (Fig. S2; see the offset between GNIP stations and LMDZ-NCEP20).Ideally, daily δ 18 O P longterm records from meteorological stations in the study region should be used as an input of MAIDENiso.Simulations from high-resolution regional circulation models, such as REMOiso, which has a 0.5 • × 0.5 • (∼ 55 km) horizontal resolution (Insel et al., 2013;Sturm et al., 2005Sturm et al., , 2007)), may produce reliable local δ 18 O P values.Such a dataset has proven to be quite helpful with MAIDENiso in the Fontainebleau forest (France) (Danis et al., 2012).However, up to now, measured or REMOiso δ 18 O P datasets in our regions of study do not exist, which is the case for most regions of the world.Moreover, early data (1970 from GNIP stations may have been compromised by pan evaporation and therefore isotopic enrichment.Therefore, we recommend that daily GNIP stations are set up in various forested ecosystems, that an effort is accomplished to homogenize older GNIP time series, and that high-resolution simulations of δ 18 O P are performed in wider regions. The modelling of δ 18 O TR values based on the estimation of δ 18 O P is relatively more accurate for northeastern Canada than for western Argentina (Fig. 3).As the mean levels of the measured δ 18 O TR values are high at the western Argentinian sites (mean value of about 30 ‰ ), the Bayesian optimization tends to increase the biochemical (ε 0 ) and kinetic (ε k ) fractionations as well as the coefficient a, while reducing the dampening factor (f o ) to reach more representative mean levels of the δ 18 O TR simulation.But still, these levels are too low in comparison with the observations (about 2.5 ‰ lower; Fig. S4).When the posterior value of a calibrated parameter is limited to the upper bound of the prior range of plausible values, as it is the case at the western Argentinian sites for a, b and ε 0 (Fig. S3), it means that either the prior range is too narrow, the model is inadequate, or some important process is not considered in the model.Here, the estimation of the prior ranges of both coefficients a and b were based on observed (GNIP stations) and simulated (GCMs) δ 18 O P values.Therefore, we expect their respective ranges to be consistent with local processes.When the prior range of a is extended to higher values in the optimization process, observed and simulated δ 18 O TR mean levels in western Argentina are better matching.However, in this case, the distribution of δ 18 O P values is shifted toward higher values, advocating for unrealistic estimated δ 18 O P variations.
One other possibility is that the prior range of ε 0 is too narrow.In accordance with DeNiro and Epstein (1981), Sternberg (1989) and Yakir and DeNiro (1990), the biochemical fractionation ε 0 is assumed here to be lower than 30 ‰.However, a recent study has demonstrated that this parameter, nearly constant between 20 to 30 • C, increases at lower temperatures to values of 31 ‰ (Sternberg and Ellsworth, 2011).During the growing season, maximum temperatures can reach 20 • C in western Argentina and 30 • C in northeastern Canada, which suggests that the high mean δ 18 O TR levels in N. pumilio may be due to biochemical fractionation higher than 30 ‰ due to temperature generally lower than 20 • C.However, when the prior range of ε 0 is extended to 31 ‰ in the optimization process, the mean δ 18 O TR levels of N. pumilio are still too low in comparison with the observations.These results advocate for the existence of other processes, which can explain this offset in mean levels in Argentina.For example, higher soil water evaporation than modelled by MAIDENiso should lead to less negative δ 18 O SW (and therefore δ 18 O XW ), which could explain the high mean levels of δ 18 O TR in Argentina.Caution should be exercised with such an interpretation since other species living in similar conditions as N. pumilio in western Argentina show comparatively lower mean δ 18 O TR levels than N. pumilio (i.e., Fitzroya cupressoides; see Lavergne et al., 2016).The ongoing monitoring and evaluation of isotopic processes based on synchronous measurements of vapour, precipitation, soil water and xylem water will certainly help understanding the high mean levels observed in Argentina, as well as increasing the representation of the involved processes in MAIDENiso.
The better fit between observed and simulated δ 18 O TR values obtained with specific forms of synthetic distributions of daily GPP for northeastern Canada and western Argentina (Fig. 3) suggests differential limiting factors in the two regions.The synthetic bimodal distribution of daily GPP with maxima in spring and autumn, as simulated in western Argentina, is often observed in a diversity of ecosystems such as in the Mediterranean environments (Baldocchi et al., 2010;Gea-Izquierdo et al., 2015).After the activation of the photosynthesis in early spring, increasing temperatures tend to be optimal for tree growth.However, in a modelling study, Lavergne et al. (2015) have shown that the influence of temperature on N. pumilio's growth becomes negative once a temperature threshold (soil moisture) is exceeded.Therefore, we assume that after reaching a threshold of temperature and soil moisture summer conditions, tree growth is inhibited, leading to a decrease of primary productivity.However, when temperature starts to decline and soil water supply tends to increase with increasing precipitation events, tree growth increases again until the end of the growing season.In contrast, because precipitation is more abundant in summer (June to September) in northeastern Canada (Naulier et al., 2014), high summer temperatures should be always beneficial to tree growth if enough soil water is available.Therefore, in agreement with GPP-derived eddy covariance data from the Fluxnet network (see Gennaretti et al., 2017a), a better fit between observations and simulations is observed when using a unimodal rather than a bimodal GPP distribution.Monitoring of tree physiology, environmental conditions and wood cell formation will provide a more detailed representation of the complex biological and ecological processes operating in Patagonia, allowing us to run the MAIDENiso model with better constraints.
4.3 What is the main origin of the temperature signal recorded in δ 18 O TR ?
The investigation of the relative contributions of the isotopic composition of the source (xylem) water and of the 18 O enrichment of the leaf water by transpiration on the simulated δ 18 O TR reveals that the variability in the former has a weaker influence on δ 18 O TR variations than that of the latter in North and South America.Therefore, the temperature signal recorded in δ 18 O TR series more likely reflects the effect of temperature on isotopic enrichment of the leaf water rather than on the isotopic composition of the source water.At the leaf level, air temperature has a strong effect on the relative humidity and therefore on the vapour pressure deficit (VPD), i.e. the difference between the saturation vapour pressure and the actual vapour pressure, which modulates the transpiration (Barbour, 2007).Thus, the imprint of the ambient air temperature on the fractionation processes occurring during transpiration is preferentially recorded in the tree rings of the two species.Furthermore, both the isotopic signature of the xylem water and of the fractionation processes occurring at the evaporation sites of the leaves have comparatively higher influence on δ 18 O TR in P. mariana than in N. pumilio.This is probably due to the lower amplitude of the day-by-day variations in the relative humidity in western Argentina (SD = 5 %) versus in northeastern Canada (SD = 16 %), which translates into a weaker influence of h air variations and therefore of leaf-level isotopic fractionation processes on δ 18 O TR values in western Argentina than in northeastern Canada.These results highlight the potential of MAIDENiso model to better refine the origin of the climatic signal recorded in the oxygen isotopic signature in the tree rings of different species.

Conclusions
Here, by using MAIDENiso model, we have provided a mechanistic overview of the climatic and biological processes controlling oxygen isotopic fractionation in two North and South American temperature-sensitive tree species.First, we have shown that using regression-based rather than model-based δ 18 O P estimates as inputs increases the predictive skills of our simulations, although this may be at the price of error compensations.Second, our study reveals that the variability in the isotopic composition of the source (xylem) water has a weaker influence on δ 18 O TR variations than that of the 18 O enrichment of the leaf water by transpiration.Last, these findings suggest that the imprint of temperature recorded in δ 18 O TR of the two species is likely related to the effect of temperature on isotopic enrichment of the leaf water.The isotopic monitoring of water within the soilvegetation-atmosphere compartments in future work will certainly provide the input and control data necessary to better constrain MAIDENiso.Our study demonstrates that the eco-physiological modelling of δ 18 O TR values is necessary and likely the only approach to accurately interpret the recorded climate signal.Based on the calibrations of MAID-ENiso presented here, the next step involves inverse modelling approaches to perform palaeoclimatic reconstructions in North and South America that are less biased by the complex and nonlinear interactions between climate, CO 2 concentrations and tree growth as recommended by Boucher et al. (2014).
Code and data availability.The code of the model can be found at https://doi.org/10.6084/m9.figshare.5446435.v1(recent version developed by Gennaretti et al., 2017).The daily δ 18 OP data from the MUGCM model were extracted via the SWING project webpage (SWING project, http://paos.colorado.edu/~dcn/SWING/database.php).The daily climatic data used for Québec were retrieved from the webpage of the Natural Resources Canada products (Natural Resources Canada, http://cfs.nrcan.gc.ca/projects/3/ 4), while those for Argentina were extracted from the Earth Sys-

Figure 2 .
Figure 2. Dependence of the correlation coefficients between observed and simulated δ 18 O TR series (a) and of the mean simulated δ 18 O TR levels (‰) (b) as a function of the range of calibrated parameters a, b, f o , ε 0 and ε k for the 50 simulations performed.The tests sites from Québec are in black and the Argentinean ones are in red.The vertical lines are the values of a plausible block of parameters retained in the MCMC optimization.The horizontal dashed lines are their respective 90 % confidence interval calculated with 50 simulations (see Table1).The horizontal dot lines in (b) are the mean values of the observed δ 18 O TR .

Figure 3 .
Figure 3.Comparison between observed (red or green) and simulated (grey) δ 18 O TR chronologies in Québec (a and b) and Argentina (c and d), respectively, using GPP (in g C m −2 day −1 ) simulated by MAIDENiso for each day of the year (DOY) (a and c) or synthetized for maximizing correlations (b and d).The simulations are based on estimated δ 18 O P series.The 50 different simulations inferred from the Markov chain Monte Carlo (MCMC) chains are in dark grey.The ±1 root mean square error (RMSE) range is represented in light grey.The mean correlation coefficients are significant at 99 % level ( * * ).

Figure 4 .
Figure 4. Density distributions of the coefficient of correlation (R) between observed and simulated δ 18 O TR chronologies in Québec and Argentina when the simulations are based on δ 18 O P series estimated by the regression model or from the MUGCM and LMDZ-NCEP20 models.

Figure 5 .
Figure5.Density distributions of the coefficients of determination (R 2 ) between the reference simulations and the (1) XW source experiment simulation (δ 18 O V and h air set as constant, black) and (2) leaf water enrichment-driven experiment simulation (δ 18 O XW set as constant, green) in Québec (bold line) and Argentina (dashed line).
Tree ring δ 18 O time series (‰) at the three sites in Argentina (NUB, ALM and CHA in dark grey) and two sites in Québec (L01 and L20 in dark grey; single trees in light grey).The bold black lines are the averaged values.The mean inter-site correlation coefficients are r = 0.60, p< 0.05 and r = 0.80, p< 0.01 in the South and North American sites, respectively.
In particular, it estimates, at a daily time step, δ 18 O values of soil water and xylem water, the isotopic fractionation occurring in the leaves due to evapotranspiration and the biochemical fractionation during cellulose formation.As input it uses daily maximum and Clim.Past, 13, 1515-1526, 2017 www.clim-past.net/13/1515/2017/ A. Lavergne et al.: Modelling tree ring cellulose δ 18 O variations

Table 2 .
Climate input data for all tested simulations.The prior ranges of the a and b coefficients in the equation of the daily δ 18 O P (Eq.2) were selected in order to get δ 18 O P values for each site consistent with the measured monthly local values from the nearest stations of the Global Network of Isotopes in Precipitation (GNIP), and with the simulated daily values from the LMDZ-NCEP20 model and from the MUGCM model (see Table