All models, from the simplest to the most complex, contain parameters that are often set by knowledge of the underlying processes and/or from information contained in the literature. When this information is not available, it is necessary to directly estimate the parameters by fitting the model simulations to data. This is the calibration step of modeling. The choice of an applicable calibration technique is necessary for successful model parameter estimation. Techniques range in complexity from manual calibration to algorithms such as global optimization, MCMC and data assimilation techniques. Manual calibration consists in a user inspecting the performance from a model and comparing it to the observations in order to determine the best parameters (Misson, 2004). This method is often time-consuming and sometimes difficult to implement. Global optimization techniques, which attempt to find the most accurate set of parameters with respect to data, are more automatic. Optimization is achieved by maximum likelihood approaches where the model parameters are adjusted by maximizing a likelihood function depending on the data and model parameters. See for example the PEST method in hydrology (Doherty, 2003). This provides accurate, reproducible parameters for a model, but lacks the ability to take into account the uncertainties of the model and data. Bayesian methods are useful to circumvent this problem. MCMC offer the opportunity to estimate the (a posteriori) parameter distribution once the calibration has been done, starting from their predefined (a priori) statistical distributions (Gelman et al., 1995; Andrieu et al., 2001). Such approaches can nevertheless lead to unrealistic results, in the event of inappropriate a priori information or non-relevant parameter distribution extractions. To bypass this difficulty, other kinds of procedures such as Particle Filtering (PF) have been recently proposed on the basis of genetic algorithms (Doucet et al., 2001; Pelletier et al., 2006). Unlike MCMC approaches, the procedure still consists in going though the parameter-space, but with a high number of initial choices, called particles, that are defined by their parameter values (their coordinates).

We present here the comparison of three methods, namely the PEST, MCMC and PF algorithms (Gaucherel et al., 2008a,b), to estimate 11 parameters of the MAIDEN model (Misson, 2004). The study is applied to Pinus halepensis ring width and xylem density series and climatic series from the Aix-en-Provence and Istres weather stations, all from Southern France (Nicault, 1999; Rathgeber et al., 2000b). The process-based model MAIDEN uses daily precipitation, minimum and maximum temperatures as input and was implemented for the 1950–1994 period. The model output of interest is the biomass allocated to the tree bole which is assumed to be comparable to the annual increment index based on wood density and ring width. The 11 parameters were calibrated on the 1962–1994 period using the Aix-en-Provence data and verified on the independent 1950–1961 period using the Istres station. In fact, several parameters were calibrated with daily physiological measurements, but we focus here on the dendrochronological calibration.

Figure 1 shows the interest of having a method able to provide uncertainties on the parameters. Some parameters have a Gaussian shape, which indicates that the parameter has an optimum value with a reasonable uncertainty, but the others are not really sensitive as they have a uniform distribution over a large interval: this means that these parameters can have relatively different values without really changing the simulation of the bole biomass increment. The PEST optimum values, the MCMC and PF probability modes of the eleven parameters are used in MAIDEN to estimate the bole biomass increments which are compared to the measured values. Figure 2 compares the three estimated curves with the observed curve for both the calibration period and the verification period. The agreement between the three methods and the observations is good, except for the 1956-1958 interval. In fact February 1956 was extremely cold and this seriously injured the trees which took 3 years to recover. This process is not included in the model so the 1956-1958 interval has been removed from the verification. The PEST optimization gives a r² = 0.46. This optimum fit is in good agreement with the MCMC final fit distribution, as the maximum MCMC fit value reached is r² = 0.52, while its modal value is 0.37. PEST needed only 16 iterations to reach this value. For the PF algorithm, the final fit gives a r² = 0.25 ± 0.02 with 1000 particles, a fit lower than that for the other calibration techniques. When applied to independent meteorological data, the three calibration parameter sets lead to rather similar bole biomass increment variations. Residuals between the standardized final simulations and the observations over the 11 years of validation show that the MAIDEN model has difficulty in simulating some years: namely the freezing years 1956 and 1985 and following, the rainy year 1963 and years 1993–1994 (Figure 2). When considering the residuals between the three final bole increment series and the corresponding observations, it appears that MCMC and PEST present rather similar behavior, distinct from that of the PF algorithm.

PEST is an extremely fast method, but its main disadvantage is that it does not provide any information on the uncertainty of the parameters. The main disadvantage of PF is that the particle number has to be high (from 10² to 10⁶) for the algorithm to be efficient, which is computationally costly. However PF has several advantages over MCMC and PEST (Gaucherel et al., 2008a): (i) it is easy to develop and to implement; (ii) it does not need any burning period (i.e., period necessary to reach equilibrium between inputs and outputs) nor parameter distribution filtering to calculate the a posteriori distributions; (iii) it can be parallelized; (iv) almost no tuning algorithm is necessary. As a sequential approach, the PF algorithm has the ability to indicate the highly constraining years, i.e., those which poorly fit simulations and observations. The final advantage of PF is that as it is based on the law of large numbers, it obtains smoother and more regular (Gaussian) parameter distributions, making it easier to analyze and interpret.

In a second step, using the parameter set obtained from MCMC, the model was applied to climate simulations provided by the ARPEGE climate model driven by the IPCC-B2 scenario (Gibelin and Deque, 2003) after statistical downscaling (Gaucherel et al., 2008b). The CO₂ atmospheric concentration has an influence on the climate and is taken into account by the climate model simulation, but it also has a direct effect through fertilization (already mentioned in section 2). Figure 3 shows the projection of Pinus halepensis growth for the next century with and without taking into account this fertilization effect (Gaucherel et al., 2008b). It appears clearly that ecophysiological models such as MAIDEN can be very sensitive to this effect. The main reason is that the model does not include any nutrient limiting factor, particularly to link the influence of nitrogen and photosynthetic capacity. Without the fertilization effect, tree growth will decrease because of lack of water (the precipitation deficit is amplified by more evapotranspiration), but if the latter effect is included, tree-growth will increase. As mentioned in the literature review in the introduction section, this fertilization effect is not always observed on living trees and is subject to controversy (Gedalof and Berg, 2010; Silva and Horwath, 2013).

Inverse modeling is a statistical method that can be used to estimate variables that are directly or indirectly related to the measured quantity. An inverse model differs from forward modeling in that it uses dependent variables (e.g., bole increment) to constrain input climate variables rather than using climate to predict dependent variable distributions. It is not an analytical inversion, but an iterative procedure which converges progressively toward the climate probability distribution compatible with the observations. The climatic space is randomly sampled to produce a large variety of climatic scenarios which are introduced in the vegetation model to simulate the corresponding vegetation composition and productivity. In dendroecological applications, the simulated bole increments from the process-based inversion are compared to the tree-ring data and those matching reasonably well are retained. These scenarios are used to build histograms, which are estimates of probability distribution functions of a climate able to generate such bole increments. This implies a measure of coherence between outputs of the model and tree-ring proxies.

Inversion of the model is done under the paradigm of the Bayesian theory (Robert and Casella, 1999). The main concepts used are prior and posterior probability distributions. The prior is the information, summarized in the form of a distribution, which is available before data analysis. The posterior is the information that will be deduced from the combination of data and the hierarchical model. In that respect, the hierarchical model is not restricted to the tree-growth model, but it is the function which relates the prior to the posterior. In statistical terms, it is the probability of climate C conditional to tree-ring data (Y). It is noted, according to Bayes' law:

where p(C) is the prior distribution of climate C and f(Y | C) is the likelihood function of tree-ring variables given the climate (as provided by the model from the climate). The prior is an initial guess of the probability distribution of the climate. It can be given by the knowledge we have from modern climate variability or, if we work on periods that are very different from the present one, from the knowledge which has been accumulated in palaeclimatology. The distribution law is then a uniform law defined on the range so defined. The likelihood function is provided by comparing tree-growth model simulations given a large set of climate scenarios to observed tree-growth data.

The above-mentioned inversion approach has been developed in the context of a multi-proxy reconstruction in Fontainebleau forest, near Paris, France (Boucher et al., 2014). At this location, different tree-ring proxies (ring widths, carbon and oxygen stable isotopes) have been measured and their sensitivity to climate has been demonstrated (Etien et al., 2008). Furthermore, from the modeling perspective, MAIDENiso was successfully parametrized and developed (in the forward mode) for this specific site (Danis et al., 2012). Thus, running MAIDENiso in the inverse mode allowed (Boucher et al., 2014) to retrieve, for each year since 1850 AD, a probability distribution of climatic inputs. Additionally, since MAIDENiso equations are constrained not only by meteorological conditions, but also by atmospheric CO₂ concentrations (e.g., via stomatal conductance modeling), it was possible to isolate the contribution of CO₂ to the overall reconstruction, and to point out a possible divergence effect associated with a CO₂ enriched atmosphere since the beginning of the industrial era. The idea behind the inversion performed at Fontainebleau was quite simple. For each year in which proxies were measured (but where climate was not available), we modified a reference year to fit the tree ring proxies as well as possible. The reference year (here 1984) was chosen to be as similar as possible to an “average year.” The modification of the reference year was performed using the delta method. Deltas are additive (temperature) or multiplicative (precipitation) values that modify the reference year, producing alternative meteorological scenarios that are iteratively injected into MAIDENiso. It is precisely the posterior distribution for these deltas that needs to be found by the inversion algorithm, for each year of the reconstruction.

As a first step in the inversion, a uniform prior is set on both temperature and precipitation delta values (Figure 4). A uniform prior means that only the bounds of paleo-temperature [−5,5] and precipitation [0.25,2] variability are fixed and all the values have equal probability. In a second step, using the MCMC procedure, we iteratively modify the delta values (called proposals) within the bounds of the uniform distribution. Proposals are accepted if two main conditions are met: (1) the proposal must be within the bounds of the prior distribution, and (2) the proposal must result in some gain in simulation accuracy (here, measured by a simple Euclidian distance between simulations and observations). A negative gain is accepted if it is above a small negative threshold. Otherwise, the proposal is rejected. The MCMC algorithm stops when the chain reaches a stable state, meaning that the successive accepted proposals are not correlated with one another (stationarity). All the accepted proposals are used to build a posterior distribution that represents the range of possible delta values required to simulate tree rings that resemble real-world proxies.

For each year of the reconstruction, delta posterior distributions differ because proxies, themselves, are different (Figures 4B,C). For example, the years 1960 and 1990 yield quite different summer temperature and precipitation estimates (black curves). Summer 1960 was relatively cold (median temperature delta = −3.4C) and wetter than normal (median precipitation delta = 1.35, i.e., +35%). In contrast, summer 1990 was warmer than normal (median temperature delta = 1.9C), while precipitation was close to average (median precipitation delta = 1.15, i.e., +15%). However, it is interesting to note that posterior distribution also changes depending on the CO₂ scenario used to constrain the inversion. In the preceding simulations, a realistic CO₂ scenario was used (Robertson et al., 2001), meaning that CO₂ variations are coherent with the increasing trend observed on the global scale [280 to 400 ppmv]. However, if CO₂ concentrations are fixed at preindustrial levels for the whole reconstruction period (1850–1999), then overall temperature reconstructions differ considerably in the low-frequency domain (red curves). This effect is partly visible on Figure 4C, where it is clear that the summer temperature reconstructed for year 1990 is colder by about 2C if pre-industrial CO₂ concentration (280 ppmv) is used instead of actual 1990 concentrations (370 ppmv). This difference becomes even more apparent when the full reconstructions are compared Figure 5. When pre-industrial CO₂ levels are used for the whole time period (red), a colder climate is reconstructed, and most importantly, the warming trend clearly observed in climatic records at Fontainebleau cannot be reproduced.

Although, as we have previously discussed, isolating the net effect of CO₂ on the forest response is likely to be more complex than the effect we modeled, our results suggest that CO₂ is an important driver for tree growth in Fontainebleau forest, and that the constraint of inversion with a realistic CO₂ scenario makes the reconstruction better converge toward observed climate variability. Abnormally low pre-industrial CO₂ concentrations force the model to compensate by too cold a climate. In the case of Fontainebleau forest, this compensation effect by a colder climate can be explained by links and feedbacks between climate, stomatal conductance and photosynthesis. With lower atmospheric CO₂ concentrations, stomata have to maximize carbon gains by increasing their area to 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.

The inversion of mechanistic dendroecological models offers a great framework for dendroclimatologists to incorporate and deal with non-stationarity (i.e., the so-called divergence D'Arrigo et al., 2008) and with interactions between climate, tree growth and CO₂ concentrations, among other things. This well-known ”divergence” problem which is related to the fact that the climate-tree-growth relationship was modified after the 1960s or later, is possibly an artifact of inadequate detrending (D'Arrigo et al., 2008; Esper and Frank, 2009) or an effect of CO₂ fertilization or change in water availability (Daux et al., 2011). The problem of tree-growth trend is often discussed in dendrochronology (Melvin and Briffa, 2008; Nicault et al., 2010). Trends in tree-growth time-series may be caused by long-term variations in climate, long-term CO₂ or nitrogen increases, and various factors related to the tree physiology and wood structure. Numerous techniques have been used to try to separate the climate trend from other factors, with more or less success. (Melvin, 2004) showed that the extensive use of concepts taken from tree-growth models and based on the assumption that common external growth forcing operates through its influence on photosynthesis may solve the problem. The example of model inversion shown in the previous section confirms the relevance of this approach.

Increasing atmospheric CO₂ concentration enhances biomass production in MAIDEN as in most ecophysiological models (Berninger et al., 2004). Under elevated CO₂ concentrations, stomatal area may be reduced without reducing internal carbon concentrations (Keenan et al., 2013). Partial stomatal closure increases water retention and increases biomass production per unit of water absorbed by the plant (better WUE). So, less precipitation is required to sustain growth rates. This rate of increase in WUE might be downgraded by various acclimation processes to long term CO₂ increases (Drake et al., 1997): a decrease in Rubisco, a rise in carbohydrate solubility, increases in light use efficiency, or lower rates of dark respiration. These processes are not really integrated in most of the models discussed. Clearly, additional research is required to better understand these processes so that they can be included and accounted for in simulations in both forward and inverse mode. On the one hand, the meta-analysis of Medlyn et al. (1999) showed that current models were adequate to model photosynthesis in elevated CO₂. On the other hand, (Berninger et al., 2004) concluded that additional processes such as those that control acclimation to CO₂ in the long term as well as other limiting factors (e.g., nitrogen availability) need to be integrated in upcoming model generations to improve the realism of climate reconstructions. A recent study (Li et al., 2014) that CO₂ fertilization over the past 50 years is too small to be distinguished in the ring-width data given ontogenetic trends and interannual variability in climate.

Finally, what is clear is that ecophysiological models of tree-growth still need more development, but that their use instead of oversimplified statistical methods (response functions or transfer functions) has enabled significant progress. It is clear that our understanding of mechanistic dendroecological processes and their integration in deterministic models is a continually evolving and rapidly changing area of research. For that reason, it would be unwise to apply the approach to a wide range of environments and species, without first going through a thorough examination of underlying dendroecological processes and achieving a clear assessment of the model's ability to simulate at least the most fundamental ones.

Data assimilation sensu stricto has been used extensively in various disciplines (meteorology, oceanography …), where many observational datasets, both in-situ and remotely sensed, are available. A spectacular application can be found in climatology with the so-called climate reanalysis (Kalnay et al., 1996). As ecology cannot call on such massive datasets, more appropriate techniques are those based on Bayesian statistics. Dubinkina et al. (2011) introduced particle filters in paleoclimatological modeling. We propose here to use a similar approach in dendroecology with MAIDENiso.

Models need to state the initial conditions to start the runs. MAIDENiso (Danis et al., 2012) needs several years of daily climate data (minimum and maximum temperature, precipitation). For simplicity, we consider the current and previous years. The scheme of the method is presented in Figure 6. We start the simulations at the year t₀ + 1, using climate at time t₀ + 1 and t₀. By introducing small perturbations in this 2-year climate vector, we generate an ensemble of N simulations of tree proxies (tree-growth and isotopes) for year t₀ + 1. These simulated proxies are the particles which need to be compared to the observed proxies. A weight w_i, proportional to the likelihood of the model state (simulation of the proxies), conditional to the observations, is attributed to them according to their similarities to the observations (in Figure 6 the closest particle has the largest bullet). A posterior probability distribution is thus generated (2):

where p(c | d) is the posterior distribution of the climate c given the proxies d, i.e., the probability distribution of the reconstructed climate; k is a normalization constant; H is the measurement model (MAIDENiso simulation); is the observation error; c is the input climate. H(c_i) is the modeled proxies which depend on the climate c_i; the weight w_i measures the distance between the i^th particle and the observed value. The posterior distribution is resampled, providing M new equal weight particles by propagation using the model. These particles are compared to the observed proxies at year t₀+2. The procedure is repeated until the end of the observed proxies series.

This is the theory, but some difficulties remain to be solved before the method is operational. First the input climate is a large dimension vector. The model needs daily values of three variables, so that the number of possible combinations of the parameters is much larger than the number of particles. In reality, daily values can be derived from a relatively small number of parameters, namely the parameters of the probability distribution of the seasonal temperatures and precipitations. These large dimension vectors can also be reduced by a few principal components. Another problem is the comparison of the simulation H(c_i) to the proxies d. In reality we have here three proxies (bore increment, carbon and oxygen isotopes). They must be combined in Equation 2 and they must be simulated in the same units by the model. For the isotopes, this is the case, but the bore increment—related to the area and the density of the ring in a section, i.e., which is not a mass—must be compared to the annual biomass increment simulated by the model. In a previous section, this was solved by standardizing simulated and observed values. This is not the only solution and maybe not the best one.

This scheme of data assimilation must be considered as an avenue of research. To be useful, this approach must be applicable to long series, to provide climate reconstructions over several centuries of the past and projections to the future. As at each step we produce climate reconstructions and these reconstructions are used in the following step, there is a danger of drift. It is not certain that the link between observed and simulated proxies is sufficiently strong to foil the possible drifts. Studies on the stationarity of the process must then be conducted. A more complex but promising approach is given by coupling a climate model like LOVECLIM (Goosse et al., 2010) to a tree-growth model like MAIDENiso (Danis et al., 2012) and by assimilating data to this coupled system. In this case, the observed proxies are used to constrain the climate simulations to make them more realistic. As the climate model is global, this needs a network of proxies and not only a single site, and also a downscaling method (Vrac et al., 2007) to adapt the climate simulations to the dimension of the studied forests.

In contrast to empirical-statistical models, forward models of tree-growth can explicitly account for non-climatic factors, which reduces the risk of attributing to climate some tree-growth changes that are in fact related to internal processes. This extends their applicability and reduces the risk of violating the uniformitarian principle (which assumes stationarity of the relationships between tree-rings and abiotic factors). These forward models have another advantage: they are excellent tools to simulate proxies in various abiotic conditions and hence to better understand their sensitivity to various factors and evaluate biases and errors associated with reconstructions of past climate. These approaches do not imply that present empirical-statistical methods should be adandoned, but rather that our toolkitsshould be extended to analyze the proxies. The approaches used in dendroclimatology were developed 30 or 40 years ago, and were creative and productive responses to the situation we faced then. We now have a much richer set of records, a more diverse and powerful set of tools for analyzing and using them, and benefit from great advances in understanding tree functioning and climate variations on the long term.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.