Modelling the climatic drivers determining photosynthesis and carbon allocation in evergreen Mediterranean forests using multiproxy long time series

. Climatic drivers limit several important physiological processes involved in ecosystem carbon dynamics including gross primary productivity (GPP) and carbon allocation in vegetation. Climatic variability limits these two processes differently. We developed an existing mechanistic model to analyse photosynthesis and variability in carbon allocation in two evergreen species at two Mediterranean forests. The model was calibrated using a combination of eddy covariance CO 2 ﬂux data,

bon allocation within plants (Breda et al., 2006;Niinemets, 2007;Chen et al., 2013). Many factors affect the different physiological processes driving forest performance. Among them, the net effect of the rising CO 2 mixing ratio ([CO 2 ]) and climate change is meaningful when determining the forests' capacity of acclimation to enhanced xericity (Peñuelas et al., 2011;Keenan et al., 2011;Fatichi et al., 2014). Forest process-based models have been developed to mimic these mechanisms. They can include different levels of complexity but generally implement calculations of leaf photosynthesis upscaled to the canopy and carbon allocated to different plant compartments (Le Roux et al., 2001;Schaefer et al., 2012;De Kauwe et al., 2013). Although there is evidence that the tree performance depends to some extent on stored carbohydrates (Breda et al., 2006;McDowell et al., 2013;Dickman et al., 2015), these models have received some criticism when used to understand plant performance in response to climate change. This is in part because they are C-source oriented, therefore can exhibit certain limitations to represent the C-sink hypothesis (i.e. that growth rates are limited by environmental factors such as water stress, minimum temperature or nutrient availability rather than by carbohydrate availability) and address dysfunctions related to the tree hydraulics (Millard et al., 2007;Breshears et al., 2009;Sala et al., 2012;Körner, 2013;McDowell et al., 2013;Fatichi et al., 2014).
Complex process-based models benefit from multiproxy calibration, particularly when such data are applied on different spatio-temporal scales (Peng et al., 2011). The temporal scale can be approached using time growth series of dendrochronological data. However, the analysis of the past always adds uncertainties related to the influence of unknown stand conditions to properly scale productivity. Flux data, including stand productivity, can be estimated using the eddy covariance technique . These data overcome many of the limitations of dendroecological data (e.g. intra-annual resolution, control of stand conditions and scaling of net productivity), but they lack their spatial and temporal coverage. Thus, CO 2 flux data can be used to implement unbiased models of canopy photosynthesis and can then be combined with dendroecological data to study how carbon is allocated to stem growth as a function of environmental forcing (Friedlingstein et al., 1999;Chen et al., 2013, McMurtrie andDewar, 2013).
Mechanistic models can also be used to analyse the environmental factors determining instability in the climategrowth response (D'Arrigo et al., 2008). Different processbased models have been applied with dendroecological data used either in forward or inverse mode (see Guiot et al., 2014, for a review). Among these models, the process-based model MAIDEN (Modeling and Analysis In DENdroecology).  was originally developed using dendroecological data. The model explicitly includes [CO 2 ] to calculate photosynthesis (hence its influence on carbon allocation) and includes a carbohydrate storage reservoir, this be-ing one of its strengths compared to other models Sala et al., 2012;Guiot et al., 2014). It has been previously employed to analyse growth variability in one temperate and two Mediterranean species Gaucherel et al., 2008) and recently, in inverse mode (also including C and O stable isotopes), to reconstruct past climate . However, it requires further development to ensure that it provides unbiased estimates of forest productivity and assesses uncertainties in the response of trees to climatic variability on a greater spatial scale at the regional level. In particular, its parameterization would need improvement if the model is applied to assess how climate modulates forest performance and the pattern of C allocation within plants (Niinemets and Valladares, 2004;Fatichi et al., 2014).
In this study we use multiproxy data to develop a processbased model and investigate how evergreen Mediterranean forests have modified stand photosynthesis and carbon allocation in response to interacting climatic factors and enhanced [CO 2 ] in the recent past. The first objective was to develop a process-based model based on MAIDEN . Within the new version of the model, photosynthesis, carbon allocation, canopy turnover and phenology are now calculated using climate-explicit functions with a mechanistic basis. The model is adapted to give unbiased estimates of canopy photosynthesis and stem growth using instrumental data. Specifically, within the new model formulation, (1) photosynthesis is penalized by prolonged water stress conditions through reductions in leaf area index (LAI) and maximum photosynthetic capacity; (2) the pattern of carbon allocation is directly determined by soil water content (i.e. water stress) and temperature through nonlinear relationships; (3) these relationships can be contrasting for different phenophases and affect photosynthesis and the pattern of C allocation independently. Once the model was developed, a second objective was to analyse how [CO 2 ] and climatic variability affect the temporal instability in annual forest productivity, water use efficiency and carbon allocation. We hypothesize that they will exhibit differences in their long-term variability in relation to recent climate change driven by different functional acclimation processes within trees.

Study sites and climatic data
The study sites were two evergreen Mediterranean monitored forests in Southern France where CO 2 , water vapour and energy fluxes are measured using the Eddy covariance technique   (Rambal et al., 2004;Limousin et al., 2012). Both forests grow on rocky and shallow soils that have a low retention capacity and are of Jurassic limestone origin. The climate is Mediterranean, with a water stress period in summer, cold or mild winters and most precipitation occurring between September and May. Meteorological data were obtained from the neighbouring stations of St. Martin de Londres (for Puechabon) and Aubagne (for Fontblanche). According to those data Puechabon is colder and receives more precipitation than Fontblanche (Table 1). Meteorological data showed a decrease in total rainfall since the 1970s in Puechabon but no trend in Fontblanche. Both sites exhibit a positive trend in temperatures more evident for the maximum values (Fig. A1). We assumed that GPP (gross primary productivity) is driven by the top pine and/or oak layers and that the percentage of LAI related to the understory shrub layer will behave like that of the oak species (evergreen, shrubby). For Fontblanche we considered a maximum leaf area index (LAI max ) of 2.2 m 2 m −2 (3 m 2 m −2 plant area index, PAI), composed of 70 % pine and 30 % oak (Simioni et al., 2013). For Puechabon we considered an LAI max of 2.0 m 2 m −2 (2.8 m 2 m −2 PAI) monospecific to Q. ilex (Baldocchi et al., 2010;Limousin et al., 2012). The specific leaf area (SLA) considered was 0.0045 m 2 g −1 for Q. ilex and 0.0037 m 2 g −1 for P. halepensis (Hoff and Rambal, 2003;Maseyk et al., 2008).

The model
We used MAIDEN , a stand productivity mechanistic model driven by a number of functions and parameters representing different processes. The model inputs are precipitation, maximum and minimum temperature, and [CO 2 ] with a daily time step. This model has been previously implemented for monospecific forests, including two oaks and one pine species, using dendroecological chronologies of growth and, when available, stand transpiration estimates from sap-flow sensors Gaucherel et al., 2008). However, the model has never been compared to actual CO 2 flux data to ensure that it provides unbiased estimates of forest productivity. In this study, the model was further developed to match ground-based observations and generalize model use by modifying the photosynthesis and allocation modules (including the different phenophases) in relation to climatic drivers. To properly scale model outputs and get unbiased estimates of stand productivity, we used CO 2 eddy covariance fluxes . Different parameters were calibrated to different data sources, including some species-dependent and some site-dependent parameters, as follows. The transpiration rate (E) of day i is calculated us-ing a conductance approach: E(i) = g s (i) × VPD(i)/P atm (i), where P atm is atmospheric pressure and g s and VPD are stomatal conductance and vapour pressure deficit, respectively, as described below . The other equations used to calculate micrometeorological covariates, soil humidity and photosynthetic active radiation, as well as those functions describing the water cycle (including soil evaporation and plant transpiration) are explained in the original model formulation by . Therefore, they will not be described here. The rest of the model was modified as follows.

Modelling the effect of climatic forcing on photosynthesis
Leaf photosynthesis (A n ) is calculated based on the biochemical model of Farquhar et al. (1980). A n is a function of the carboxylation (V c ), oxygenation (V o ) and leaf dark respiration rates (R d ): where photosynthesis during the day i is limited by either the rate of carboxylation when Rubisco is saturated (W c ) or when it is limited by R d was considered a fixed function of A c (0.006 × A c ) because this formulation performed better than an exponential function of temperature (Sala and Tenhunen, 1996;De Pury and Farquhar, 1997;Bernacchi et al., 2001). Following De Pury and Farquhar (1997): where C i is the CO 2 intercellular concentration, is the [CO 2 ] compensation point for photosynthesis in the absence of dark respiration, and K c and K o are the kinetic Michaelis-Menten constants for carboxylation and oxygenation, respectively. V cmax and J max are temperature-dependent parameters, as outlined below. Photosynthesis is known to respond to the carbon concentration within chloroplasts C c rather than to C i . Throughout the paper we retain the notation presented here in Eqs. (1) and (2) but discuss below how mesophyll conductance is taken into account empirically in relation to water stress when calculating g s and acknowledge the possible limitations of our approach (Reichstein et al., 2002;Grassi and Magnani, 2005;Flexas et al., 2006;Sun et al., 2014).
Climate influences leaf photosynthesis calculations through the temperature dependence of different parameters (Bernacchi et al., 2001;Nobel, 2009). , K c and K o were modelled using Arrhenius functions of daily mean temperature (T day , in • C) with parameters from De Pury and Farquhar (1997). We modelled J max as a fixed rate of V cmax (J max (i) = J coef · V cmax (i)) after comparing it with different temperature-dependent formulations (De Pury and Farquhar, Table 1. Characteristics of mean annual gross primary productivity, climatic (annual means) and growth data. Standard deviations are shown in parentheses. Precipitation: mean annual precipitation; T max : annual mean of mean daily maximum temperature; T min : annual mean of mean daily minimum temperature; length: chronology year replicated with more than 5 radii; RW: mean annual ring width; Rbs: mean correlation between series; AR: mean autocorrelation of raw series; MS: mean sensitivity; EPS: mean expressed population signal. Rbs, AR, MS and EPS are classical statistics to characterize growth chronologies, and they follow Fritts (1976 Maseyk et al., 2008). The model behaviour was better when the temperature dependence of V cmax was modelled using a logistic function (Gea-Izquierdo et al., 2010) rather than an exponential function as in : V max , V b and V ip are parameters to be estimated, with V max being the asymptote and V ip the inflection point. θ p is a soil water stress function dependent on the soil moisture conditions of the previous year. It takes into account the downregulation of photosynthesis in response to protracted drought through its impact on the photosynthetic capacity of active LAI in evergreen species caused by constraints in V cmax , in turn produced by irreversible photoinhibition, modifications in leaf stoichiometry and/or the aging of standing foliage through lower leaf replacement rates in response to long-term water stress (Sala and Tenhunen, 1996;Niinemets and Valladares, 2004;Niinemets, 2007;Vaz et al., 2010).
where p str is a parameter to be estimated and SWC 180 is the mean soil water content (mm) from July to December of the previous year. Photosynthesis is coupled to the calculation of stomatal conductance, which is estimated using a modified version of the Leuning (1995) equation: where g 1 and VPD 0 are parameters; VPD(i) is daily vapour pressure deficit; C s is the leaf surface [CO 2 ]; θ g is a nonlinear soil water stress function calculated as soil b and soil ip are parameters; and SWC(i) is daily soil water content (mm). θ g accounts for variability in gas exchange under drought conditions which cannot be taken into account through stomatal control alone; thus, the variability can also be related to, e.g., mesophyll conductance or stomatal patchiness. Therefore, with this empirical expression, we partly represent the effect of CO 2 fractionation during mesophyll conductance under water stress, acknowledging that this will likely be more complex under environmental stress (Reichstein et al., 2002;Grassi and Magnani, 2005;Flexas et al., 2006;Sun et al., 2014). The coupled photosynthesis-stomatal conductance system of equations was estimated separately for sun and shade leaves. Canopy photosynthesis was integrated using LAI, divided into its sunlit and shaded fractions (De Pury and Farquhar 1997). Transmission and absorption of irradiance was calculated following the Beer-Lambert law as a function of LAI, with LAI sun = (1 − exp(−LAI)) · K b (K b is the beam light extinction coefficient, which was set to 0.8) and LAI shade = LAI − LAI sun . In the mixed stand (Fontblanche), photosynthesis was calculated separately for Q. ilex and P. halepensis and then integrated to obtain stand estimates of forest productivity.

NSC
--C roots Figure 1. Outline of the different phenological phases (P1 to P5) and carbon allocation in the model within a given year. A n : net daily carbon assimilation; NSC: storage (non-structural carbohydrates); GDD: growing degree days; GDD l = parameter determining shift from P2 to P3 (see text); C: carbon allocated either to the stem, canopy or roots; d: day of year. Solid arrows correspond to allocation within the plant, whereas dashed arrows correspond to litterfall (canopy or roots). f 3 and f 4 are nonlinear functions of soil water content and temperature, determining carbon allocation to different compartments (see text for more details).

Modelling the effect of climatic forcing on carbon allocation
The model allocates daily carbon assimilated either to the canopy, stem, roots or storage of non-structural carbohydrates (NSC) to mimic intra-annual carbohydrate dynamics Dickman et al., 2015). Although trees can store carbon within different above-ground and belowground compartments (Millard et al., 2007), carbon storage is treated as a single pool within the model. Tree autotrophic respiration (R a , in addition to R d ) is modelled as a function f (i) of daily photosynthesis and maximum daily temperature (T max ; Sala and Tenhunen, 1996;Nobel, 2009): where p respi is a parameter. Net photosynthesis is calculated for day i as A N (i) = A n (i) − R a (i). This assumption means that respiration would be considered zero when there is no photosynthesis; hence, maintenance respiration would not be taken into account those days. Although this could bias the overall carbon balance, we assume that this effect will be very reduced in the studied forests because they present photosynthetic activity all year round (see results). The model simulates several phenological phases during the year (see Fig. 1): (P1) winter period where all photosynthates assimilated daily, A N (i), are allocated to the storage reservoir (NSCs) but there is no accumulation of growing degree days (GDD).
(P2) winter period where all A N (i) are allocated to storage (i.e. the same as in (P1)), but in contrast to (P1) there is active accumulation of GDD, which defines the threshold GDD 1 to trigger the next phenophase (P3) (budburst, leaf flush).
(P3) budburst, where carbon-available C T (i) = A N (i)+C bud (C bud is daily C storage utilized from buds, a parameter) is either allocated to the canopy, to roots or to the stem.
(P4) once the canopy has been completed in (P3), the next phenophase (P4) starts; in this period, daily photosynthates A N (i) are allocated either to the stem or to storage; (P5) the last phenophase (P5) starts when the photoperiod (parameter) crosses a minimum threshold in fall. In this phase, root mortality occurs. Otherwise (P5) is similar to (P1) and (P2) in the sense that all A N (i) is used for storage until next year's (P3) starts.
The allocation of carbon to different plant compartments is complex because it can be decoupled from photosynthetic production depending on different factors, some of them climatic, acting on different temporal scales (Friedlingstein et al., 1999;Sala et al., 2012;Chen et al., 2013;McMurtrie and Dewar, 2013). In this new version of the model, we set the different allocation relationships as nonlinear functions of temperature and soil water content, (P3) and (P4) following the functional relationships described in Gea-Izquierdo et al. (2013). This means that now we take into account homeostatic acclimation processes at the canopy level related to LAI dependence on water availability (Hoff and Rambal, 2003;Sala and Tenhunen, 1996;Reichstein et al., 2003). LAI is negatively related to long-term drought because litterfall is negatively linked to water stress Misson et al., 2011) and bud size depends on the climate influencing the period of bud formation (Montserrat-Marti et al., 2009). Therefore, the actual carbon that can be allocated to the canopy in (P3) of year j (AlloC canopy (j )) was set as a function of the previous year's moisture conditions (θ LAI (j )) and the maximum carbon that can be allocated to the canopy (MaxC canopy ). MaxC canopy is calculated from LAI max and SLA, and AlloC canopy (j ) = θ LAI (j ) × MaxC canopy , where p LAI is a parameter to be calibrated representing the threshold over which θ LAI (j ) = 1 and SWC 250 is mean soil water content for May-December of the previous year. Leaf turnover is variable within years and partly related to water availability (Limousin et al., , 2012. We considered a mean leaf turnover rate of 3 years for pines and et al. (2008) and Limousin et al. (2009; Fig. 1). C allocation to the canopy (i.e. including primary growth) in (P3) is calculated as C canopy root/leaf ; Ratio root/leaf was fixed to 1.5 for both species Ourcival, unpublished data), and where p 3moist , p 3temp and p 3sd are parameters representing the scale of the SWC and the optimum and dispersion of the T max functions, respectively. The carbon allocated to the stem with h 3_1 (i) as in Eq. (9); st 3moist , st 3temp and st 3sd_temp are parameters as in h 31 (i). The carbon allocated to roots in (P3) is set complementary to that of the other compartments to close the carbon budget within the tree, i.e. C roots ( Finally, in (P4) carbon-assimilated daily A N (i) is allocated either to stem growth or to storage until changing to (P5). In (P4), the amount of carbon to be allocated to stem growth is also set as a function of climatic forcing: where st 4temp and st 4sd_temp are parameters.

Eddy covariance CO 2 flux and dendrochronological data
The process-based model was calibrated using daily gross primary productivity (GPP), dendrochronological data and inventory data. To develop the model, those functions used to model daily stand photosynthesis (i.e. Eq. 1 to 9) were first calibrated against GPP values. GPP estimates were obtained from half-hourly net CO 2 flux measurements (NEP). GPP was obtained as the difference between measured net ecosystem productivity and calculated ecosystem respiration (Reichstein et al., 2005). Negative GPP values were corrected following Schaefer et al. (2012). Half-hourly GPP data were integrated to obtain daily estimates for the period 2001-2013 (Puechabon, methods detailed in Allard et al., 2008) and 2008-2012 (Fontblanche , Table 1). In the second step, those functions used to model how carbon assimilated and/or storage is allocated to the growth of the tree stem (i.e. Eqs. 10 and 11) were developed using calculated annual stem biomass increment time series. Stem biomass increment chronologies were built combining dendroecological data and forest inventory data collected at each site. We built one chronology for Q. ilex in Puechabon, a second for Q. ilex in Fontblanche and a third one for P. halepensis at Fontblanche (Fig. 2). For pines, two perpendicular cores were extracted using an increment borer from 25 trees in fall 2013, whereas for oaks we used cross sections. In Fontblanche, 15 oak stems were felled and basal sections collected in spring 2014. A total of 17 oak stems from Puechabon were logged in 2005 and 2008. The age and diameter distributions of the studied forests are depicted in Fig. A2.
All samples were processed using standard dendrochronological methods (Fritts, 1976). Annual growth (RW) was measured using a stereomicroscope and a moving table connected to a computer. RW cross-dating was visually and statistically verified. RW estimates were transformed to basal area increments (BAI, cm 2 yr −1 ). Mean BAI chronologies were obtained by averaging individual tree BAI time series. , 12, 3695-3712, 2015 www.biogeosciences.net/12/3695/2015/ In Fontblanche, BAI during the period 1987-1995 was standardized relative to the mean calculated after excluding that period (Fig. 2). BAI data were standardized because we did not find a climatic explanation for the abrupt growth peak observed in Fontblanche during that period (Fig. 2). Therefore, we assumed that it had been caused by a release event (i.e. reduction in competition) produced by the death of neighbours as a consequence of winter frost during 1985and 1987(Vennetier, personal communication, 2014. These two frosts were reflected by the presence of characteristic frost rings in most individuals from Fontblanche. To scale BAI chronologies to the same units as annual stem biomass (which is an output of the model), we used plot inventory data collected around the flux towers at the two sites. Inventory data included stem diameter for all trees and tree height collected for a subsample every 2 years during 2007-2011 in Fontblanche as well as annual diameter estimates for the period 1986-2011 for Puechabon. Individual annual biomass increments were estimated by subtracting the stem biomass of one year from that of the next; then, stand stem biomass increments (SBIs, g C m −2 yr −1 ) were calculated by integrating plot data. Stem biomass was calculated using allometric functions. For pines, we calculated stem biomass using diameter and estimated stem height assuming that the tree bole follows a paraboloid shape (Li et al., 2014). For oaks, stem biomass was calculated following Rambal et al. (2004). Once SBI had been estimated for the years when we had available inventory data, BAI chronologies were correlatively scaled to SBI units (g C m −2 yr −1 ). We built two mean stand SBI chronologies, one for each site, meaning that we analysed carbon allocation within stands, not differentiating between species in Fontblanche. These two SBI chronologies were used to calibrate sitewise Eqs. (10) and (11).

Model development and analyses
Parameters were selected according to the ecological characteristics of the species, exploring the model using comprehensive sensitivity analysis to sequentially optimize groups of parameters. In a first step, a group of common parameters was selected using GPP data from Fontblanche ( Table 2). The species-dependent parameters selected for Q. ilex in this first step (those parameters in Table 2 which are common to the two sites) were independently validated when applied in Puechabon. In a second step, a subset of site-dependent parameters was calibrated against GPP and SBI data. Four parameters from Eqs. (6) and (9) were calibrated using GPP data, and five parameters in Eqs. (10) and (11) were calibrated using stem biomass increment data ( Table 2). The local parameters were calibrated constrained to an ecologically realistic range and using a global optimization algorithm and maximum likelihood principles (Gaucherel et al., 2008).
To compare model output with stem biomass chronologies as estimated from dendroecological data, we used only the period for which we had available daily meteorological data , which was also a period that did not include juvenile years with increasing BAI (BAIs reached an asymptote after increasing for the first 15-20 juvenile years; Fig. 2). The model does not take into account how size differences in allometry or ontogeny affect carbon allocation (Chen et al., 2013). We tried to keep the model as simple as possible also because we had no such data to calibrate ontogenic effects. Hence, the model is designed for non-juvenile stands with canopies that reached a steady state with asymptotic LAI max . For the same reasons it does not take into account how changes in management affect carbon allocation. The model was analysed in terms of goodness of fit. Additionally, for the period for which we had available daily meteorological data, we simulated time series of GPP, ecosystem water use efficiency (WUE = GPP/ET, with ET being actual evapotranspiration) and intrinsic water use efficiency of sun leaves (iWUE = A N /g s ) calculated following Beer et al. (2009).

Results
The studied evergreen forests exhibit a bimodal pattern in GPP with maxima in spring and autumn (Fig. 3) as often observed in Mediterranean ecosystems (e.g. Baldocchi et al., 2010). GPP was above 0 almost every day of the year, including in winter, particularly at the milder site, Fontblanche (Table 1). This means that there is active photosynthesis all year round in these evergreen forests, including during both periods of climatic stress, i.e. those with low temperature and short photoperiod in winter and with low moisture availability in summer (Fig. 3). Mean annual GPP was 1431.4 ± 305.4 g C m −2 yr −1 and precipitation 642.7 ± 169.7 mm in Fontblanche, whereas it was 1207.3 ± 206.7 g C m −2 yr −1 and 1002.6 ± 328.2 mm in Puechabon (see Table 1 for more details). Mean GPP was higher at Fontblanche because carbon assimilation was greater in the low-temperature winter period but similar the rest of the year (Fig. 3). Stem growth did not show any longterm (decadal) growth trend for any of the species studied (Fig. 2).
The model accurately represented the low-frequency response of daily GPP: both the seasonal variability in GPP within years and variability in GPP among years (Fig. 4). The model explained over 50 % of the annual biomass growth variance, and 46 and 59 % of daily GPP in Fontblanche and Puechabon, respectively (Fig. 4). This means that we were able to mimic the daily, seasonal and long-term trends in stand productivity with unbiased estimates but also to model how carbon is allocated to stem growth throughout the year during the different phenophases described. The model assumed species-specific carbon allocation responses set to the different plant compartments as nonlinear functions of temperature and soil moisture. These relationships were biologically meaningful in the sense that photosynthesis and carbon www.biogeosciences.net/12/3695/2015/ Biogeosciences, 12, 3695-3712, 2015 Table 2. Model parameters. Those parameter differing between sites were optimized either with GPP data (photosynthesis and allocation module) or with growth-based biomass increment chronologies (allocation module). The rest were common parameters for both sites and selected while developing the model in the first step for Fontblanche using GPP data (represented in the "Cal" column by "-"   allocation could be decoupled to some extent as a function of climatic variability. Once the canopy had been formed in spring, the model allocated more carbon to the stem and less to storage when less severe climatic stress occurs, i.e. with decreasing temperatures and more humid conditions (Fig. 5).
Both sites exhibited an increase in temperature particularly evident in the maximum values, but only Puechabon suffered a decrease in annual precipitation between 1960 and 2012 (Fig. A1). In the model, the studied forests acclimated to changing conditions in the last decades, coupling different physiological traits, and simulated annual GPP largely followed the overall trends in precipitation observed. In Fontblanche, which is milder and receives less precipitation, GPP has remained stable since the 1960s and presented no apparent long-term trend (Fig. 6). In contrast, at the coldest and rainiest site (Puechabon), the model simulated a decrease in GPP (Fig. 6), which was driven by the prevailing decrease in precipitation observed since the 1970s (Figs. A1; 6). This reduction in GPP was partly a consequence of decreased LAI in response to enhanced long-term water stress  The modifier [0,1] is a function of soil water content (SWC) and maximum temperature (Tmax); multiplied with available daily carbon, it gives the distribution of daily carbon allocated to secondary growth and storage.
( Fig. A3; Limousin et al., 2009;Misson et al., 2011). Simulated long-term decadal trends in mean annual stomatal conductance were similar and decreased at the two sites with greater water stress as a consequence of enhanced temperatures (Fig. 6). The two species studied showed a long-term increase in simulated iWUE (Fig. 7) following the decrease in simulated g s (Fig. 6). The interannual variability in WUE and iWUE were highly and positively correlated (Fig. 7). However, in the long-term they followed a different pattern, particularly in Puechabon, where there was a recent decline in WUE (not observed in iWUE) forced by trends in ET and GPP (Fig. 7). This means that the recent reduction in simulated GPP was proportionally greater than that of simulated ET (Figs. 6; A3).

Linking photosynthetic production to carbon allocation as a function of climate
The model calculates stand productivity and carbon allocation to stem growth in response to climate and [CO 2 ] with realism. It is particularly well suited to mimic the effect of water stress in plant performance by the explicit assessment of different acclimation processes at the canopy level, including changes in stomatal conductance and photosynthetic capacity (Sala and Tenhunen 1996;Reichstein et al., 2003;Limousin et al., 2010;Misson et al., 2011). Additionally, the model simulates carbohydrate storage dynamically as a function of environmental variability. Climate can affect differently the carbon dynamics and pattern of C allocation to different tree compartments at different phenophases. In the model the storage reservoir is an active sink for assimilated carbon during some periods of the year and a source in spring to be used in primary and secondary growth (Fig. A5). Additionally stem growth is limited by climatic constraints (in (P3) and (P4)) rather than just by the amount of available carbohydrates (Millard et al., 2007). This means that water stress and optimum temperature directly affect the modelled processes, assuming that cell-wall expansion in the xylem can be related to climatic variability differently from photosynthetic production (Sala et al., 2012). The model showed C limitation (for primary growth) in the years when LAI max was not achieved (i.e. a limitation in LAI is driven by limitations in the C supply in spring), e.g. all years in Puechabon for the period shown in Fig. A5 (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) but only those years in Fontblanche when the minimum value considered as a threshold was reached. Therefore, both C-source (photosynthesis) and C-sink (just related to growth; other sinks such as volatile organic compounds or root exudates are not explicitly included in the model) limitations can be assessed in different years within one site and even at different periods within the same year (Millard et al., 2007;Sala et al., 2012;Chen et al., 2013;Fatichi et al., 2014). This hypothesis seems plausible as drought stress affects both C-source (e.g. through reduced stomatal conductance) and C-sink limitations (e.g. cell water turgor, hydraulic performance; McDowell et al., 2013). Whether the pattern of C storage simulated is realistic is something that needs to be validated against actual data. However, the flexible way in which stored C is modelled has much potential to improve ecosystem models that only view a carbon-source limitation (Sala et al., 2012;Friend et al., 2014). Water stress is generally considered the greatest limitation for Mediterranean ecosystems, driving a close relationship between precipitation and both growth and photosynthesis (Breda et al., 2006;Pereira et al., 2007;Baldocchi et al., 2010;. Our results show that a long-term decrease in precipitation triggered a decrease in simulated GPP at the rainier and coldest site. However, this decline was not expressed in the growth trends. This means that long-term productivity and the allocation of C to secondary growth were decoupled and did not match (Sala et al., 2012;Chen et al., 2013;Fatichi et al., 2014). The existence of trade-offs between carbon assimilation and allocation in relation to environmental variability suggests exercising caution when using growth as a direct proxy to investigate stand productivity dynamics (e.g. Piovesan et al., 2008;Peñuelas et al., 2008;. GPP was greater at the site receiving less precipitation, which could be related to differences in soil retention capacity. However, both soils are calcareous, shallow and stony, and differences in GPP were, to a large part, explained by less limitation for carbon assimilation of low winter temperatures at the warmest site (Fontblanche). They can also be a result of a different species composition (oak vs. pine oak). LAI is greater at the site yielding higher annual GPP. Nonetheless, had this factor been responsible for the observed differences in winter photosynthesis, there would also have been differences in spring photosynthesis, which was not the case (Fig. 3). A better understanding of the underlying processes determining carbon allocation will benefit process-based models (Sala et al., 2012;Fatichi et al., 2014). Model parameters were within the range found in the literature, bearing in mind that using a daily time step to study differential processes or not distinguishing between leaf ages will affect the scaling of parameters such as J max , V cmax or R d (De Pury and Farquhar, 1997;Grassi and Magnani, 2005;Maseyk et Vaz et al., 2010). Daily climatic data are readily available on a greater spatial scale than data with a higher temporal resolution, which increases the applicability of daily models. Model performance could be improved by addressing respiration changes related to ontogeny, allometry and nutrient limitations (e.g. N/P) on photosynthesis or by including more complex upscaling of leaf-level photosynthesis (Niinemets et al., 1999;Niinemets, 2007;Chen et al., 2013;McMurtrie and Dewar, 2013). However, it is difficult to find suitable data to calibrate such processes. Similarly, it would be challenging to include allocation to reproductive effort in the carbon budget. This is because, even if it is influenced by water stress in the studied forests (Pérez-Ramos et al., 2010), there is still great uncertainty regarding the causal factors driving multi-annual variability in fruit production (Koenig and Knops, 2000). Addressing stand dynamics would also help to generalize model applicability. Stand disturbances modifying stand competition can leave an imprint on growth for more than a decade whereas they do not seem to affect stand GPP over more than 1 or 2 years if the disturbance is moderate (Misson et al., 2005;Granier et al., 2008). In response to changes in competition, the trees modify the carbon allocation pattern or keep the root : shoot ratio constant to enhance productivity on a per-tree basis but up to an asymptotic stand GPP. Still, the model behaviour was good compared with other studies that addressed ontogenic changes in the carbon allocation response to photosynthesis (Li et al., 2014) and similar or better than that of other mechanistic approaches calibrated to standardized dendroecological data Evans et al., 2006;Gaucherel et al., 2008;Tolwinski-Ward et al., 2011;Touchan et al., 2012).

Forest performance in response to recent climate change and [CO 2 ] enhancement
Few studies under natural conditions have observed a net increase in growth rates in response to enhanced [CO 2 ] levels since the late 1800s, meaning that other factors, such as water stress and/or N or P, were more limiting for photosynthesis and/or allocation to growth than [CO 2 ] (Niinemets et al., 1999;Peñuelas et al., 2011;McMurtrie and Dewar, 2013;Lévesque et al., 2014). Yet the forests have increased their iWUE. This can be partly a passive consequence of enhanced [CO 2 ], but higher iWUE observed in more water-stressed sites suggests that climate is co-responsible for an active acclimation of physiological plant processes (Keenan et al., 2013;Leonardi et al., 2012;Saurer et al., 2014). These processes would include a higher stomatal control, like in our results, where, in turn, we did not observe any increase in longterm carbon assimilation. The mean annual stomatal conductance simulated was driven by climate but also decreased simultaneously in time with increasing [CO 2 ] (Fig. A4). Furthermore, there is debate on whether there has been an increase in ecosystem WUE in response to recent changes in [CO 2 ] under a warming climate (Beer et al., 2009;Reich-stein et al., 2002;Keenan et al., 2013). In our results the high frequency of WUE followed that of iWUE, but there was some mismatch between the two traits in the low frequency. We observed no dominant time trends in simulated annual WUE for the period 1980-2000 at the site where precipitation remained stable, whereas a decrease in WUE following that in GPP was particularly evident at the site experiencing a drier climate in recent years. This trend was not observed in iWUE, which means that reductions in GPP and g s were proportionally greater than those in ET (Figs. 6, 7, A3).
Higher [CO 2 ] concentrations enhance photosynthesis with the equations used to calculate leaf photosynthesis in biochemical models (e.g. Gaucherel et al., 2008;Keenan et al., 2011;Leonardi et al., 2013;Boucher et al., 2014). Thus, the absence of a long-term increase in GPP and growth would not mean that enhanced [CO 2 ] was not beneficial for model outputs (particularly in the case of C-source limitation) but that the net control exerted by other factors such as climatically driven stress was more limiting than that of [CO 2 ] availability: growth and photosynthesis would have been lower had we used constant [CO 2 ] with the same model parameters. The absence of any modification in the growth trends, even if there are changes in WUE, would express a sink limitation mostly related to hydraulic constraints (Peñuelas et al., 2011;Sala et al., 2012;Keenan et al., 2013). Often, the trees show a growth decline at those sites where there is an enhancement in long-term water stress that dominates species performance (e.g. Bigler et al., 2006;Piovesan et al., 2008;. In contrast, it has been observed that, under certain conditions, trees have increased growth with warming since the 1850s (Salzer et al., 2009;. These studies suggest the existence of a positive effect of warming, rather than of CO 2 fertilization, upon growth in forests where water stress is not the most limiting factor. Our study sites are located at the northern limit of the Mediterranean region, meaning that the two species studied occupy drier and warmer areas further to the south. The two species have different functional characteristics, e.g. oaks are anisohydric, whereas pines tend to be isohydric. This confers different capacities of adaptation to climate change on them, which means that they should play different roles in future stand dynamics. Our results show the existence of trade-offs in response to climate at different phenological periods. This is important since synergistic environmental stresses acting at different periods can trigger tree mortality (McDowell et al., 2013;Voltas et al., 2013). Model sensitivity analysis could be performed to discuss the influence of specific factors, such as climate or [CO 2 ], causing instability in the climate-growth response (D'Arrigo et al., 2008;Boucher et al., 2014). However, [CO 2 ] enhancement and climate warming are mixed in analysis performed using data from field studies, which can make the isolation of their effect problematic. The model can be applied using abundant dendrochronological data used to determine the site-dependent parameters. This would provide much flexi-bility for investigating growth trends and forest performance in response to global change on a larger scale.

Conclusions
By developing an original process-based model with carbon allocation relationships explicitly expressed as functions of climate, we accurately simulated gross primary productivity and the allocation of carbon to secondary growth in evergreen Mediterranean forests. Different processes were modelled as functions of environmental variability, including [CO 2 ] and climate. The studied forests showed trade-offs in carbon allocation to different plant compartments in response to stress in different seasons: with low temperatures and a short photoperiod in winter and with moisture shortage in summer. We modelled a decreasing time trend in stomatal conductance, which would suggest a partly active increase in iWUE in the forests studied. Interannual variability in WUE followed closely that in iWUE. However, WUE exhibited a decreasing trend at the site where we simulated a decrease in LAI and GPP in response to a decrease in annual precipitation since the 1970s. Long-term GPP has remained at similar levels in the last 50 years in just one stand, whereas it declined in the forest suffering a reduction in precipitation. This suggests different acclimation processes at the canopy level and in the pattern of allocation in response to enhanced xericity and increasing [CO 2 ] levels; these acclimation processes could not counterbalance the negative effect of warming only at the site where there was a simultaneous decrease in precipitation. Tree growth was partly decoupled from stand productivity, highlighting that it can be risky to use growth as a direct proxy for GPP. The model is flexible enough to assess both C-source and C-sink limitations and includes a dynamic estimation of stored C. These features would improve ecosystem models with a fixed C-source formulation. By calibrating a limited number of parameters related to carbon allocation, there is great potential for using the model with abundant dendroecological data to characterize past instability in the growth response in relation to environmental variability and to simulate future forest responses under different climatic scenarios.