Ecophysiological modeling of photosynthesis and carbon allocation to the tree stem in the boreal forest

A better understanding of the coupling between photosynthesis and carbon allocation in the boreal forest, with implicated environmental factors and mechanistic rules, is crucial to accurately predict boreal forest carbon stocks and fluxes, which are significant components of the global carbon budget. Here we adapted the MAIDEN ecophysiological forest model to better consider important processes for boreal tree species, such as non-linear acclimation of photosynthesis to temperature changes, canopy development as a function of previous year climate variables influencing bud formation, and 15 temperature dependence of carbon partition in summer. We tested these modifications in the eastern Canadian taiga using black spruce (Picea mariana (Mill.) B.S.P.) gross primary production and ring-width data. MAIDEN explains 90% of the observed daily gross primary production variability, 73% of the annual ring width variability and 20-30% of its high frequency component (i.e. when decadal trends are removed). The positive effect on stem growth due to climate warming in the last decades is well captured by the model. In addition, we illustrate the improvement achieved with each introduced 20 model adaptation and compare the model results with those of linear response functions. This shows that MAIDEN simulates robust relationships with the most important climate variables (those detected by classical response-function analysis), and is a powerful tool for understanding how environmental factors interact with black spruce ecophysiology to influence presentday and future boreal forest carbon fluxes.


Introduction
Photosynthetic production is the primary factor affecting growth of trees and other vegetation.However, empirical studies have shown that the correlation between photosynthetic production and the diameter growth of trees is far from perfect (Gea-Izquierdo et al., 2014;Rocha et al., 2006;Berninger et al., 2004).This imperfect correlation is due to the fact that plant hydraulics (e.g., turgor pressure) and thermal limitations during very short periods of time can be more important than carbon (C) availability for secondary tree growth (Kirdyanov et al., 2003;Rossi et al., 2016;Zweifel et al., 2016;Fatichi et al., 2014; secondary growth is the increase in the girth of the plant roots and stems).These factors influence the proportion of net primary productivity allocated to stem growth each year, dampening the correlation between gross primary production (GPP) and growth.A better understanding of these factors and of carbon allocation mechanisms is needed when studying forest dynamics, forest carbon balance and the impact of climate change on forests.Carbon allocated in different tree components (e.g., canopy, stem or roots) has a specific function and is stored for a different length of time (Moorcroft, 2006).

F. Gennaretti et al.: Ecophysiological modeling of photosynthesis and carbon allocation
The varying roles of allocation and photosynthetic production are integrated in ecophysiological models (Li et al., 2014).Such models are important tools for analyzing the direct influence of climate and other environmental factors (e.g., CO 2 concentration) on tree growth and biogeochemical processes in forest ecosystems (Li et al., 2016).Climategrowth relationships have traditionally been assessed using empirical response functions based on linear relationships, thus considering the underlying processes as a black box.In contrast, ecophysiological models are built on mechanistic rules and allow for consideration of non-stationarity and nonlinearity in tree responses to environmental variables as well as their interactions (Vaganov et al., 2006).Ecophysiological models may be refined using model-data fusion approaches and optimization techniques (Guiot et al., 2014).
Different models with a different degree of ecophysiological complexity and/or spatiotemporal resolution have already been used to investigate the influence of climate and weather on tree growth in the boreal forest.Some studies focused on the drivers of photosynthetic capacity.For example, Mäkelä et al. (2004) proposed a model to study the influence of temperature on the seasonal variation in photosynthetic production of Scots pine through a delayed dynamic response.Other studies focused on the drivers of carbon allocation.For example, in Manitoba, Canada, a model related GPP and carbon allocation to absorbed photosynthetically active radiation as a function of environmental constraints (Girardin et al., 2008).Another model, called CASSIA (Schiestl-Aalto et al., 2015), was developed to investigate how environmental factors and the ontogenetic stage of tree development influence the annual course of carbon sink-source dynamics in Scots pine stands.However, despite recent progress, few models have been able to simultaneously simulate the meteorological control on daily photosynthetic production and the meteorological and phenological controls on daily carbon allocation for temperature-limited boreal forest ecosystems.Such models should be able to simulate the following observed phenomena: (i) the delayed response of photosynthesis to temperature (Gea-Izquierdo et al., 2010;Mäkelä et al., 2004), (ii) the influence of preceding season conditions on current-year canopy development (Salminen and Jalkanen, 2005) and (iii) a strong positive relationship between wood biomass production and temperature (Cuny et al., 2015).
Here, we try to fill this gap by adapting the MAIDEN forest ecophysiological model, developed for temperate and Mediterranean environments (Misson, 2004;Gea-Izquierdo et al., 2015), to mimic how weather and climate influence photosynthesis, phenology and carbon allocation in the North American boreal forest on a daily basis.MAIDEN offers an ideal framework to analyze the impact of introducing relevant processes for carbon assimilation and allocation in temperature-sensitive boreal trees into the model.Indeed, the model simultaneously simulates the course of photosynthesis and sets different phenological phases to determine the allocation of carbon to different plant compart-ments in a dynamic manner.In this study, we first test and optimize new model features on GPP and growth data from black spruce (Picea mariana (Mill.)B.S.P.), the dominant tree species across the North American boreal biome.Second, we show the impact of single processes in the model runs and the improvements achieved with the new model adaptations.Last, we compare the simulated GPP and stem growth results with those obtained with conventional empirical linear response functions.This comparison allows us to verify that the process-based ecophysiological model satisfactorily reproduces the variability in the observed data and that its simulations keep robust relationships with the most significant climate variables.
2 Materials and methods 2.1 The MAIDEN model MAIDEN (Misson, 2004;Gea-Izquierdo et al., 2015) can consider the influence of several environmental factors on forest water and carbon cycles.Starting from daily minimum-maximum air temperature, precipitation and CO 2 atmospheric concentration (these are the minimum required input variables; radiation, relative humidity and wind speed are included when additional meteorological data are available; Misson, 2004), MAIDEN models the phenological and meteorological controls on GPP and carbon allocation (Fig. 1; see also flowcharts in Misson, 2004, andGea-Izquierdo et al., 2015).The model explicitly allocates carbon to different pools (storage, canopy, roots and stem) on a daily basis using phenology-dependent mechanistic rules.The model has already been successfully optimized for Quercus petraea (Matt.)Liebl.and 12 Mediterranean species, including several Pinus spp.and Quercus spp.(Gaucherel et al., 2008a, b;Danis et al., 2012;Misson, 2004;Misson et al., 2004;Boucher et al., 2014;Gea-Izquierdo et al., 2015, 2017).Thus far, the model has never been used to simulate forest growth under boreal conditions.
MAIDEN requires the definition of species-and sitedependent parameters (Misson, 2004;Gea-Izquierdo et al., 2015), such as soil texture and depth and the root-to-leaf mass fraction in the studied trees.The parameters that could not be set for the studied black spruce sites were analyzed with sensitivity analysis, and the most influential of them were estimated with Bayesian optimization algorithms (Robert, 1996) using observed time series (daily GPP and annual ring width) as a reference.In total, 6 parameters influencing the GPP for black spruce and 12 parameters controlling the carbon allocation to the stem (Dstem) were optimized (they are described in the following paragraphs and in Table 1).The optimization was based on Markov chain Monte Carlo (MCMC) sampling, which, through its iterations, only retains combinations of parameters satisfying some conditions (Supplement S1; Fig. S1  .MAIDEN-simulated phenology (blue), water (black) and carbon (red) fluxes.AN: net photosynthesis corresponding to net primary production.Cstored, Cstem, Ccanopy and Croots: carbon allocated daily to stored non-structural carbohydrates, stem, canopy and roots, respectively.DOY: day of the year (1-365).GDD: growing degree days.f 3 and f 4 : functions determining carbon allocation in phases 3 and 4. Cbud: amount of storage carbon that is used each day by the plant in phase 3. ment).Among the retained blocks of parameters, one block of six parameters controlling GPP ("plausible block GPP") and one block of 12 parameters controlling Dstem ("plausible block stem") were selected to illustrate the results with likely parameter values (Supplement S1).The robustness of the parameters' posterior distributions was tested on a crossvalidation exercise (Supplement S1).

Modeling the GPP of boreal forests
In MAIDEN, the daily stand GPP (g C m −2 day −1 ) is derived from the modeling of the coupled photosynthesis-stomatal conductance system.Leaf photosynthesis is calculated following De Pury and Farquhar (1997), while stomatal conductance is estimated using a modified version of the Leuning equation (Leuning, 1995;Gea-Izquierdo et al., 2015).The photosynthesis-stomatal conductance system is estimated separately for sun and shade leaves based on the photosyn-thetic photon flux density they receive.The partition of leaf area index (LAI) in its shaded and sunlit fractions and the transmission and absorption of photosynthetically active radiation (PAR) are computed as explained by Misson (2004), following De Pury and Farquhar (1997).After a sensitivity analysis, and as stated in the literature for boreal forests (Gea-Izquierdo et al., 2010;Mäkelä et al., 2004Mäkelä et al., , 1996)), we found that the modeling of assimilation/photosynthesis for black spruce is very sensitive to the parameters controlling the temperature dependence of the maximum carboxylation rate (Vcmax; µmol C m −2 of leaves s −1 ), the water stress level (θg) influencing the stomatal conductance and consequently the intercellular CO 2 concentration.The computations of Vcmax and θg used here are identical to those of the prior formulation of MAIDEN (Gea-Izquierdo et al., 2015).The Vcmax is modeled as follows: θg is a logistic function that varies from 0 (maximum stress) to 1 (no stress) at day i depending on the soil water content (SWC; mm).soilb and soilip are the slope and the inflection point of θ g i , respectively.With its already published MAIDEN configuration (Gea-Izquierdo et al., 2015), the model overestimated black spruce GPP in spring.This overestimation is due to the fact that the model has been developed for temperate and Mediterranean trees where no time delay between the recovery of photosynthesis and temperature increase in spring (i.e., no temperature acclimation) can be assumed.However, such a delay is common in boreal trees (Gea-Izquierdo et al., 2010;Mäkelä et al., 2004).For this reason, we modified MAIDEN by including an extra function and an extra parameter (τ ) to take into account the acclimation of photosynthesis to temperature.We replaced Tday in Eq. ( 1) by a temperature transformation (S), which responds smoothly with a determined time lag to temperature variations.S of day i was computed from the following differential equation (Mäkelä et al., 2004), which was solved with Euler's method: The new parameter τ is a time constant interpretable as the number of days needed by the photosynthetic apparatus to acclimate to changing temperature.
2.1.2Modeling carbon allocation to the stem (Dstem) in boreal forests MAIDEN allocates the daily available carbon from photosynthesis and stored non-structural carbohydrates to all plant compartments (stem, roots, canopy and storage) using functional rules specific to each of the five phenological phases characterizing a year (see Fig. 1).Although we maintained the original MAIDEN structure, we modified some previously used functional rules from Gea-Izquierdo et al. (2015) to consider significant processes for the boreal forest.Below, we describe the functional rules controlling Dstem according to phenological phases.
During "winter period 1" (phase 1), few processes are active.However, at the beginning of each year, the model defines the maximum amount of carbon that the canopy can potentially contain that year (AlloCcanopy j ; g C m −2 of stand) as a function of previous-year climate variables.Based on previous studies on black spruce forests (Girardin et al., 2016;Ols et al., 2016;Mamet and Kershaw, 2011), we modified the model to consider the effect of the previous-year April precipitation and July-August temperature that likely influence the length and the thermal-hydraulic stress of the previous growing season, respectively.Previous-year climate conditions of specific months are known to influence the shoot extension of boreal trees likely because they control the accumulation of resources in the buds (Salminen and Jalkanen, 2005).Here, we calculated the carbon potentially allocated each year to the canopy with the following equations: where Temp j −1 is the previous-year mean July-August temperature (detrended and transformed to z scores), Precip j −1 is the previous-year April precipitation (detrended and transformed to z scores) and MaxCcanopy is the absolute maximum canopy carbon reservoir determined based on forest traits, diameter distributions and previously published allometric equations (Chen, 1996;Bond-Lamberty et al., 2002a, b).CanopyT and CanopyP are two parameters that were optimized and represent the slopes of the relationships between CanopyMult (i.e., the overall climate dependence) and Temp j −1 or Precip j −1 , respectively.In this way, AlloCcanopy j may vary between 70 and 100 % of Max-Ccanopy, as in the previous version of the model (Gea-Izquierdo et al., 2015).
During "winter period 2" (phase 2), growing degree days (GDD) start to accumulate.We computed the accumulation of GDD by summing the mean daily temperature values over 3 • C (Nitschke and Innes, 2008;Man and Lu, 2010).MAIDEN simulates budburst (i.e., the transition from phenological phase 2 to 3) either when the GDD sum threshold is reached (parameter GDD1) or when a selected day of the year related to photoperiod is passed (parameter veg-phase23).With this model configuration, the start of the growing season overreacted to GDD yearly variations.To correct for this simulated bias, we modified MAIDEN by adding a mechanism reducing the interannual variability of budburst dates.This mechanism simulates the acclimation of the plants to varying GDD sums from year to year.The yearly time series of days of the year corresponding to budburst (determined by GDD and photoperiod) is smoothed at the beginning of each simulation with an n-year cubic smoothing spline.The integer number n was called day23_flex and optimized similar to the other parameters.
The "budburst phase" (phase 3) starts with budburst and ends when AlloCcanopy j is reached or when the carbon in the storage reservoir (i.e., stored non-structural carbohydrates) is lower than a minimum value (Misson, 2004).This phase was set to be shorter than 51 days based on available spruce budburst and shoot elongation data (Lemieux, 2010).During this phase, the daily available carbon (CT i ) comes from photosynthesis and the mobilization of storage carbon.The parameter Cbud, which was optimized, is the amount of storage carbon that is used each day by the plant.The total CT i amount is then allocated to the canopy, the roots or the stem following some functional rules.In the previous version of MAIDEN (Gea-Izquierdo et al., 2015), these rules were functions of daily soil moisture and air temperature.In this case, these rules did not improve the simulated results, and we retained a simpler version independent of the climate: where Cstem i is the portion of CT i allocated to stem and h 3 is a parameter to be defined in the range between 0 and 1.The rest of CT i is allocated to the canopy or the roots, with a prescribed 1.65 root-to-canopy mass ratio for black spruce (Czapowskyj et al., 1985;Jenkins et al., 2003).
During the "growth and accumulation phase in summer" (phase 4), CT i comes from only photosynthesis and is allocated either to stem growth or to storage as a function of climate forcing.For water-limited sites in the previous version of MAIDEN (Gea-Izquierdo et al., 2015), the allocation rule used a combination of daily soil moisture and air temperature as predictors.In this case, for temperature-limited sites, we used only temperature and set soil moisture at null (i.e., always equal to 1; note that for more water-limited boreal sites, this water stress dependence can be used): where Tmax i is the daily maximum temperature and st4temp is a parameter that corresponds to the inflection point of the function.The value 0.8 was chosen to force a minimum threshold of C allocation to the stem in this phase (at least 20 %) and to guarantee the correspondence between the inflection point and the temperature where approximately 50 % of CT i is allocated to the stem.
The transition from phase 4 to the "fall phase" (phase 5) is determined by either the parameter photoper (threshold of duration of daylight in hours) or by the occurrence of negative minimum daily temperature values after 1 September.During the "fall phase", all photosynthetic products are allocated to the storage reservoir, and mortality of fine roots occurs.No specific functional rule influences Dstem during this phase.
The equation controlling partial carbon losses from the canopy (i.e., litterfall), which influences the photosynthetic capacity through modifications of the total leaf area in the studied evergreen species, runs all year round.This equation is adapted from Maseyk et al. (2008): where outCcanopy i is the carbon loss from the canopy at day i and is influenced by parameters PercentFall, OutMax and OutLength (to be optimized), which determine the yearly canopy turnover rate, the day of the year with maximum losses and the length of the period with losses, respectively.

Model evaluation
The proportion of the observed variability explained by MAIDEN was evaluated with the coefficient of determination (R 2 ), which compares the performance of simulated time series relative to that of straight horizontal lines centered on the data:  S2) to optimize the six parameters influencing stand GPP simulated by MAIDEN for the studied species.

Ring width data from the northern Quebec taiga
We assumed that the yearly Dstem is proportional to tree-ring growth to use ring width data to optimize MAIDEN (12 influential parameters).A regional chronology (RW) and a detrended regional chronology (RWhighF) were obtained from 46 black spruce trees sampled in the riparian forests of five lakes in the eastern Canadian taiga (Gennaretti et al., 2014; the coordinates of the central point are 54.26 • N, 71.34 • W; see Fig. S3, data set S1 and Supplement S2).RWhighF was used as a reference for the optimization of the MAIDEN parameters, while the observed and simulated low frequencies were compared after the optimization of the model parameters.MAIDEN outputs were simulated for the central point of the source area of ring width data over the 1950-2010 period.

Climate data
MAIDEN needs daily climate data as inputs.These data were obtained from the gridded interpolated Canadian database of daily minimum-maximum temperature and precipitation for 1950-2015 (Hutchinson et al., 2009;   model.We used linear response functions to analyze the relationships between observed daily GPP at EOBS and the daily mean, maximum and minimum temperatures or weekly precipitation (explored time lag from 0 to 30 days before; in the case of precipitation, lag n indicated the sum of the daily precipitation of the week ending in day n).In this analysis, we excluded the winter days (days of the year between 15 November and 1 April) where GPP is zero.The 10 predictors most strongly correlated with GPP (and not highly correlated with each other; pairwise r ∈ [−0.8, 0.8]) were retained for the analysis.All linear response functions, resulting from a combination of these 10 predictors, were tested and classed according to their Bayesian information criterion (BIC).
We also used linear response functions to analyze the relationships between RWhighF and climate variables (same methodology as for GPP).We tested all monthly temperature and precipitation values of the previous and current years as predictors.Time windows of 31 days were used to obtain the time series of monthly data (over the 1950-2010 period) for each day (central day), averaging the values of each window and each year.These climate time series were also detrended.

GPP and tree-ring growth variability explained by MAIDEN
The optimized model (see parameters' posterior distributions in Figs.S4 and S5) explained a large proportion of the observed GPP daily variability (90 %; r = 0.95, df = 2918, p < 0.001; Fig. 2a).Although the model was optimized with daily data, the GPP time series also reproduced the annual variability of the observed data quite well (Fig. 3).
As expected, the ring growth variability at our sites was more linked to temperature than to precipitation variables (see Fig. 4a and Gennaretti et al., 2014;Mamet and Kershaw, 2011;Nicault et al., 2014).The model reproduced this correlation pattern (Fig. 4b) and explained approximately 20-30 % of the observed yearly RWhighF variability, corresponding to correlations of 0.58-0.66(df = 59, p < 0.001; Figs.2b and S6).This result is good because the simulated detrended annual GPP values (i.e., photosynthetic assimilation before any carbon allocation) had only a negative R 2 with RWhighF (Fig. 2c; meaning performance was worse than a straight line centered on RWhighF) and much lower correlations (Figs.S6 and S7).The variance explained by the model increased when the time series of stem growth were analyzed with their trends (R 2 = 0.73 and r = 0.86, df = 59, p < 0.001; Fig. 5b).The positive trend in response to the warming of the last few decades was well captured by the model simulations of stem increments, which included some CO 2 fertilization contribution (Fig. S8).

Mechanistic and regression-based diagnostics
The modeled impact of temperature on the maximum rate of Rubisco-catalyzed carboxylation (Vcmax) is shown in Fig. S9.This figure was obtained using Eqs.( 1) and (3) with the parameters of plausible block GPP and using actual temperature data.The obtained Vcmax values (up to 30 µmol C m −2 of leaves s −1 ; Fig. S9c) were comparable to those obtained for another mature black spruce forest in Saskatchewan, Canada (Rayment et al., 2002).Furthermore, the impact of soil water content on the water stress level (θg) influencing the stomatal conductance is shown in Fig. S10.Simulated GPP values were sensitive to all single parameters controlling Vcmax or θg, except soilb (Fig. S11-S15).The temperature transformation (S) introduced here in MAIDEN also influenced the simulation results (Fig. 6).With no time delay between photosynthesis and temperature increases (i.e., τ = 1 and S = Tday), MAIDEN overreacted to temperature variations in spring, and the GPP annual cycle was antedated (the start in spring and highest summer values were too early).In contrast, the use of S with τ values between 10 and 15 days synchronized the   GPP annual cycle with observations.This result means that black spruce photosynthetic capacity needs approximately 10-15 days to acclimate to a higher daily temperature (e.g., τ equal to 12.43 days was selected for plausible block GPP).This time delay is a little longer than that previously found for black spruce but comparable to values found for Scots pine (Mäkelä et al., 2004;Gea-Izquierdo et al., 2010, 2014).
We modified important processes for carbon allocation to adapt MAIDEN to black spruce.For example, previous-year precipitation and temperature values influenced the potential maximum amount of carbon that the canopy could contain during the growing season, as illustrated in Fig. 7a (see Eq. 4).If both previous April precipitation and July-August temperature indexes are negative, then the potential amount of carbon simulated by the model would be maximum, otherwise it would be minimum.This result was consistent with the correlations shown in Fig. 4.
Another important process is the start of the growing season.According to our simulations, the start could not occur later than 17 June (Figs.S5d and S16; Table 1) and was influenced by the GDD sum and the photoperiod, which are known to be relevant for black spruce budburst along with the tree provenance (Rossi and Bousquet, 2014).However, because we added a mechanism to smooth yearly variations (see the day23_flex parameter), more years were needed by the plants to acclimate to variable GDD accumulations in winter-spring.With the selected parameters to simulate stem growth, the median onset of the growing season was 10 June (similar to observations for black spruce in northern Mani- toba, Canada; Bronson et al., 2009), with a standard deviation of 7.8 days.If the smoothing term was excluded, then the standard deviation increased to 9.4 days (see Fig. S16a).
The inclusion of the smoothed mechanism also decreased the correlation between the simulated detrended annual Dstem and May average temperature from 0.70 to 0.59 (df = 58, p < 0.001).Although this correlation is still high, it was closer to the correlation between RWhighF and May temperature (r = 0.27, df = 58, p < 0.05; Fig. 4).These results show how the new model configuration decreased the yearly variability of the growth onset and helped achieve more plausible correlations with climate variables.According to the simulations, the onset of the growing season shifted by 7 days from 14 June to 7 June between the 1950-1970 and 1990-2010 periods (Fig. S16b-c).This result is consistent with the study of Bronson et al. (2009) on the effect of warming on black spruce budburst.During phase 3, which corresponds to budburst, a portion of the available carbon simulated by MAIDEN comes from stored non-structural carbohydrates that are from the current and previous years (parameter Cbud; see Table 1).In our case, Cbud was quantified as approximately 1.69 g C m −2 day −1 (Fig. S5f), and this remobilization improves the correlations between Dstem and RWhighF (Fig. S17).During phase 3 of our simulations, almost all available carbon was allocated to the canopy and roots (h 3 ≈ 0.9905; Eq. 5; Fig. S5g; Table 1).For this reason, the previously used soil moisture and temperature dependences, determining the portion of carbon allocated to the stem in Mediterranean evergreen woodlands (Gea-Izquierdo et al., 2015), did not improve the results and could be excluded here.The partition of carbon during the growth and accumulation phase in summer (phase 4) was instead modeled as a function of temperature (Eq.6).The simulations were highly sensitive to the st4temp parameter (Fig. 8c-d), and warmer temperatures corresponded to a greater amount of carbon allocated to the stem and less to non-structural carbohydrates (Fig. 8a-b).These results are in line with those of Cuny et al. (2015), who showed that woody biomass production is low in the first part of the growing season for most coniferous tree species because production follows the seasonal course of temperature (highest peak in summer).The simulated accumulation of carbon to the stem ended each year when the photoperiod became shorter than approximately 13.41 h (Fig. S5i; Table 1), corresponding to 2 September.The model performance was very sensitive to this parameter, which is known to impact black spruce dormancy induction (D'aoust and Cameron, 1982).
Another important process for carbon allocation is the loss of carbon from the canopy, a process that influences the seasonal course of the photosynthetic capacity.According to the simulations, the canopy mean annual turnover rate was approximately 13-14 % (Fig. S5j; Table 1), which corresponds well to previously published values for boreal spruce species ( Ťupek et al., 2015).The simulated annual cycle of canopy losses (Fig. S18) culminated on 2 July, and 80 % of litterfall occurred between 27 May and 19 July.This cycle is also similar to published results showing that the majority of litterfall (≈ 80 %) occurs in summer during needle growth for conifer species (Maseyk et al., 2008).
The comparison between MAIDEN simulations and classic linear response functions confirmed the quality and www.biogeosciences.net/14/4851/2017/Biogeosciences, 14, 4851-4866, 2017  ) at the Quebec Eastern Old Black Spruce site.Only the τ parameter determining the S values was allowed to vary, while the other parameters were fixed to the values of plausible block GPP.(a) τ is 1 day (S same as Tday).(b) τ is 6.71 days (a middle value between 1 and 12.43).(c) τ is 12.43 days (same τ than in plausible block GPP).(d) τ is 18.14 days (a higher value than in plausible block GPP).The R 2 between observations and simulations is reported in each plot.Panels (e-h) show the impact of the respective τ values on S if the daily Tday time series corresponds to a single step of 10 • C lasting 30 days.plausibility of the simulated results with the process-based model.MAIDEN performed better than response functions in explaining the variability of the daily GPP (R 2 = 0.90 vs. 0.69; Table 2).In the case of annual radial growth, the explained variability with the best response function (50 %; Table 3) was greater than with MAIDEN (20-30 %; r ≈ 0.65, df = 59, p < 0.001; Fig. 2b).However, the MAIDEN-simulated time series maintained the relationship with the significant monthly climate variables detected in the response function analysis.Correlation coefficients of −0.39, 0.46 and 0.57 (df = 58, p < 0.01) were obtained between MAIDEN Dstem (g C m −2 yr −1 ) and previous July-August, growing-year July and growing-year May-June temperature values, respectively (Fig. 4b; these coefficients are for comparison with those in Table 3).

Discussion
In this study, the MAIDEN model was successfully modified to consider important processes for boreal tree species and to improve the simulation of the coupling between photosynthesis and carbon allocation to the stem in boreal forests.Because we used a Bayesian optimization procedure, we start the following discussion with the interpretation of the parameters' posterior distributions (Sect.4.1) and of the simulation uncertainties (Sect.4.2).Subsequently, some model predictions at the 2050 horizon are presented to identify the likely response of the studied boreal forests under future environmental change (Sect. 4.3).Finally, we conclude by illustrating factors that may potentially influence the obtained results (Sect.4.4).

Parameter interpretation
The posterior distributions of the parameters were relatively sharp (Figs.S4 and S5; Table 1; by sharpness, we mean the shrinking of the distribution relative to the prior acceptable range toward a posterior distribution with a well-defined, narrow peak).Sharp distributions with small posterior ranges relative to the prior ones indicate sensitive parameters.This result means that the model posterior probability (i.e., model plausibility) increased significantly with the specific values of the selected parameters retained by the MCMC sampling.The slightly bimodal structures of the posterior distributions of Vmax, Vb and Vip were likely a consequence of their significant cross-correlations (Table S1).However, the posterior distributions of these three parameters were robust and consistent even when the Bayesian optimization was executed on independent periods (Fig. S19).The optimization of some parameters controlling Dstem (the three related to the start of the growing season and Cbud) was sensitive to the choice of the period and the site in the cross-validation exercise (Figs.S20 and S21), likely as a result of the short length of the available observed data (61 yearly RWhighF values) and of significant cross-correlation coefficients (Table S2).However, in all cases, the uncertainties in the parameters' posterior distributions (Figs.S4 and S5) did not affect our interpretations because the MAIDEN simulations were extremely consistent irrespective of the selected block of parameters (see Figs. 3 and 5).
The interpretation of some parameters needs specific attention, such as the parameters controlling the negative impact of both previous-year April precipitation and July-August temperature values on canopy development.Warm previous Aprils with infrequent late snowfalls may accelerate snowmelt and the start of the previous growing season, allowing optimal reserve accumulation during the previous year, which then would influence tree performance the following growing year.This mechanism may be significant, especially if we do not observe high temperatures limiting soil water availability and reserve accumulation during the previous summer (Girardin et al., 2016).It has already been shown that shoot elongation of boreal conifers is determined by climate conditions during bud formation (Salminen and Jalkanen, 2005).However, for Scots pine, previous summer www.biogeosciences.net/14/4851/2017/Biogeosciences, 14, 4851-4866, 2017 temperatures are positively correlated with shoot elongation, while in our case, the opposite process was simulated, and the simulations were even more sensitive to the values of the CanopyT temperature-dependent parameter than to those of the CanopyP precipitation-dependent parameter (Fig. 7b-e).
We need more data on canopy development and shoot elongation to verify the model results.

Interpretation of model performance
The comparisons with the observed data suggest that MAIDEN as revised produces accurate simulations of GPP, of ring width and stem biomass variability and of intragrowing season dynamics.The explained variance (R 2 = 0.73 and r = 0.86, df = 59, p < 0.001; Fig. 5b) is higher than that explained by MAIDEN for Mediterranean sites (R 2 slightly above 0.5; Gea-Izquierdo et al., 2015).However, the ensembles of daily and annual time series retained by the MCMC sampling were not always centered on the observed time series (Figs. 3 and 5).Specifically, the simulated annual GPP values often underestimated the actual GPP, especially at low observed GPP.This result reflects the fact that the MCMC sampling maximized the model plausibility according to the model structure and, by doing so, retained similar blocks of parameters.Thus, the range of simulated values in Figs. 3 and 5, obtained with all retained iterations, should be interpreted as the uncertainty due to only parameter selection, while the uncertainty due to the non-perfect fit between observations and simulations was not considered.
We drew the following conclusions based on the response function analysis.First, in the case of daily GPP (Table 2), MAIDEN performed better than response functions, suggesting that it properly simulates climate-driven processes governing photosynthetic assimilation, which are well known to be a result of several nonlinear processes.Second, most of the variance explained by the response functions was due to temperature variables, reflecting the greater sensitivity of northern black spruce forests to temperature compared to drought stress (Gennaretti et al., 2014) and justifying the modeling in MAIDEN of the maximum carboxylation rate as a function of temperature.Third, only temperature vari- ables of preceding days were retained, justifying the inclusion of our acclimation function of photosynthesis into temperature to increase the influence of previous days.Fourth, the coefficient estimate for precipitation of lag 0 (i.e., week ending in day 0) was negative, while that of lag −2 was positive, even though these variables share 5 of 7 days of data.The reduction of absorbed photosynthetically active radiation associated with cloudiness during raining days could explain this result.In the case of annual stem growth (Table 3), the explained variability with the best response function was greater than with MAIDEN, suggesting that the process-based modeling can potentially be improved with additional data and by including stronger legacy effects of the year preceding ring formation (Girardin et al., 2016).Indeed, most of the variance explained by the response function was due to a negative correlation with the temperature of the previous summer.Contrasting correlations with summer temperature values of the previous and the current growing year are also visible in Fig. 4a and have already been observed for black spruce (Mamet and Kershaw, 2011;Ols et al., 2016).

Model predictions
It is possible to use the new optimized MAIDEN to predict forest growth and allocation dynamics of the studied boreal forests under future environmental change.At the 2050 horizon, daily maximum temperature, daily minimum temperature and precipitation within the study area should increase by approximately 2.3 • C, 4.3 • C and 12 %, respectively (Guay et al., 2015).If we modify the climate data used (Sect.2.2.3) with these median changes and we fix the CO 2 concentration at the 2010 level, then the median increase in the annual GPP and Dstem values simulated by MAIDEN for the studied forests is 43 and 68 %, respectively (Fig. S22).It is important to note that the ring width data used for the optimization of MAIDEN come from lake riparian trees and that these results are too optimistic for more water-limited boreal sites.

Limits and error sources of the study
Although the simulated results with MAIDEN were satisfactory, we have to consider two important limits and error sources of the study.First, for the optimization of carbon allocation, we assumed that stem biomass (or carbon) increments were proportional to ring growth.This approach was necessary because data from field plots were not available from all study sites.A recent study showed that the maximum rate of ring width increase during the growing season precedes the maximum rate of increase in wood biomass and that these processes could exhibit differential sensitivities to local environmental conditions (Cuny et al., 2015).However, Cuny et al. (2015) also highlighted that wood biomass production follows the seasonal course of temperature in coniferous forests, and this is what we observed once MAIDEN was optimized.Indeed, almost all available carbon in spring was allocated to the canopy and roots (Fig. S5g; Table 1), whereas C allocation to the stem (Dstem) in summer increased with temperature (Fig. 8).Furthermore, the ring width series used were highly correlated with July-August temperature, as expected for wood biomass production and for climate-growth analysis for the studied species.Second, some fixed parameters are present in the MAIDEN code (see Eqs. 4 and 6).These parameters may be potentially modified, but their specification is justified in Sect.2.1.2and limits additional parameter tuning.Third, we modeled GPP and carbon stem increments of a boreal tree species using mechanistic rules, which increased the capability of MAIDEN to reproduce observed variations.However, our choice of mechanistic rules was subjective to some extent and depended on previous physiological knowledge and on model-data comparisons.Such model refining is an important step of all model-data fusion approaches (Guiot et al., 2014) and increases our understanding of ecosystem functioning and responses.Nevertheless, the proposed mechanistic rules should F. Gennaretti et al.: Ecophysiological modeling of photosynthesis and carbon allocation be verified in the future with additional data from a wider boreal area.

Conclusion
In this study, we adapted a process-based forest ecophysiological model developed for temperate and Mediterranean forests to simulate gross primary production and stem biomass increment for black spruce, the dominant species across the North American boreal biome.The model used, MAIDEN (Misson, 2004;Gea-Izquierdo et al., 2015), has the specificity to simultaneously simulate the course of photosynthesis and phenological phases characterized by specific allocation rules dependent on climatic conditions.The model represented the tree-ring interannual variability even though detrended radial growth was poorly explained by the simulated annual GPP (Fig. 2b-c), which suggests that the relationship between GPP and wood production is complex and nonlinear (Rocha et al., 2006).Significant simulation improvements were obtained, introducing important processes for temperature-sensitive boreal forests into the model, such as (i) the acclimation of photosynthesis to temperature over several days (see Gea-Izquierdo et al., 2010;Mäkelä et al., 2004), (ii) the influence of previous-year climatic conditions affecting bud formation on the potential amount of carbon allocated to the canopy each year (see Salminen and Jalkanen, 2005) and (iii) the positive relationship between temperature and the carbon allocated to the stem in summer (see Cuny et al., 2015).Although we used black spruce data from the northern Quebec taiga to test and optimize the model, the new model modifications have the potential to work within other boreal regions and tree species.The effects of the introduced functions can be amplified, reduced or canceled in the Bayesian optimization procedure according to the relevance of specific processes in the studied forest.Boreal ecosystems are crucial carbon stores that must be quantified and preserved (Bradshaw et al., 2009).Their future evolution is extremely important for the global carbon budget.Development of process-based models, such as the one used and improved here, combined with continuous field data acquisition will help determine the role of the different environmental factors and underlying mechanisms in present and future boreal forest carbon fluxes.In this context, we believe that our study helps to explain how boreal forests assimilate and allocate carbon depending on weather/climate conditions.
Figure1.MAIDEN-simulated phenology (blue), water (black) and carbon (red) fluxes.AN: net photosynthesis corresponding to net primary production.Cstored, Cstem, Ccanopy and Croots: carbon allocated daily to stored non-structural carbohydrates, stem, canopy and roots, respectively.DOY: day of the year (1-365).GDD: growing degree days.f 3 and f 4 : functions determining carbon allocation in phases 3 and 4. Cbud: amount of storage carbon that is used each day by the plant in phase 3.

Figure 2 .
Figure 2. Variance explained by the model.(a) R 2 between observed and simulated GPP daily values.R 2 (computed on data transformed to z scores) between the mean of the detrended series of black spruce ring growth (RWhighF) and simulated yearly detrended C allocation to the stem (b) or GPP (c).Vertical dashed line is the mode, and the blue line is the value with plausible block GPP (in a) or with plausible block stem (in b, c).All probability density functions are based on 50 simulations.

Figure 3 .
Figure 3.Comparison between observed GPP values and MAIDEN-simulated values at the Quebec Eastern Old Black Spruce site.(a) Daily values (units are µmol C m −2 day −1 ).In the scatterplot (R 2 = 0.90; r = 0.95, df = 2918, p < 0.001), observations are compared with the values obtained with plausible block GPP.In the annual cycle plot, black lines are the medians and the 5th and the 95th percentiles of the simulated values from all iterations retained by the MCMC sampling.(b) Annual values (units are µmol C m −2 yr −1).In both plots, observations are compared with the values from all iterations retained by the MCMC sampling.In the scatterplot, the R 2 of the data is 0.31 (r = 0.76, df = 6, p < 0.05).

Figure 4 .
Figure 4. Correlation between monthly climate variables of the study area (precipitation in blue and mean, maximum and minimum temperature in red, violet and orange, respectively) and the mean of the detrended series of black spruce ring growth (RWhighF; a) and the simulated detrended annual carbon allocation to the stem (Dstem; b).For the climate variables, time windows of 31 days are used to obtain time series of monthly data (over the 1950-2010 period) for each day (central day), averaging the values of each window and each year.These climate time series are then detrended.Thresholds of significance (p < 0.05) are shown by horizontal dashed lines.

Figure 5 .
Figure5.Comparison between the observed mean series of black spruce ring growth (unitless growth indexes) and MAIDEN-simulated carbon allocation to the stem (Dstem; g C m −2 yr −1 ).(a) Detrended series (in the scatterplot the R 2 computed on data transformed to z scores is 0.24; r = 0.62, df = 59, p < 0.001).(b) Series with their trends (in the scatterplot the R 2 is 0.73; r = 0.86, df = 59, p < 0.001).In all plots, observations are compared with the values from all iterations retained by the MCMC sampling.

Figure 6 .
Figure6.Influence of the temperature transformation (S) on the modeled annual cycle of GPP daily values (µmol C m −2 day −1 ) at the Quebec Eastern Old Black Spruce site.Only the τ parameter determining the S values was allowed to vary, while the other parameters were fixed to the values of plausible block GPP.(a) τ is 1 day (S same as Tday).(b) τ is 6.71 days (a middle value between 1 and 12.43).(c) τ is 12.43 days (same τ than in plausible block GPP).(d) τ is 18.14 days (a higher value than in plausible block GPP).The R 2 between observations and simulations is reported in each plot.Panels (e-h) show the impact of the respective τ values on S if the daily Tday time series corresponds to a single step of 10 • C lasting 30 days.

Figure 7 .
Figure 7. Temperature and precipitation dependence of CanopyMult (a; unitless multiplier; Eq. 4; CanopyT and CanopyP are those of plausible block stem), which determines the yearly canopy potential amount of carbon.Previous-year mean July-August temperature indexes are on the x axis, and previous-year April precipitation indexes are on the y axis.Black dots are observed values in the central point of the region with ring width data.Panels (b-e) show how CanopyT and CanopyP (varying over their prior acceptable ranges) impact the correlations between simulated Dstem (g C m −2 yr −1 ) and observed ring width data (RW or RWhighF; unitless indexes), when all other parameters are fixed to the values of plausible block stem.The vertical lines are the selected values for plausible block stem, and the horizontal lines are the 90 % confidence intervals based on the parameters' posterior densities (Fig. S5).

Figure 8 .
Figure 8. Temperature dependence of the daily partition of carbon in phase 4 (growth and accumulation phase in summer) when MAIDEN is run with the parameters of plausible block stem at the center of the region with ring width data in the northern Quebec taiga.(a) Probability density of daily maximum temperature values in summer.(b) Relationship between maximum temperature values and portion of carbon allocated to the stem (Eq.6).The vertical dashed lines show the range of maximum temperature values.Panels (c-d) show how the parameter st4temp influencing this process impacts the correlations between simulated Dstem and observed ring width data (RW or RWhighF) when all other parameters are fixed to the values of plausible block stem and st4temp varies over its prior acceptable range.The vertical line is the selected value for plausible block stem, and the horizontal line is the 90 % confidence interval based on the parameter's posterior density (Fig.S5).

Table 1 .
Definitions, symbols and prior and posterior ranges of calibrated parameters.Small posterior ranges relative to the prior ones indicate sensitive parameters.

Table 2 .
ANOVA table for the best response function (here, a combination of 4 of the 10 tested predictors minimized the BIC) with daily GPP at EOBS as dependent variable (excluding days between 15 November and 1 April).All F values are highly significant (p < 0.001).For precipitation, lag n indicates the sum of the daily precipitation of the week ending in day n.

Table 3 .
ANOVA table for the best response function (here, a combination of 3 of the 10 tested predictors minimized the BIC) with the observed mean detrended ring width series (RWhighF) as the dependent variable.