An inverse modeling approach for tree-ring-based climate reconstructions under changing atmospheric CO 2 concentrations

Over the last decades, dendroclimatologists have relied upon linear transfer functions to reconstruct historical climate. Transfer functions need to be calibrated using recent data from periods where CO 2 concentrations reached unprecedented levels (near 400 ppm – parts per million). Based on these transfer functions, dendroclimatologists must then reconstruct a different past, a past where CO 2 concentrations were far below 300 ppm. However, relying upon transfer functions calibrated in this way may introduce an unanticipated bias in the reconstruction of past climate, particularly if CO2 has had a noticeable impact on tree growth and water use efficiency since the beginning of the industrial era. As an alternative to the transfer function approach, we run the MAIDENiso ecophysiological model in an inverse mode to link together climatic variables, atmospheric CO 2 concentrations and tree growth parameters. Our approach endeavors to find the optimal combination of meteorological conditions that best simulate observed tree ring patterns. We test our approach in the Fontainebleau Forest (France). By comparing two different CO2 scenarios, we present evidence that increasing CO2 concentrations have had a slight, yet significant, effect on the reconstruction results. We demonstrate that realistic CO2 concentrations need to be inputted in the inversion so that observed increasing trends in summer temperature are adequately reconstructed. Fixing CO 2 concentrations at preindustrial levels (280 ppm) results in undesirable compensation effects that force the inversion algorithm to propose climatic values that lie outside from the bounds of observed climatic variability. Ultimately, the inversion approach has several advantages over traditional transfer function approaches, most notably its ability to separate climatic effects from CO2 imprints on tree growth. Therefore, our method produces reconstructions that are less biased by anthropogenic greenhouse gas emissions and that are based on sound ecophysiological knowledge.


Introduction
With climatic change and rising CO 2 concentrations, the need to place recent trends in a multisecular perspective has led to an unprecedented expansion of the field of paleoclimatology.Tree rings have played a central role in that development, offering the possibility to reconstruct various climatic fields at high temporal resolution (annual to intra-annual) (e.g., Fritts, 1976;Schweingruber, 1996) and over large areas (regional (e.g., Cook et al., 2004), hemispheric (e.g., Esper et al., 2002;Briffa et al., 1998a;Mann et al., 1999) or global (e.g., Mann et al., 2008)).To perform such reconstructions, dendrochronologists have typically relied upon various black-box approaches (also called transfer functions) (Fritts, 1976) that approximate relationships between tree rings and climate by a time-invariant function.However, evidence for nonstationarity and divergence in the biological response to climate have called into question the transfer function approach (Vaganov et al., 2006).Time-varying tree-ring-to-climate responses may have various origins (D'Arrigo et al., 2008): thresholded, complex and nonlinear responses, changes in resources and nutrient availability (Reich and Hobbie, 2013;Silva and Anand, 2013), É. Boucher et al.: Inverse modeling of tree rings alteration of climate-growth relationships by CO 2 and sustained changes in water-table levels.If these phenomena are not accounted for in the reconstruction model, discrepancies between real and reconstructed paleoclimatic variability can be expected (Briffa et al., 1998b).
Since carbon dioxide is well mixed and homogeneously distributed throughout Earth's atmosphere (Myhre et al., 1998), CO 2 might be one of the most widespread sources of divergence, yet its biological impact on climate-to-treegrowth relationships remains imprecise (D'Arrigo et al., 2008).Since the beginning of the industrial period, CO 2 concentrations have increased by about 150 ppm (parts per million), reaching 400 ppm in 2013 (Tans and Keeling, 2013).Whether or not this increase has had a significant effect on tree growth however remains vigorously debated.On one hand, several studies that aimed at detecting and attributing the various signals contained in tree rings argued against a possible long-term fertilization effect (e.g., (Gedalof and Berg, 2010;Silva et al., 2010;Girardin et al., 2011)).On the other hand, a whole set of empirical experiments (see reviews by Norby et al. (1999), Norby et al. (2005) and more recently the landmark paper by Keenan et al. (2013)), though conducted on smaller timescales, have concluded that CO 2 has had important effects on tree productivity and water use efficiency.However, eddy-covariance measurements cannot be directly compared with tree ring series to precise the role of CO 2 on overall productivity and water use efficiency, since the former does not allow discriminating between trees, understory vegetation, and soil respiration/evapotranspiration.From a paleoclimatological perspective, this debate is fundamental.If tree ring chronologies have been imprinted by a strong CO 2 signal since the late industrial period (LaMarche et al., 1984), the appropriateness of transfer models needs to be called into question, as it cannot be assumed that recent tree-ring-to-climate relationships reflect those of the preindustrial period, a period where CO 2 concentrations were lower than at present.
To circumvent the possible biases caused by divergent relationships, several alternatives exist, most of which make use of extensive process-based modeling.Data-assimilation approaches such as the one proposed by Goose et al. (2010) aim at constraining climate models in order to identify plausible and coherent climatic scenarios that match the variations of natural proxies.Data-assimilation techniques allow testing various forcing scenarios, but in the case of tree rings, the linkages between climate simulations and proxies remain imprecise as signal and noise components cannot be clearly distinguished.In other fields of paleoclimatology, several authors have put forward the idea that deterministic models can be reversed so that the set of rules and mathematical equations representing the many processes that occur in ecosystems can be used in an inverse mode to retrieve probabilistic estimates for input meteorological conditions (Chuine et al., 2004;Hughes and Ammann, 2009;Guiot et al., 2009Guiot et al., , 2000;;de Cortázar-Atauri et al., 2010;Evans et al., 2013).This ap-proach has been tested by Garreta et al. (2010) for pollen data, yet no dendroecological model was inverted to perform a reconstruction.Tolwinski-Ward et al. (2011) have recently accomplished a first step towards the reconstruction and used a simplified version of the Vaganov-Shashkin model (VSlite) that mimics the principle of a limiting factor to simulate thresholded ring-width responses.However, while a simple water bucket model is efficient computationally, some of the most important ecophysiological processes responsible for tree growth in nature (e.g., stomatal conductance, photosynthesis, resource allocation and partitioning of carbon, reuse of the preceding year's photosynthetates, etc.) are omitted, with uncertain consequences for the inversion.Additionally, in VS-lite, the focus is on a single proxy (ring widths) while it has been demonstrated that a multiproxy approach is the most efficient way to enhance long-term environmental signals (McCarroll and Loader, 2004;Mann, 2002).Variables such as oxygen and carbon isotopes in tree-ring cellulose are now routinely measured; so inverse modeling approaches need to account for these proxies that may present complementary signals and different sources of noise (McCarroll and Loader, 2004).
Here, we present a reconstruction performed from the inversion of the MAIDENiso model (Danis et al., 2012), an updated version of the MAIDEN model (Misson, 2004) that accounts for carbon and oxygen isotope fractionation.MAID-ENiso is a full ecophysiological model, meaning that the links between meteorological (input) variables and (output) tree ring variables (e.g tree ring widths, carbon and oxygen isotopes) are not empirically parametrized as in a regression, but explicitly represented as mathematical equations within the model.The main objective of this paper is to demonstrate the scientific and conceptual advantages that inverse dendroclimatological modeling may have for the study of past climates.With this objective in mind, we will address the following questions: (1) how good are the reconstructions performed from the inversion?(2) Are they better than traditional transfer function techniques based, for example, on linear modeling?(3) Is the approach able to take into account and eventually isolate the effect of CO 2 on tree ring growth and thus attest of its impact on past climates reconstructed from tree ring series?This paper will be divided in two parts.In the first part, we rapidly review the model's capability as well as its parametrization.We then present our inversion approach explicitly.In the second part, we apply the inversion technique to reconstruct past climatic variations in the Fontainebleau Forest (near Paris, France), where Danis et al. (2012) have recently calibrated the model.

MAIDENiso
MAIDENiso is a process-based, dendrogeochemical model that can simulate various tree growth parameters useful in  (Gaucherel et al., 2008;Misson et al., 2004), as well as Belgium (Misson, 2004).Recently, the model was completed by Danis et al. (2012) to capture the main variations of δ 18 O and δ 13 C of tree ring cellulose found in sessile oak (Quercus petraea Matt.).In the forward mode, it uses relatively simple meteorological inputs: daily precipitation (cm d −1 ), minimum and maximum temperature ( • C).However, other variables are also required (on a daily basis) such as atmospheric concentration of CO 2 (ppm), the content of δ 18 O of precipitation, and δ 13 C content in atmospheric CO 2 .Lastly, a limited amount of parameters need to be found in the literature or estimated directly from tree rings using optimization techniques (see Gaucherel et al. (2008) and Misson et al. (2004) for more explanations).MAIDENiso captures and models site-specific phenological as well as meteorological controls on transpiration, stomatal conductance, respiration and photosynthetic production.For an extensive review of the processes and how they are modeled, readers are referred to the original works by Misson (2004) and Danis et al. (2012).The algorithm includes several allocation rules that distribute biomass production, carbon and oxygen isotopes in various reservoirs such as leafs, bole, roots, and storage.From a dendroecological perspective, bole is the most important reservoir, as it can be readily compared to standard tree ring measurements such as tree ring widths, δ 18 O and δ 13 C.However, storage is another important property of the model as it allows simulating processes responsible for autocorrelation and persistence in tree growth parameters, therefore enabling the use of carbohydrates to be passed to the next year of simulation.This aspect has important implications in ecology but also in paleoclimatology.For example, in the case of the Fontainebleau Forest (described later) Danis et al. (2012) argued that an adequate modeling of storage reserves by MAIDENiso allowed for a proper simulation of the 1976 drought in that area.In the years following the drought, and in spite of the amelioration of climatic conditions, trees were unable to store a sufficient amount of reserves, which resulted in sustained growth reductions.From a paleoclimatological perspective, it is important that the inversion distinguishes between climatic and storage-related effects.Otherwise, prolonged ecophysiological responses could be mistaken for persistent harsh climatic conditions.

MAIDENiso's inversion procedure
The inversion approach described here (Fig. 1) allows to retrieve paleoclimatic information from ecophysiological models.The method is presented with examples from the Fontainebleau Forest (France).First, it must be stated that inverse modeling procedures do not aim at finding a backward solution to the set of rules and equations that are implemented in a process-based model.Instead, inverse modeling approaches seek for the optimal combination of input data so that outputs (i.e., simulations) are as close as possible to the observations.In the precise context of dendroecological modeling, inverting a model such as MAIDENiso implies finding the optimal combination of meteorological inputs that will best simulate tree ring parameters (tree ring widths, δ 18 O and δ 13 C).Consequently, the model has to be run repeatedly with different meteorological scenarios each time, preserving only the ones that simulate tree ring parameters that resemble observations, and eliminating the rest.
So, at first glance, a contradiction seems to emerge: the model has to be run with meteorological data to produce tree ring simulations, but, in the case of a climate reconstruction, it is precisely those conditions that are unknown.Therefore, the general strategy proposed here is to define, from the modern data set, a reference or average meteorological year (AMY) that will serve as a basis to produce alternative meteorological scenarios (Fig. 1).Identifying the AMY ensures that a certain degree of temporal coherence exists in the relationship between temperatures and precipitation.In the case of the Fontainebleau Forest, the meteorological data set originally used by Danis et al. (2012) corresponds to the average of two stations (Glandée at Villiers-en-Bière and Faissandière at Fontainebleau) located 10km north of the forest.From this data set and with a simple measure of euclidian distance, we identified 1984 as the year that is closest, in terms of average temperature and precipitation (spring and summer) to the 1953-2000 average (Fig. 2a, b). Figure 2c and d present temperature and precipitation variations that occurred in 1984 at Fontainebleau.The temperature curve presents a typical hyperbolic shape with an average temperature of 10.2 • C. The precipitation series depicts considerable variability at the daily to monthly scales, with low summer precipitation and precipitation peaks in June (4.5cm) and September (3.8 cm), respectively, and with an average of 0.27 cm precipitation per day.
Once the AMY is identified, the next step is to find a way to modify it in order to produce alternative meteorological scenarios (AMS) that will be iteratively reinjected within MAIDENiso, through the undermentioned inversion algorithm.The strategy proposed here is to modify the AMY's original temperature and precipitation series with constants (thereafter named deltas, or ) that generate AMS.For each year simulated, we defined two different : one for temperature (δ τ ) that is additive and a second for precipitation (δ π ) that is multiplicative.They both apply during the growing For each year simulated, two modify the AMY: one for temperature (δ τ t ) and one for precipitation (δ πt ).The resulting AMS is passed to MAIDENiso, along with atmospheric CO 2 concentrations and atmospheric δ 13 C to simulate tree growth parameters.Simulated tree growth is compared to observations, through the metric λ, which has to be maximized (Eq.1).The Metropolis-Hastings algorithm iteratively modifies the value of δ τ and δ π so that it converges towards stable states (Eq.2).season, which is defined here as the period between early April (Julian day 91) and late October (Julian day 304).Figure 2c and d present two alternative meteorological scenarios computed from two different δ τ and δ π values.In the first scenario (red), δ τ and δ π were assigned a value of 2, meaning that all days are incremented positively by 2 • C, and that all precipitation events, when they occur, are doubled in magnitude.In the second scenario (blue), we assigned values of −2 and 0.5 to δ τ and δ π , respectively.Therefore, 2 • were subtracted to the temperature series and the daily amount of precipitation was halved during the growing season.δ τ and δ π were implemented so that the modifications of meteorological conditions by the Metropolis-Hastings algorithm (described hereafter) are kept as simple and straightforward as possible.However, the modifications themselves must be confined within bounds that respect the properties of local climatology in order to avoid generating implausible meteorological conditions.Uniform priors are then set to δ τ and δ π , with values bounded between −5 and +5 for the former and between 0.25 and 4 for the latter.Such values produce bounds that are larger than the modern  interannual meteorological variability, but were set uniform because we assumed to have no a priori knowledge on past fluctuations within these bounds.Moreover, Gaussian priors centered around AMY would give too much weight to the "normal" conditions and are likely to be incompatible with our wish to generate sufficiently different climatic conditions from the present ones.
The Metropolis-Hastings random walk algorithm (MHA) is used to estimate, for each year t of the reconstruction the posterior probability distribution of δ τ (t) and δ π(t) .The MHA is a Markov chain Monte Carlo (MCMC) procedure and therefore sample candidates produced within the chain depend only on the current sample value.The chain progresses towards a stable state by accepting or rejecting proposed samples.Concerning our application, the posterior distribution of δ τ (t) and δ π(t) corresponds to the distribution of accepted values that best simulate tree ring parameters, at year t.
The MHA is extensively described in Hastings (1970); here, the most important principles will be presented and exemplified with our tree ring application (Fig. 1).For each meteorological year t to be reconstructed from tree ring parameters, let f ( t ) be a target density that is proportional to the desired probability distribution of at year t, P ( t ).First, f ( t ) executes a MAIDENiso run through which tree ring parameters are simulated given input AMS.Then a metric λ (Fig. 1) needs to be maximized [−∞ → 0] and correspond to the negative squared distance between simulated (s) and real-world observations (o): 1) To calculate λ, all tree ring variables are transformed into dimensionless standardized scores.From there, the MHA starts by defining current values for δ τ t (i) and δ π t (i) (hereafter grouped into t (i) ) that fall within the bounds of a uniform prior distribution q( t | t (i) ).The algorithm then defines a candidate set of deltas: * t .The chain M walks towards the next iteration, accepting or rejecting the candidate values according to the following rules: ( In other words, if f ( * t ) ≥ f ( t (i) , then the candidate set * t is more likely than the original set of deltas to produce tree ring parameters that mimic those found in nature.In that case, * t is accepted with a probability of 1 and the algorithm "walks" to the next iteration replacing t (i+1) by * t .In the case where f ( * t ) < f ( t (i) , there are two possibilities.First, there is a chance equal to α that * t is accepted and passed to the next iteration.Second, there is a 1-α chance that * t is rejected and consequently, t (i) is maintained at the next iteration.After a large number of iterations, the chain M will tend to be proportional to P ( t ) because the algorithm will "intuitively" visit high-density regions more often than it does for low-density regions.
The inversion scheme presented earlier (Fig. 1) enables sampling in the posterior distribution P ( t ) one year at a time.However, tree ring series usually present a certain degree of autocorrelation.Such a persistence may be caused by climatic variations, but may also result from persisting insufficient carbon storage.Fortunately, MAIDENiso can account for storage and remobilization of carbon reserves within the tree, but to do so, it needs to simulate more than a single year at a time.In order to correctly model autocorrelation, we opted for a 4-year block simulation (Fig. 3).Thus, to obtain a distribution for t , we fed MAIDENiso with t−1 , t−2 , and t−3 , which are passed from the preceding simulations (dashed lines on Fig. 3).Hence, all four respectively modify the AMY (year 1984) to produce 4-year-long daily temperature and precipitation series that are used by MAIDENiso.Accordingly, simulated tree ring series exhibit a 4-year dependance structure that should be similar to the one observed in real world samples.Median δ τ (t) and δ π (t) values are passed to the subsequent year and so forth (Fig. 3).It is important to underline that the optimization of the metric λ (Eq. 1) by the MHA is exclusively performed for the last year of the simulation, irrespective of the number of years included in the simulation.

Application of the inversion algorithm to the
Fontainebleau Forest (France)

Tree ring data
All tree ring data (Fig. 4) used for the inversion are latewood measurements previously published by Etien et al. (2008).Thirty living, dominant oaks were sampled in the Fontainebleau Forest (two stations, N = 15/station).On each tree, three cores were sampled at 1.3 m height.Each series was measured and cross-dated by the phytoecology team of the National Institute of Agronomical Research (INRA, Nancy, France) using a master chronology constructed from more than 400 oaks sampled in the nearby area.The living tree chronology was standardized using the adaptive regional growth curve technique (Nicault et al., 2010) to produce an average chronology (Fig. 4a).Tree rings were then cut using a scalpel to separate latewood from earlywood.Samples were pooled and then milled with an 80 µ sieve.α-cellulose was extracted using the Soxhlet method (Leavitt and Danzer, 1993).Isotopic compositions (Fig. 4b, c) were determined with a Carbo Erba elemental analyzer coupled to a Finnigan MAT252 mass spectrometer (at LSCE, Gif/Yvette, Fr).Correction of cellulose δ 13 C series for the Suess effect was not necessary since, in MAIDENiso, atmospheric δ 13 C is used as  For each year t where the inversion is conducted, AMS from the three preceding years is also given as input to MAIDENiso, so the simulation is performed on a total of 4 years (bold arrows).Once the best combination of is found by the Metropolis-Hastings algorithm, the mean δ values are passed to the next iteration.An example is given by red dashed lines.s found for year t − 2 are also used to simulate t − 1 and t.
an input and the model directly yields depleted series comparable to those found in tree rings.
The sensitivity of tree ring proxies to climate was examined by Etien et al. (2008).Briefly, Fontainebleau's oaks respond to both summer temperature and precipitation, however at different levels (Table 1).Latewood widths respond positively to summer precipitation while the correlation with summer temperature is not significant (Etien et al., 2008).This feature was also highlighted by Kelly et al. (2002) and recalls that, in temperate areas, oak growth is mainly controlled by moisture availability.δ 13 C and δ 18 O variability are both positively correlated to summer temperature and negatively related to summer precipitation (Etien et al., 2008).Correlations with δ 18 O are slightly stronger, but the fact that both isotopic proxies are sensitive to the same climatic variable suggest that, at Fontainebleau, stomatal closure is an important response mechanism to drought.

Atmospheric data
In addition to tree ring and meteorological data, three additional variables were prescribed as inputs to MAIDENiso: atmospheric CO 2 concentrations, atmospheric δ 13 C content -Atmospheric CO 2 concentrations.We retained two scenarios that will help us identify CO 2 effects and associated biases within the reconstruction (Fig. 5).
The first scenario ("A1") mimics the reality and depicts a progressive increase in CO 2 concentrations towards modern values.The nearest CO 2 record (Gifsur-Yvette) covers the 2000-2007 period.Only the mean annual cycle was retained from that short series and was superimposed over the long-term trend in CO 2 concentrations extracted from ice cores (Robertson et al., 2001).The latter reflects global rather than local trends in CO 2 .Therefore, both series were combined into a single one and simply extrapolated at the daily time step as in Danis et al. (2012).Scenario A1 was also shifted by +15 ppm by comparison to the ice core series to take into account observable differences between the two series during their common period.The second scenario ("A2") represents a fictive situation into which CO 2 concentrations have remained stable (at preindustrial levels, i.e., 280 ppm) over the full reconstruction period.
-Atmospheric δ 13 C content.Atmospheric δ 13 C data (Fig. 5) was retrieved from published ice core series (Francey et al., 1999) and covers the full reconstruction period.Modern annual cycles were extracted from the Schauinsland data set (Black Forest, Germany) (Schmidt et al., 2003;Danis et al., 2012) that extends back to ca. 1976.As for CO 2 trends, the Schauinsland annual cycle was superimposed over the long-term trend extracted from the ice core series.The resulting series was shifted by 1 ‰ towards lighter ratios to match values observed during the common period at Schauinsland.δ 18 O of precipitation.In MAIDENiso, δ 18 O p is modeled statistically using linear regressions, with daily temperature and precipitation used as regressors (Danis et al., 2012).Regression parameters were calibrated using 42 years of δ 18 O p (dependant variable), temperature and precipitation data extracted from the REMO(REgional MOdel)iso mesoscale climate model (Sturm et al., 2005).The regression equation was directly incorporated into MAIDENiso's code so that the modeling of δ 18 O p could be achieved at the daily time step.Such parameters are site-specific.

Method of inversion-run details
In this study, a Metropolis-Hastings analysis were performed using the R statistical software (R Core Team, 2012), most particularly through the metrop( ) function of the mcmc package (Geyer and Johnson, 2012).For each year t simulated by the MHA, 500 burn-in samples were generated.Also, to ensure that the chain M is stable and composed only of successively uncorrelated and independent candidates, we picked only one t candidate at each 100 accepted samples.
In total, we sampled 1000 candidate values for t .We finally verified the chain convergence by checking that the lag-1 correlation between successive candidate values remained below 0.35.

Results and discussion
Here we present a paleoclimatic reconstruction for the Fontainebleau Forest area, based on the inversion of the MAIDENiso model.The reconstruction extends back to ca. 1850 and focuses on summer (JJA) temperature and precipitation.The A1 CO 2 scenario is used, unless the contrary is mentioned.

1960-1999 period
A comparison between reconstructed and observed summer climatic conditions is presented in Fig. 6 for the 1960-1999 period.The visual correspondence between both variables is fairly good, but slightly better for precipitation, both in terms of high-frequency and low-frequency variability.Reconstructed temperatures correlate well (r = 0.53) with observations (Table 1, left column), but precipitation correlates even better (r = 0.67).As is the case for observed chronologies (Table 1, right column) the temperature signal is weaker than the precipitation signal.In other words, the climatic signal recorded by trees seems to be well extracted by our reconstruction.Moreover, as the model was already parametrized by Misson (2004) for tree rings and by Danis et al. (2012) for C and O isotopic fractionation, we argue that this climate reconstruction is coherent with the main ecophysiological processes that govern tree growth at Fontainebleau.In the present simulation, a single δ τ (t) and δ τ (t) value now modifies the AMY to generate AMS but, ultimately, the comparison is made only for the summer period, owing to the fact that all proxies were measured on latewood.Consequently, all AMY are modified identically between day 91 and day 304, with no distinction between spring and summer conditions.Accordingly, by increasing the number of parameters to be estimated one would certainly augment the resolution of the ecophysiological modeling, therefore allowing seasonal reconstructions to be performed.Seasonal integration would enable the temporal integration of spring, summer and autumn growth processes and would yield even more sound reconstructions.Paleotemperature modeling would also certainly benefit from such an approach because even more important for tree growth than absolute temperatures is the length of the growing season.For instance no parameters modify the length of the growing season, which is fixed to exactly 214 days.Further modeling efforts will be invested in order to increase the number of parameters to take into account varying growing season length as well as seasonal differences between spring and summer.However, it remains to be demonstrated that the Metropolis-Hastings algorithm can still converge when the number of parameters is increased.

1850-1999 period
We extended the reconstruction back to AD 1850 and compared it to the Climate Research Unit (CRU) gridded temperature (CRUTEM4, Jones et al., 2012) and precipitation (Hulme et al., 1998) reconstructions (Fig. 7).The general agreement between both series is fairly good, with, again, a clear visual correspondence in the high-and low-frequency domains.The correlations with the CRU data are 0.51 and 0.40, respectively, for temperature and precipitation (n = 100  .Median reconstructed values correspond to the black curve, observations in red.The quantiles between 2.5 and 97.5 fall within the gray shading. for precipitation and n = 150 for temperature).The temperature reconstruction seems to be in phase with the CRU data and the correlation is similar to that calculated for the 1960-1999 period.(Fisher r to z score = −0.15,p > 0.05).There is, however, a clear drop in the strength of the correlation for the precipitation (Fisher r to z score = 2, p < 0.05).While this drop could be interpreted as a flaw of the method, some other aspects must be considered.First, CRU data correlates poorly with precipitation at Fontainebleau over the 1960-1999 period (r = 0.37), so there is no reason to assume that both series did correlate well in the past.Second, it is clear from the visual examination of both curves that low frequencies are well reconstructed.Such low frequencies may relate to the regional trend.However, high-frequency variations of precipitation are difficult to interpolate at the local scale, especially in the Paris area where a temperate-oceanic climate dominates, with an important influence of convective precipitation during the summer.This kind of precipitation may be hard to interpolate and, therefore, divergence can be expected with observations.

Does the reconstruction perform better than a linear transfer function?
Inversions were compared to the reconstructions obtained from calibrating transfer functions (multiple linear regression) of the 1960-1999 period (Fig. 8).Although both techniques produce temperature and precipitation reconstructions that correlate well with one another, (r = 0.74 and 0.84, respectively), several differences need to be underlined.to an underestimation of the variance of temperatures for the last 50 years.Taking the CRU (temperature) variance as a reference for the pre-1950 period (0.9), the variance obtained by inversion is 0.91 and the transfer function produces a comparable series (variance of 0.89).However, for the 1950-2000 period, the CRU variance augments to 1.16.The inversion follows that trend (post-1950 variance of 1.12) while the transfer function does not and the variance remains low (0.81) after 1950.This is a very important point, as variance losses and underestimation have plagued the transfer function approach over the last decade (Bürger, 2007) and may suggest that model inversions have a slight advantage over transfer functions from this point of view.For precipitation, the two methods produce series with comparable variances for both pre-1950 (0.90 ± 0.1) and 1950-2000 periods.Nevertheless, the transfer function methods creates a slight but significant (p < 0.01, n = 150) positive trend in the amount of summer precipitation.Such a trend leads to an underestimation of summer precipitation before 1900 and to a slight overestimation from 1950 onward.Such a trend does not exist either in the inversion or in the CRU precipitation series over the last century (p > 0.05).
In addition to reproducing variance and trends that better characterize Fontainebleau's climate variability since AD 1850, the inversion approach seems to have other conceptual advantages.By contrast to traditional transfer functions, no calibration is required, except for the initial parametrization of the ecophysiological model (here, MAIDENiso).Consequently, we expect the method to be more adaptable to the modeling of nonstationary tree-ring-to-climate relationships.The inversion procedure also builds on sound ecophysiological knowledge and therefore the links between climate and tree growth are not reduced to a set of calibrated parameters, but are modeled with regards to the complexity and multiplicity of processes that control vegetation growth.From a modeling perspective, it is clear that each and every improvement made within the plant ecophysiology science community and implemented in deterministic models will directly translate to better reconstruction performances, and a greater adaptability of the model to a wide range of environments.

Is the reconstruction affected by changes in atmospheric CO 2 concentrations?
MAIDENiso is an ecophysiological model that simulates tree growth parameters, given atmospheric inputs such as concentrations in CO 2 and δ 13 C. Therefore, the model can be used as a testbed to evaluate the impact, among other things, of CO 2 increases on the reconstruction.As mentioned earlier, we performed two reconstructions, each forced with a different CO 2 scenario (Fig. 5): scenario A1, which corresponds to the typical anthropogenically modified curve; and A2, which represents the nonanthropogenic curve (CO 2 remains stable).The A1 and A2 simulations yield divergent temperature reconstructions (Fig. 9).Both reconstructions correlate well with one another (r = 0.75 for temperature), implying that high-frequency variations are well reconstructed using both CO 2 scenarios.However, the series exhibit different longterm trends.Temperatures reconstructed using the A1 scenario show a clear rise towards warmer conditions (a trend that is also observed in instrumental climatic series) while temperatures reconstructed using the A2 scenario slightly decreased over time.Contrastingly, from the point of view of precipitation, the A1 and the A2 scenarios produced comparable reconstructions at both high and low frequencies (r = 0.87).
Our results suggest that the inversion must be forced by a realistic and increasing CO 2 scenario (e.g., A1) to adequately reconstruct low frequency temperature variability.Black curves and associated 95 % uncertainty range are median reconstructed values obtained when the inversion is constrained by realistic, increasing CO 2 data.The red curve and associated 95 % confidence intervals are median reconstructed values obtained when the inversion is forced when CO 2 is fixed at preindustrial levels.Observed climate variability  is in dashed green.Linear lines with corresponding colors were added to each series to evaluate multidecadal trends in reconstructed and observed climates.
This is because, at least at the Fontainebleau Forest, MAID-ENiso is able to balance the interplay between climate, atmospheric CO 2 concentrations, tree growth and isotope fractionation, even when the interaction between these variables is not stationary through time.As a consequence, without inputting a realistic CO 2 scenario to the inversion, the MHA compensates by altering climate s in order to bring AMS outside from the bounds of observed variability and converge towards tree ring simulations that match with observations.The A2 reconstruction perfectly exemplifies this compensation phenomenon.When the inversion is forced by fixed, preindustrial CO 2 concentrations, stomata remain open in order to maximize carbon gains and simulate growth variations in a way that mimics observed tree ring series (i.e., measured from oak trees that grew in a CO 2 -enriched atmosphere).However, with stomata wide open, evapotranspiration rates increase accordingly and reconstructed temperature needs to drop significantly (Fig. 9) in order to reduce water losses at the leaf surface.It is important to recall that this compensation mechanism originates from the fact that, no matter which CO 2 scenario is chosen for the inversion, the MHA systematically "walks" towards the best way to modify climate inputs in order to keep λ as close to zero as possible.Ultimately, we could verify that, because of this compensation mechanism, simulated values for each tree ring proxy and for each scenario are strongly related to the original observations (Fig. 4, dashed lines).
The impact of CO 2 on the inversion to be discussed in light of the mechanisms that govern CO 2 exchanges at the leaf-atmosphere interface.In MAIDENiso, stomatal conductance is modeled using the Ball et al. (1987) model, later modified by Leuning (1995), which is based on the early works of Wong et al. (1979).The Ball-Berry-Leuning model expresses stomatal conductance as a function of relative humidity, photosynthesis and CO 2 concentrations at the leaf surface.The equation however implies a nearly constant ratio of internal (C i ) to atmospheric (C a ) CO 2 concentrations, which is maintained stable through stomatal aperture and closure.A direct consequence of the Ball-Berry-Leuning model is a quasi-parallel response of internal (C i ) to increasing CO 2 (C a ).For that reason, its applicability to very high CO 2 concentrations (e.g., 700-800 ppm), without accounting for acclimation mechanisms such as those mentioned in Drake et al. (1997) (decrease in Rubisco, augmentation of carbohydrate solubility, increases in light use efficiency, lower rates of dark respiration or nitrogen limitation (Berninger et al., 2004)) is questionable.Our results suggest that the Ball-Berry-Leuning model might however be appropriate for the range of CO 2 concentrations that prevailed since (and even before) the preindustrial era.A corollary interpretation is that acclimation mechanisms must not have had a large influence over the last 150 years, at least in the Fontainebleau Forest's context.Otherwise, the inversion would not have converged towards realistic climate scenarios when forced by the A1 scenario.These findings ultimately support the conclusions of Medlyn et al. (1999), which underly that in several environments, acclimation processes have not completely overridden the net carbon gains related to CO 2 .
One of the strengths of our study is that the inversion of MAIDENiso is constrained by multiple proxies simultaneously, rather than by a single one, therefore better reflecting the array of physiological processes that govern tree growth.This is a very important point, especially in a context where the inversion is used to attest of the impact of CO 2 on tree-ring-based climatic reconstructions.It is well known that each proxy responds differently to rising CO 2 concentrations, each with its own strengths and source of bias.For example, tree ring widths potentially record changes in tree productivity, but low-frequency growth trends attributable to rising CO 2 concentrations might be removed by standardization (Melvin and Briffa, 2008;Cecile et al., 2013), in which case the inversion algorithm (like any other reconstruction method) would converge towards unrealistic trends that are exempt of any CO 2 imprint.Similarly, it has been argued by Silva and Horwath (2013) that, due to an artifact of calculation, δ 13 C chronologies might systematically overestimate the response of carbon isotopes to rising CO 2 concentrations, therefore creating artificial increases in reconstructed water use efficiency.However, as mentioned by McCarroll and Loader (2004) and also by Silva and Horwath (2013) for the particular case of CO 2 effects, the solution might lie in the use of multiple proxies, as it is proposed here.Although we cannot ignore the possible source of bias associated with each proxy, our approach certainly enables the simulation of all three proxies (from reconstructed climate) in a way that mimics real world observations (Fig. 4).Indeed, a large number of climatic and CO 2 combinations can explain variations in each proxy individually, but a much smaller range of combinations can explain covariations consistent with all three proxies.At the Fontainebleau site, it was previously shown by Etien et al. (2008) that stomatal conductance has a drastic importance on the overall growth, as evidenced by negative latewood width responses and covarying δ 13 C and δ 18 O.In such water-sensitive environments, tree rings are particularly responsive to CO 2 effects (Knapp et al., 2001;Owensby et al., 1999;Huang et al., 2007), therefore reinforcing the indispensable role of CO 2 on tree growth and isotope fractionation at our site and, indirectly, attesting of the necessity to include realistic CO 2 scenarios to produce unbiased climatic reconstructions.
Ultimately, these results are appealing for the study of divergence in dendroclimatology (D'Arrigo et al., 2008) and suggest that the inversion of process-based dendroecological models is a suitable approach to reconstruct climate under changing CO 2 concentrations.As it was demonstrated earlier, our analysis shows that tree-ring-to-climate relationships are significantly affected by changes in atmospheric CO 2 concentrations.Failing to include a realistic atmospheric CO 2 scenario triggers an undesirable climatic compensation effect that forces the inversion to diverge from a realistic climatic reconstruction.Of course, the quality of the inversion will always be dependent on the quality of tree ring series used as reference.In this regard, the development of novel "signal free" standardization techniques (e.g., Melvin and Briffa (2008) or Cecile et al. (2013)) will help constrain the inversion with tree ring chronologies that contain as little bias and/or noise as possible.In addition, from an ecophysiological modeling point of view, the amelioration of equations and rules implemented within ecophysiological models from both field and experimental studies will contribute to expand the applicability of the approach to a variety of environments and species.Finally, the improvement of random-walk algorithms and the design of better parametrization strategies will directly translate into more robust and more precise reconstructions.

Conclusions
In this study, we demonstrated that the MAIDENiso ecophysiological model can be used to perform paleoclimatic reconstructions that take into account the effect of anthropogenic CO 2 concentrations on tree growth.We provided a way to inverse that model through the use of a random-walk Metropolis-Hastings algorithm that iteratively converges towards the best combination of meteorological conditions to fit observations.The approach was exemplified in a reconstruction of past summer temperatures and precipitation at the Fontainebleau Forest (France).The novel reconstruction technique seems to perform well in that environment, and presents significant conceptual advantages compared to the regression technique.First, since the inversion is constrained by an ecophysiological model such as MAIDENiso, the complex and nonlinear interactions between climate, CO 2 concentrations and tree growth are mechanistically determined, and not reduced to a set of calibration parameters.Second, an ecophysiological model such as MAIDENiso can account for the temporal dependance of tree-ring-to-climate relationships on changing CO 2 concentrations.It is therefore an appealing and promising tool for the study in divergence in dendroclimatology and for the reconstruction of climate in the presence of nonstationary tree-ring-to-climate relationships.
Our study shows that atmospheric CO 2 is an important variable at the Fontainebleau Forest, and as a consequence, it is necessary to constrain the inversion with a realistic CO 2 scenario in order to perform reconstructions that converge towards observed climate variability.Running the inversion with CO 2 concentrations fixed at preindustrial levels (280ppm) forces the model to propose climate values that are outside from the bounds of observed variability.In the Fontainebleau context, this compensation effect leads to the reconstruction of a colder climate.This compensation mechanism can be explained by links and feedbacks that exist between climate, stomatal conductance and photosynthesis, all of which are represented by MAIDENiso under the form of mechanistic equations.With lower atmospheric CO 2 concentrations, stomata need to remain wide open in order to maximize carbon gains and match real-world growth rates.However, by doing so, water losses by evapotranspiration and evaporation at the leaf surface also increase.So, in order to compensate for water losses, reconstructed temperatures need to be lower, a situation that ultimately minimizes evapotranspiration and evaporation fluxes.
This analysis remains a first step towards a more generalized, large-scale utilization of ecophysiological models in www.biogeosciences.net/11/3245/2014/Biogeosciences, 11, 3245-3258, 2014 dendroclimatology.As suggested by Tingley et al. (2012), process-based models such as MAIDENiso can link instrumental and proxy observations in a framework that accounts for nonstationary and nonlinear relationships.Because of that particularity, their integration in a Bayesian hierarchical modeling is highly desirable and can improve the quality of climate-field reconstructions while opening up new and stimulating perspectives on the modeling of the spatiotemporal dependence between proxies and observations.In that context, our work supports the ideas developed by Tingley et al. (2012) and exemplifies a case where a mechanistic model actually integrates, within the reconstruction, nonstationary tree-ring-to-climate relationships and their dependance on changing CO 2 concentrations.In years to come, the performance of the process-based inversion approaches will likely increase, as a result of the amelioration of equations and rules implemented within ecophysiological models, the improvement of random-walk and inversion algorithms, and the design of better parametrization strategies.

Figure 1 .
Figure 1.MAIDENiso's inversion flowchart.For each year simulated, two modify the AMY: one for temperature (δ τ t ) and one for precipitation (δ πt ).The resulting AMS is passed to MAIDENiso, along with atmospheric CO 2 concentrations and atmospheric δ 13 C to simulate tree growth parameters.Simulated tree growth is compared to observations, through the metric λ, which has to be maximized (Eq.1).The Metropolis-Hastings algorithm iteratively modifies the value of δ τ and δ π so that it converges towards stable states (Eq.2).

Figure 2 .
Figure 2. Yearly and daily (1984) summer temperature and precipitation series for the Fontainebleau Forest.The left panel presents temperature (a) and precipitation (b) series, with the AMY (1984) indicated by red squares.Average conditions are in dashed red.The right panel presents daily variations (in 1984) for temperature (c) and precipitation (d).Black curves represent observations in 1984 (AMY).The red curve corresponds to the AMS that results from the modification of the AMY by a δ τ of 2 and a δ π of 2. The blue curve represents an AMS that originates from the modification of the AMY by a δ τ of −2 and a δ π of 0.5.

Figure 3 .
Figure3.Flowchart of the autocorrelation modeling.Each gray box corresponds to the inversion algorithm presented in Fig.1.For each year t where the inversion is conducted, AMS from the three preceding years is also given as input to MAIDENiso, so the simulation is performed on a total of 4 years (bold arrows).Once the best combination of is found by the Metropolis-Hastings algorithm, the mean δ values are passed to the next iteration.An example is given by red dashed lines.s found for year t − 2 are also used to simulate t − 1 and t.

Figure 4 .
Figure 4. Observed and simulated tree ring proxies.Black bold lines represent observed latewood widths (LW), δ 18 O and δ 13 C of tree ring α cellulose.Black dotted lines correspond to the simulations performed from the A1 (increasing CO 2 ) scenario, while red dotted lines represent the A2 (stable CO 2 ) scenario.

Figure 5 .
Figure 5. Daily atmospheric CO 2 (a) and δ 13 C (b) data used in the inversion.The A1 (black) scenario is a realistic one while scenario A2 (red) describes a hypothetical case where CO 2 concentrations have remained at preindustrial (1850) levels.

Figure 6 .
Figure 6.Comparison with modern climatological records(1960- 1999).Median reconstructed values correspond to the black curve, observations in red.The quantiles between 2.5 and 97.5 fall within the gray shading.

Figure 7 .
Figure 7. Summer temperature and precipitation reconstructions back to ca. 1850.Reconstructed climatic conditions (black) are compared to the CRU (green) data during the common period.Observations during the 1960-1999 period are in red.The quantiles between 2.5 and 97.5 fall within the gray shading.

Figure 8 .
Figure 8.Comparison between two reconstruction techniques: the inversion of MAIDENiso (black) and the transfer function (dashed blue) calibrated for the 1960-1999 period.

Figure 9 .
Figure9.Summer temperature and precipitation reconstructions at Fontainebleau, France, obtained from the inversion of MAIDENiso.Black curves and associated 95 % uncertainty range are median reconstructed values obtained when the inversion is constrained by realistic, increasing CO 2 data.The red curve and associated 95 % confidence intervals are median reconstructed values obtained when the inversion is forced when CO 2 is fixed at preindustrial levels.Observed climate variability is in dashed green.Linear lines with corresponding colors were added to each series to evaluate multidecadal trends in reconstructed and observed climates.
dendroecological analysis such as tree ring widths, δ 18 O and δ 13 C of tree ring cellulose.Earlier versions of the model were used to simulate tree growth in northern and southern France Biogeosciences, 11, 3245-3258, 2014 www.biogeosciences.net/11/3245/2014/

Table 1 .
Correlation between (left column) observed and reconstructed summer climatic conditions at the Fontainebleau Forest for the 1960-2000 period and (right column) correlation between tree rings and observed climatic conditions.