The Role of Respiration in Estimation of Net Carbon Cycle : Coupling Soil Carbon Dynamics and Canopy Turnover in a Novel Version of 3 D-CMCC Forest Ecosystem Model

Understanding the dynamics of Organic Carbon mineralization is fundamental in forecasting biosphere to atmosphere Net Carbon Ecosystem Exchange (NEE). With this perspective, we developed 3D-CMCC-PSM, a new version of the hybrid Process Based Model 3D-CMCC FEM where also heterotrophic respiration (Rh) is explicitly simulated. The aim was to quantify NEE as a forward problem, by subtracting Ecosystem Respiration (Reco) to Gross Primary Productivity (GPP). To do so, we developed a simplification of the Soil Carbon dynamics routine proposed in DNDC [1]. The method calculates decomposition as a function of soil moisture, temperature, state of the organic compartments, and relative abundance of microbial pools. Given the pulse dynamics of soil respiration, we introduced modifications in some of the principal constitutive relations involved in phenology and littering sub-routines. We quantified the model structure related uncertainty in NEE, by running our training simulations over 1000 random parameter-sets extracted from parameters distributions expected from literature. 3D-CMCC-PSM predictability was tested on independent time series for 6 Fluxnet sites. The model resulted in daily and monthly estimations highly consistent with the observed time series. It showed lower predictability in Mediterranean ecosystems, suggesting that it may need further improvements in addressing evapotranspiration and water dynamics.


Introduction
Global concerns over increasing levels of greenhouse gas concentrations, particularly carbon dioxide (CO 2 ), have pushed research efforts to better investigate biogeochemical carbon (C) flux dynamics and patterns between atmosphere and biosphere, and to upscale C flux estimates from site-specific to regional, continental and global scales.Increased atmospheric concentration of CO 2 , combined with increasing temperatures and size variations of ecosystem C pools, are responsible for year-to-year terrestrial ecosystem carbon flux perturbations, through the variation of both photosynthetic and respiration rates [1].
In the last decades, the eddy covariance (EC) technique has provided long-term continuous measurements of net carbon ecosystem exchange (NEE), water vapor and energy, within the global network of EC flux towers (FLUXNET) distributed over major terrestrial ecosystems.The availability of EC measure of NEE contributed to quantify and to determine seasonal and inter-annual variability of ecosystem C budgets at EC tower site-specific scale [2][3][4].Observed NEE does not directly quantify the two major components of ecosystem C flux balance represented by ecosystem respiration (R eco ) and gross primary productivity (GPP).Thus, flux partitioning algorithms have been developed to partition eddy covariance NEE into photosynthetic uptake and respiratory release [5,6].At the same time, EC flux measurements provide key information for the parameterization, calibration and validation of process-based forest ecosystem models (FEMs) contributing to large-scale estimates of main ecosystem C pools.
The implementation of both forest process-based models (PBMs) [7][8][9][10][11][12] and functional-structural tree models (FSTMs) [13][14][15][16], based on the widely used light use efficiency (LUE) approach [17], has contributed to understanding and upscaling the main physiological processes supporting ecosystem C uptake.Although most of forest ecosystem models provide reliable estimates of forest growth, limitations for NEE estimates are related to the uncertainty in R eco estimation [18,19].Hence, the implementation of biogeochemical models integrating soil respiration models and FSTMs is a great opportunity to reliably estimate NEE [20][21][22].
Although soil respiration and soil organic matter (SOM) decomposition depend mainly on abiotic factors, such as temperature and soil moisture [23][24][25], a key role is played by soil organic carbon (SOC) stock size, microbial pools [26,27], and dynamics of SOC supply to soil with littering [28].Leaf and fine root senescence is of primary importance in determining dynamics in heterotrophic respiration [29,30], and exhibits seasonal patterns and intra-seasonal pulses [31,32].These pulses are mostly driven by phenological transitions through stages of dormancy, active growth, and senescence.For example, supply of dead leaves in soil is strongly dependent on a tree's leaffal strategy.This is also true when predicting tree respiration, which is directly related to their growth and nitrogen content, which depends on spring phenology.Unfortunately, processes involving budburst and senescence are still partly obscure.PBMs usually represent these processes simplistically; these simplifications may lead some terrestrial ecosystem models to result in biased predictions [33][34][35].
For these reasons, we developed the 3D-CMCC-PSM (3D-CMCC-Phenology and soil Model), a new version of the hybrid process-based model 3D-CMCC FEM proposed by Collalti et al. [36,37], where (1) spring phenology was directly taking into account tradeoffs between growth and Non Structural carbon (NSC) demand; (2) phenological transitions and supply of Fresh organic matter (FOM) to soil were more explicitly represented; and (3) heterotrophic respiration (R h ) was explicitly simulated to dynamically quantify stock changes of 7 different SOC pools mediated by the amount of active microbial C pools, following the rationale proposed in the DNDC [38].The aim of this study was to: (1) test the performance of the modified model version, comparing model NEE estimates against independent time series for 6 Fluxnet sites, representing different forests in different climatic areas, distributed over a wide latitudinal gradient amongst European EC sites; and (2) quantify uncertainty associated to 3D-CMCC-PSM constitutive relations structure and parameterization.

Study Area and Data
Eddy covariance data were collected from FLUXNET [39].We chose the 6 sites to represent a climatic and longitudinal transect through Europe (Figure 1), so that the model could be tested on Forests 2017, 8, 220 3 of 24 different critical boundary conditions (Table 1).We used EC data from different time series from daily to annual, processed using the method described in [40].Gross primary production (GPP) and ecosystem respiration (R eco ) were partitioned using [6].Information about forest structure and total SOC at the beginning of the simulation was collected from literature (e.g., [41,42]) and PIs information.Sites were chosen to represent 3 diverse forest ecosystems, dominated by different species composition from deciduous broadleaved, DBF (i.e., Fagus sylvatica L.), evergreen broadleaved, EBF (i.e., Quercus ilex L.), and evergreen needle leaved, ENF (i.e., Pinus sylvestris L. and Picea abies L.), representing the most common European forest species from boreal to mediterranean ecoregions across Europe.

Model Description
The 3D-CMCC FEM is a stand-scale process-based model (PBM) designed to simulate C and water cycle in natural and managed forest ecosystems (for a full description see [36,37]).Several eco-physiological processes were modeled at species-specific level, and at a variable temporal scale (from daily to annual) depending on the process to simulate (Figure 2).Model outputs were generally represented at hectare scale, while processes were simulated at different spatial scales from cellular (e.g., stomatal conductance), to canopy (e.g., transpiration), to individual tree, up to stand level.
Carbon assimilation was modeled for sun and shaded leaves using the LUE approach [43]: potential C assimilation was constrained by environmental and stand structural (e.g., tree age) scalars [44].Autotrophic respiration (R A ) was explicitly modeled as the sum of growth and maintenance respiration (R G and R M , respectively).The first was computed as a fixed ratio of new growth tissues (30%) and the latter was based on nitrogen content in stems, branches, leaves, fine and coarse roots, non-structural carbon (NSC), and fruit tree pools.carbon allocation among these pools was controlled by species-specific parameters, phenology, light and water availability.Water cycle was modeled calculating the daily balance between precipitation, canopy transpiration, evaporation, soil evaporation, and runoff.Meteorological variables used to force the model were: global solar radiation (MJ m −2 day −1 ), maximum and minimum air temperature ( • C), relative humidity (%), and daily cumulated precipitation (mm day −1 ).To be initialized, the model required knowing stand structural characteristics such as: species composition, stand density, diameter at breast height (DBH), tree height and age.Soil initialization required the estimation of total organic carbon (TOC) in the different SOC pools, as described in the following section.calculating the daily balance between precipitation, canopy transpiration, evaporation, soil evaporation, and runoff.Meteorological variables used to force the model were: global solar radiation (MJ m −2 day −1 ), maximum and minimum air temperature (°C), relative humidity (%), and daily cumulated precipitation (mm day −1 ).To be initialized, the model required knowing stand structural characteristics such as: species composition, stand density, diameter at breast height (DBH), tree height and age.Soil initialization required the estimation of total organic carbon (TOC) in the different SOC pools, as described in the following section.

Soil Carbon Dynamics
The most recent 3D-CMC FEM model version (v.5.1, [37]) lacks in representing SOC dynamics, preventing any estimation of NEE.With that perspective, we developed a simplified version of the method described in [38] to quantify, dynamically, changes of 7 different SOC pools mediated by the amount of active microbial C pool (Figure 3).

Soil Carbon Dynamics
The most recent 3D-CMC FEM model version (v.5.1, [37]) lacks in representing SOC dynamics, preventing any estimation of NEE.With that perspective, we developed a simplified version of the method described in [38] to quantify, dynamically, changes of 7 different SOC pools mediated by the amount of active microbial C pool (Figure 3).The litter C decomposed by microbial activity is partly mineralized as CO2, partly stored into microbial metabolic biomass (labile), partly in structural microbial biomass (resistant), and partly transformed in other organic compounds [38].SOC stability is also related to its chemical recalcitrance, its accessibility, and interaction with clays [45].We divided SOC in 7 pools to take these differences into account.Humic pool (we use the term Humads, to be consistent with [38]) was divided into a more labile (labile Humads, which stands for Humic acid), an intermediate (resistant Humads, which stands for Fulvic acid) and a more resistant sub-pool (Humus, which stands for Humine).
We assumed that microbial C use efficiency was different for different pools, but constant within each pool [44].Microbial activity is strictly related to micro-environmental conditions too.Given these premises we modeled C dynamics among soil pools using the following two partial differential The litter C decomposed by microbial activity is partly mineralized as CO 2 , partly stored into microbial metabolic biomass (labile), partly in structural microbial biomass (resistant), and partly transformed in other organic compounds [38].SOC stability is also related to its chemical recalcitrance, its accessibility, and interaction with clays [45].We divided SOC in 7 pools to take these differences Forests 2017, 8, 220 6 of 24 into account.Humic pool (we use the term Humads, to be consistent with [38]) was divided into a more labile (labile Humads, which stands for Humic acid), an intermediate (resistant Humads, which stands for Fulvic acid) and a more resistant sub-pool (Humus, which stands for Humine).
We assumed that microbial C use efficiency was different for different pools, but constant within each pool [44].Microbial activity is strictly related to micro-environmental conditions too.Given these premises we modeled C dynamics among soil pools using the following two partial differential equations: with ∂t being the CO 2 efflux produced from a specific C pool decomposition, C L (t) is the amount of new carbon entering the specific soil pool, C p (t) the amount of carbon in that C pool, B(t) the microbial biomass competing for C(t).α and β respectively represented the microbial turnover and SOM consumption factors.α was treated as a constant value, as in [38].
The consumption factor β was estimated as a function of both SOM stability and micro-environmental conditions: where µ T represented the temperature factor, µ M the moisture factor, µ A the accessibility, k pool the recalcitrance, and µ D a clay dependent depth factor [46].The k pool factor was treated as a specific parameter depending on each compartment's biomass.T represented temperature in Celsius degrees, θ stands for water moisture, z the mean depth of the soil, clay % is the percentage of clays in soil.Humads were decomposed using the same rationale, but had slower rates.Inert organic matter (IOM) was calculated following [47].

Deciduous Phenology
Similarly to [37] the phenology scheme was constrained in 3 and 5 sub-phases respectively for evergreen and deciduous species (Table 2).These phases were driven by photoperiod, thermal sum, and maximum leaf biomass (resulting from maximum attainable Leaf Area Index, LAI, m 2 m −2 ).
3D-CMCC FEM represents leaves development as a by-product of leaf biomass pool dynamics [48].Despite being a reasonable simplification, such method leads to non-negligible excesses in growth respiration and NSC consumption during budburst.Such unrealistic demand could eventually consume all the available carbon, causing the death of the tree.For this reason, we proposed a modification of budburst phenology, to explicitly simulate the dynamic tradeoffs between demand in NSC, increase in foliar biomass, and maturation of progressively self-sufficient leaves.We based the tradeoff function on the hypothesis that new leaf demand of NSC is higher during budburst, and gets progressively lower when maturing [50].Idealistically, the total NSC mass (B R ) demanded by maturing shoots and leaves could be represented by the linear differential equation (ODE): where R NSC (t) is the instantaneous proportion of NSC demand.We assumed that C request per leaf exponentially decreases with maturity, while the total demand increases by the increasing number of leaf primordia.For simplicity, we assumed that the two components resulted in the linear reduction of R NSC (t).Being a fraction, the maximum value of the integral of R S&F&R (t) is equal to 1. Being the first vegetative day t 0 , and the last day of budburst (BB T ) the last possible one to reach complete leaf development, the domain of the function was [t 0 , BB T ]: resolving the ODE, and substituting R S&F&R (t) it gives: where BB T is the parameter used in [37], e  it gives: where is the parameter used in [37], is the biomass dependent and the maturity dependent factors (expressed in days).Graphically, the equation represents a skewed function of the amount of NSC allowed to be used by the trees of the specific class; the faster the leaves reach maturity, the more daily specific allocation is allowed (Figure 4). .The shorter BBT, the higher the maximum NSC fraction.
Another major modification we introduced involves fall leaf yellowing and senescence.Falling leaves contribute to half of annual litter production.Correctly estimating their timing has a strong impact on both GPP and Rh dynamics [32].The 5.1 version of 3D-CMCC-FEM represents senescence by linearly decreasing leaves biomass of a predefined fraction until their pool is emptied.However, such method systematically overestimated leaf turnover at the beginning of the leaffall phase, and is poorly flexible in calculating its duration.For these reasons, linear loss of leaf biomass was replaced with a sigmoid function.Assuming for hypothesis that all leaves fall by the end of senescence season, Another major modification we introduced involves fall leaf yellowing and senescence.Falling leaves contribute to half of annual litter production.Correctly estimating their timing has a strong impact on both GPP and Rh dynamics [32].The 5.1 version of 3D-CMCC-FEM represents senescence by linearly decreasing leaves biomass of a predefined fraction until their pool is emptied.However, such method systematically overestimated leaf turnover at the beginning of the leaffall phase, and is poorly flexible in calculating its duration.For these reasons, linear loss of leaf biomass was replaced with a sigmoid function.Assuming for hypothesis that all leaves fall by the end of senescence season, the sigmoid function was: where α(t), β(t), γ(t) are three parameters with biological or physical meaning.In fact: where L 0+δt stands for LAI value at peak of green, t 0 (s) is the first day of senescence (triggered by a species-specific day-length parameter), ∆t(s) is the length of the senescence period (days).∆t(s) was calculated using [49], as a function of temperature and photoperiod: where T Max is the maximum temperature at which senescence is effective, T(t) is daily average Temperature, D L (0) photoperiod at the first day of senescence, D L (t) photoperiod at the ith day of the year.

Evergreen Phenology
Evergreen canopy turnover was modified from [37].3D-CMCCFEM v.5.1 assumes that evergreen leaf turnover is constant throughout the year, and that annual leaf turnover is equal to leaf biomass produced the year before [37,48].However, leaf turnover seems to be concentrated in specific seasonal windows: (1) consistently to leaf emergence (spring); and (2) approaching of photosynthetic inefficient season (fall) [51].For this reason, this approach may affect the ability in estimating leaf biomass and GPP intra-annual variability.To better represent leaf turnover dynamics, we developed a new framework where competition for light dynamically affected leaf turnover.We assumed the canopy to be a population of leaves optimized to intercept the highest amount of light.Since leaves cannot move from their position in the canopy, they get partially shaded by new emergent leaves when aging.Conceptually, leaves can be assumed to be in competition for one resource, light, and their turnover can be predicted by using a competition for one resource scheme as in [52]: In this model, R and R i * are respectively the concentration of available light in J m −2 day −1 (i.e., resource that leaves are competing for), and the light needed to survive in a progressively more shaded canopy.r i is the max photosynthetic rate (gC m −2 day −1 ), m i is the maintenance respiration in gC m −2 day −1 , Y i carbon yield.We assumed for hypothesis that: Forests 2017, 8, 220 9 of 24 (1) Older leaves live in the shaded portions of the canopy, where light transmitted is reduced following Lambert Beer's exponential decay equation.For this reason, we expect an exponential reduction in absorbed photosynthetically active radiation (APAR); (2) An age dependent quasi-exponential decay in leaf quantum yield efficiency [53].These decays impact on the reduction of r i ; (3) Nitrogen content in older leaves is often lower than in young ones, because of its transfer from portions of the crown with low productivity to portions more exposed to light [53].Since maintenance respiration is proportional to nitrogen content, we expect an exponential reduction in m i ; (4) Y i was assumed to be constant as in [17] because of the joint effect of reduction in respiration rate and quantum yield efficiency.
We assumed that the three components of the Equation ( 9) may have the following shapes: where t is time, the independent variable.m 0 , α 0 , λ 0 are MR, quantum yield and APAR at first day of budburst, the second year of life.m m , α m , λ m are the theoretical minimum values of MR, quantum yield and APAR to let a proportion of leaves survive in a non-shading context.k η , k α , k λ are the exponential parameters for the three functions.
Based on these hypotheses Equation ( 10) becomes: Assuming that k η = k α = k λ , and that the denominator of Equation ( 9) has to be greater than zero for the leaves in the ith generation to survive, R* is a sigmoid positive function.Knowing R*, we can calculate the amount of live leaf biomass of the ith generation ( (i) ) as the inversion of R* (i.e., S * i (t)).We simplified the theoretical model using the following function: where t is the number of days since leaves of the ith generation have emerged.According to this model the theoretical maximum age of each generation (BF LS(i) ) should correspond to the year in which the R* almost reaches its asymptote.We used Equation ( 12) to quantify leaf turnover for each generation.About 60% of leaf turnover happened in early spring, when new leaves emerged: the amount of biomass lost every day was proportional to new leaf biomass production the same day.The rest of annual turnover happened in fall, when each leaf biomass generation was reduced of a constant value calculated from Equation (12).No foliar re-sprouting was simulated in fall, even though there are evidences of it for Quercus ilex [54].Leaf biomass reduction was determined by linearly decreasing each Bi to the quantity predicted by the specific parabolic decay for the end of the year.

Production of Fresh Organic Matter
At a relatively fine temporal scale (i.e., daily time step) timing of litter formation and FOM production may be fundamental in correctly estimating Rh [55].Littering for woody tissues followed the rationale of BIOME-BGM family [48].Partitioning of leaves and fine root turnover followed [56].FOM coming from any plant biomass pool to the soil was added to the litter pool.Litter pool was consisting of three sub-pools: metabolic very labile C, structural labile C, and structural resistant C. FOM was divided into the three litter sub pools according to its original C:N ratio, as in [38]: where C lt is the FOM entering litter from each structural C compartments of the plant.FOM C and N were distributed to the litter sub pools with closest C:N.C i+1 , represents the pool with higher recalcitrance.When CN lt was higher than any litter sub pool, all the new C was added to the structural resistant pool; otherwise, if CN lt was lower than the CN of the metabolic pool, all its C and N were added to the very labile sub pool.Litter C dynamically moved from a pool to another.Microbes absorbed and partially immobilized litter C in their biomass, and released it again in the soil, during the humification processes [57].

Optimization
We introduced a new calibration scheme to provide an optimized parameterization, and quantify uncertainty related to the equations used in 3D-CMCC-PSM to simulate NEE, and choice of a subset of 45 species-specific physiological parameters (Table S1).Calibration was performed on each site independently, using a training dataset composed by 3 years of EC daily NEE time-series.We chose 2000-2002 for DEHai, 2001-2003 for FIHyy, ITCpz, and FRPue; 2001-2004 for ITCol; 2006-2008 for ITRen.Because of the strong autocorrelation characterizing time-series, we didn't randomly split our dataset, and decided to use sub-time series of contiguous years [58].We preferred to use the beginning of each simulation for training, before the global recession following the Financial Crisis of 2007-2008.We decided to do so to (1) facilitate comparison among different sites' simulations; (2) have the calibration start when vegetation structure was known from literature; and (3) reduce the risk for parameters to be fitted over a changing environment instead of eco-physiological properties.
To sample the parameterization space (the realistic values of each physiological parameter of the species simulated in each site), we randomly extracted 1000 parameter-set combinations from prior distributions.Prior distributions were assumed to be the same among individuals of the same species, different across species.We assumed each parameter to follow a truncated normal distribution, to avoid any possibility to have non-realistic negative values.Average and variance were estimated by using values found in literature, as in [37].We used the same averaged value as in [37] for those parameters whose observations in literature were less than 3, because we didn't have enough information to calculate sample standard deviation.The optimization was performed by choosing the parameterization set maximizing the objective function QF through: where  To evaluate the model efficiency, we calculated daily, monthly, and seasonal: (1) R; (2) NSE; (3) root mean square error (RMSE); and (4) mean absolute bias (MAB).Each statistic was considered differently informative [59] as summarized in Table 3.The model's ability in representing observed anomalies was determined by analyzing inter annual variability (IAVs) following [60] and [37].

Statistics
Formulation Use and Ranges Estimation of model's measure of correlation with EC data [0;1] Nash Sutcliffe efficiency

Evaluation of Daily, Seasonal, and Annual NEE Estimations
To evaluate 3D-CMCC-PSM NEE predictions, we compared predicted (MD) daily and monthly NEE time series to EC daily data.The analyses were performed only on the test data (i.e., portions of the series which have not been used for calibration) to avoid any effect of overfitting.The model showed high correlations with observed EC data at all sites for both daily and monthly fluxes, apart from ITCpz site (Table 4).Excluding ITCpz, R ranged at all sites from 0.65 to 0.84 for daily, and 0.59 to 0.97 for monthly scale.Beech-dominated Deciduous Forests (DBF) performed better than Conifer species (ENF) and evergreen Mediterranean broadleaved forests (EBF).ENF and EBF in FRPue performed similarly on daily scale, for all the statistics used.However, ENF predictability significantly increased on monthly scale (R ranging between 0.92 and 0.97), while EBF performed worse (R 0.42 in ITCpz, and 0.59 in FRPue).RMSE on average was 1.92 gC m −2 day −1 .MAE ranged between 0.96 and 1.78 gC m −2 day −1 , and on average it decreased almost twice on monthly timescale.MAB showed similar behavior for DBF and ENF.It ranged between 0.39 and 0.56 gC m −2 day −1 (0.50 on average) for daily time series.Mediterranean forests resulted as the ones with highest MAB, and showed no significant reduction when predictions were aggregated on monthly scale.Differently from the other simulations, even NSE just improved slightly for ITCpz, and even reduced for FRPue simulation.Daily results were aggregated in seasonal series to evaluate seasonal predictability.Daily NEE were averaged to give a time series of mean seasonal NEE (gC m −2 day −1 ), one value per season, for the duration of the test dataset.Seasonal aggregations showed that 3D-CMCC-PSM poorly performed in predicting seasonal fluxes.NSE was generally negative in summer, with the exclusion of DEHai and FRPue.3D-CMCC-PSM generally best reproduced NEE dynamics in fall (R ranging between 0.22 and 0.89).ENF ecosystems showed consistently higher correlation in spring predictions, with R of 0.65 and MAB of 0.62 gC m −2 day −1 on average.In the case of evergreen stands, 3D-CMCC-PSM consistently showed poor performance in summer.Expectedly, DBF performed the worst in winter (Table 5).NSE on average resulted positive only in fall for both DBF and EBF, and spring, for ENF stands.We used Taylor diagrams [61] to graphically summarize how closely Daily, Monthly and Annual NEE patterns matched EC observations (Figure 5).3D-CMCC-PSM performance was generally satisfactory.Daily simulations resulted in all sites but ITCpz being within the ±1 normalized standard deviation region.Monthly scale predictions were more consistent with EC data, especially for BDF and ENF sites.It resulted in all 4 simulations falling within ±0.5 normalized standard deviation from the reference point, and R > 0.9.Again, 3D-CMCC-PSM performed worst in EBF, with FRPue still Forests 2017, 8, 220 13 of 24 inside ±1 normalized standard deviation region, and ITCpz falling outside the ±1.5 normalized SD region.The consistently worse predictability in ITCpz and FRPue confirm a systematic weakness in 3D-CMCC to represent fluxes for these sites as already described in [37].Model performance on annual scale showed a different pattern, mostly because of some sites' consistent biases in seasonal NEE, and the difference in NEE magnitude.Delay in spring phenology, and the consistent underestimation of summer NEE, resulted in significant underestimation and scarce predictability of ITCol annual NEE (R < 0.2).ITCPz and FRPue resulted among the sites with higher annual predictability, partly because of the low seasonal variance in NEE, partly because winter and spring bias tends to compensate each other.
We used Taylor diagrams [61] to graphically summarize how closely Daily, Monthly and Annual NEE patterns matched EC observations (Figure 5).3D-CMCC-PSM performance was generally satisfactory.Daily simulations resulted in all sites but ITCpz being within the ±1 normalized standard deviation region.Monthly scale predictions were more consistent with EC data, especially for BDF and ENF sites.It resulted in all 4 simulations falling within ±0.5 normalized standard deviation from the reference point, and R > 0.9.Again, 3D-CMCC-PSM performed worst in EBF, with FRPue still inside ±1 normalized standard deviation region, and ITCpz falling outside the ±1.5 normalized SD region.The consistently worse predictability in ITCpz and FRPue confirm a systematic weakness in 3D-CMCC to represent fluxes for these sites as already described in [37].Model performance on annual scale showed a different pattern, mostly because of some sites' consistent biases in seasonal NEE, and the difference in NEE magnitude.Delay in spring phenology, and the consistent underestimation of summer NEE, resulted in significant underestimation and scarce predictability of ITCol annual NEE (R < 0.2).ITCPz and FRPue resulted among the sites with higher annual predictability, partly because of the low seasonal variance in NEE, partly because winter and spring bias tends to compensate each other.Simulations with R ≥ 0.9, and difference in SD with EC NEE less than 0.5 5 gC m −2 day −1 showed very good performance.Simulations with 0.75 ≤ R < 0.9, and difference in SD with EC NEE between 0.5 and 1 gC m −2 day −1 showed good performance.Simulations with 0.35 ≤ R < 0.75, and difference in SD with EC NEE between 1 and 1.5 gC m −2 day −1 showed sufficiently good performance.

Anomalies and Parameters Related Uncertainty
Figure 6 shows uncertainty associated to random choice of parameters.Overall, uncertainty was expectedly higher in summer and fall.Such increase was particularly clear for deciduous forests, which not only showed wider NEE standard deviation, but also had optimal modeled NEE falling outside standard deviation area.
DBF sites showed also high uncertainty in estimating the first vegetative day, suggesting that a better representation of winter dormancy effects on budburst dates may significantly improve model's predictability.Uncertainty was generally lower in Mediterranean sites, probably because model's performance was generally lower, and fitting parameters to data would have little effect on performance.3D-CMCC-PSM uncertainty was generally low for ENF for most of the year, but was generally high when temperatures were higher.Higher uncertainty for warmer days was generally found in DBF sites too, suggesting that 3D-CMCC-PSM was expectedly sensitive to high temperatures for both photosynthesis [62] and respiration [23]; because of such variability, calibrating parameters on data resulted in a significant boost in model's predictions.
NEE inter-annual variability was generally underestimated by the model.Nevertheless, 3D-CMCC-PSM correctly reproduced 81% of the sign of the anomalies, and residual difference in magnitude was usually less than 0.3 gC m −2 day −1 .Highest difference in magnitude occurred in ITCol (difference in residuals higher than 0.5 gC m −2 day −1 in 5 years out of 12).Highest difference was shown in ITCpz, where the sign was correctly reproduced only once out of 8 years, and having more than 1gC m −2 day −1 of residual difference (Figure 7).and 1 gC m −2 day −1 showed good performance.Simulations with 0.35 ≤ R < 0.75, and difference in SD with EC NEE between 1 and 1.5 gC m −2 day −1 showed sufficiently good performance.

Anomalies and Parameters Related Uncertainty
Figure 6 shows uncertainty associated to random choice of parameters.Overall, uncertainty was expectedly higher in summer and fall.Such increase was particularly clear for deciduous forests, which not only showed wider NEE standard deviation, but also had optimal modeled NEE falling outside standard deviation area.S2).First column represents DBF sites, the second ENF, the third EBF.
DBF sites showed also high uncertainty in estimating the first vegetative day, suggesting that a better representation of winter dormancy effects on budburst dates may significantly improve model's predictability.Uncertainty was generally lower in Mediterranean sites, probably because model's performance was generally lower, and fitting parameters to data would have little effect on performance.3D-CMCC-PSM uncertainty was generally low for ENF for most of the year, but was generally high when temperatures were higher.Higher uncertainty for warmer days was generally found in DBF sites too, suggesting that 3D-CMCC-PSM was expectedly sensitive to high temperatures for both photosynthesis [62] and respiration [23]; because of such variability, calibrating parameters on data resulted in a significant boost in model's predictions.
NEE inter-annual variability was generally underestimated by the model.Nevertheless, 3D-CMCC-PSM correctly reproduced 81% of the sign of the anomalies, and residual difference in magnitude was usually less than 0.3 gC m −2 day −1 .Highest difference in magnitude occurred in ITCol (difference in residuals higher than 0.5 gC m −2 day −1 in 5 years out of 12).Highest difference was shown in ITCpz, where the sign was correctly reproduced only once out of 8 years, and having more than 1gC m −2 day −1 of residual difference (Figure 7).S2).First column represents DBF sites, the second ENF, the third EBF.

Comparison with the 5.1 Version of 3D-CMCC-FEM
We compared the results of 3D-CMCC-FEM 5.1 [37] and 3D-CMCC-PSM for daily and monthly GPP for ITRen, FRPue, DEHai and FIHyy.We used the same non-optimized parameterization set for both versions.Except for FR-Pue, 3D-CMCC-PSM showed lower RMSE and higher R for daily GPP.

Daily and Monthly Reco
We evaluated Reco (ecosystem respiration) by comparing modeled (MD) and observed (EC) daily and monthly time series.Daily R ranged between 0.45 in ITCpz and 0.9 in FIHyy (Table 7).RMSE was of 1.28 gC m −2 day −1 on average, and MEB ranged between 0.43 and 1 gC m −2 day −1 .Most of the bias happened in summer, where Reco was generally overestimated, especially in ITCol and ITCpz.NSE was positive in any case but ITCpz.It was generally lower in DBF (0.32), and higher in FIHyy (0.60) and FRPue (0.57).Reco at monthly timescale was strongly improved in predictability, especially for ITRen, whose R increased to 0.67 and NSE to 0.84.Monthly predictions showed improvements for the other simulations too, improving R and NSE of about 0.07.Since most of the bias occurred in summer, monthly predictions showed no dramatic improvements for neither RMSE nor MEB.The only exception was ITRen, whose RMSE reduced of about 0.4 gC m −2 day −1 , suggesting that daily Reco may be noisy.

3D-CMCC-PSM Predictability in Estimating NEE
In general, the inclusion of a simplistic SOC routine resulted in a reliable estimation of daily and monthly NEE patterns.While daily and monthly patterns are consistent with EC data, seasonal patterns showed non-negligible misrepresentations, which resulted in negative NSE in most of the cases.This inconsistency may be driven by the strong seasonality in both R eco and GPP [63], which positively affects correlation between EC data and MD results.
NEE patterns during summer and fall were much more consistent with measured ones, than winter and spring patterns.During these seasons the biases appeared mostly affected by estimation of R eco .The scarcity of the model in representing EBF C fluxes was especially attributable to GPP predictions.3D-CMCC-PSM and 3D-CMCC FEM inability in predicting GPP in ITCpz and FRPue sites, denoted the necessity to better represent the relations between Mediterranean forests and environmental factors [64,65].In FRPue the model well reproduced spring, summer and fall NEE.On the other hand, it showed a bias of around 1 gC m −2 day −1 in winter, suggesting it was missing some particularly important seasonal processes.For example, evergreen phenology still didn't consider secondary or continuous growth.Thus, species like Quercus ilex, which exhibit secondary gem sprouts in fall [54], have fresh leaves and mild temperatures to guarantee photosynthetic activity in fall and winter, partly explaining 3D-CMCC FEM and 3D-CMCC-PSM systematic underestimation.
ITCpz showed the same pattern.However, differently from FRPue, it poorly performed also in early spring and summer, especially for GPP patterns.This poor performance was expected, because of the physical characteristics of the site.In fact, 3D-CMCC FEM soil water dynamics routine was still simplistic, and to date, such other similar models have not included any effect of water table dynamics.On the contrary, ITCpz is characterized with the presence of a shallow groundwater table, which seems to reduce water stress in early summer [37].Moreover, we used average daily meteo data collected from the Eddy Covariance stations, and initialized simulations using average structural information found in literature.Uncertainty about those data was potentially high, and could have dramatically affected 3D-CMCC-PSM results.
Summer NEE misrepresentation in DBF was probably affected by the assumption that LAI and photosynthetic capacity reach their maximum in early summer, at the same time.On the contrary, maximum photosynthetic capacity may be reached in late summer, and varies across the canopy.Without taking this into account, GPP could be overestimated up to 40% [66].Notwithstanding, comparing model outputs with published works [67,68], these defects are common also for other PBMs.
Seasonal patterns showed that the model consistently misrepresented NEE in winter, suggesting that R eco still needs to be improved.Especially for DBF sites (e.g., DEHai), winter R eco was mostly driven by RH.RH was exponentially affected by soil temperatures and especially moisture [69], which are calculated by the model, and could be over-fluctuating in winter.Moreover, EC data are prone to random noise [70], whose relative impact on performance metrics may be relatively larger.
Interestingly, annual predictions suggested reasonably high performance of 3D-CMCC-PSM, despite these seasonal inconsistencies.This suggests that biases are usually consistent within a season, but have different signs across seasons (Figure 8), resulting in a compensating effect at coarser time scale.
It was not always possible to individually validate the different components of respiration.Since EC Reco was not measured but inferred, evaluation metrics should be interpreted as a general ability of 3D-CMCC-PSM in predicting it.Despite daily Reco being noisy, 3D-CMCC-PSM could reproduce respiration processes well enough, at least at monthly timescale.
NPP:GPP ranged between 0.37 and 0.62, and was consistent with literature [71][72][73].It was generally higher in DBF (0.49 in DEHai, 0.55 in ITCol), lower in EBF (0.38 in FRPue, 0.41 in ITCpz).These results matched with those of [71] who showed that the ratio between NPP and GPP (CUE) is generally 0.53, ranging from 0.23 to 0.83.CUE was relatively low in FIHyy, where Reco predictions were overestimated: therefore, poor predictability was probably ascribable to excesses in RA.
It was not always possible to individually validate the different components of respiration.Since EC Reco was not measured but inferred, evaluation metrics should be interpreted as a general ability of 3D-CMCC-PSM in predicting it.Despite daily Reco being noisy, 3D-CMCC-PSM could reproduce respiration processes well enough, at least at monthly timescale.
NPP:GPP ranged between 0.37 and 0.62, and was consistent with literature [71][72][73].It was generally higher in DBF (0.49 in DEHai, 0.55 in ITCol), lower in EBF (0.38 in FRPue, 0.41 in ITCpz).These results matched with those of [71] who showed that the ratio between NPP and GPP (CUE) is generally 0.53, ranging from 0.23 to 0.83.CUE was relatively low in FIHyy, where Reco predictions were overestimated: therefore, poor predictability was probably ascribable to excesses in RA.Ecosystem respiration was overestimated during summer, causing NEE systematic overestimation at each simulation, especially in ITCol.Except for FRPue and ITCol, RA grew about 1 gC m −2 day −1 higher than RH, suggesting that it may be the principal driver of biased summer Reco (Figure 9).This misbehavior may be related to the method used to estimate maintenance respiration, an exponential relationship between respiration, moisture, and temperature [61].
Lack of data to validate the SOC dynamics reduced the spectrum of speculations, which could be statistically analyzed.SOC didn't change its quantity in ten years; this result was consistent with the theoretical stability of the SOC, an indicator which rarely changes within 10 years if no strong disturbance events (e.g., land use changes) have occurred [74].Litter C was highly fluctuating within a year, but its quantity was stable if compared at the end of each year.This suggested that the model realistically represented litter turnover and decomposition, since residues were degraded into humus labile substances about within a year [75].Microbial biomass was highly variable, as expected.However, the magnitude of change was too broad throughout the simulations.These results may be related to the use of 5% as the initial active microbial biomass for each site, a value that may be far from the equilibrium for different soils.Moreover, tradeoffs within microbial growth and between the environmental conditions may be scarcely represented.As a matter of fact, 3D-CMCC-PSM simulated the soil as having the same physical-chemical structure throughout the profile.This implied that microbes could find the same amount of C, O 2 , and living space, with no depth limitation.
1 gC m −2 day −1 higher than RH, suggesting that it may be the principal driver of biased summer Reco (Figure 9).This misbehavior may be related to the method used to estimate maintenance respiration, an exponential relationship between respiration, moisture, and temperature [61].Lack of data to validate the SOC dynamics reduced the spectrum of speculations, which could be statistically analyzed.SOC didn't change its quantity in ten years; this result was consistent with the theoretical stability of the SOC, an indicator which rarely changes within 10 years if no strong disturbance events (e.g., land use changes) have occurred [74].Litter C was highly fluctuating within a year, but its quantity was stable if compared at the end of each year.This suggested that the model realistically represented litter turnover and decomposition, since residues were degraded into humus labile substances about within a year [75].Microbial biomass was highly variable, as expected.However, the magnitude of change was too broad throughout the simulations.These results may be related to the use of 5% as the initial active microbial biomass for each site, a value that may be far from the equilibrium for different soils.Moreover, tradeoffs within microbial growth and between the environmental conditions may be scarcely represented.As a matter of fact, 3D-CMCC-PSM simulated the soil as having the same physical-chemical structure throughout the profile.This implied that microbes could find the same amount of C, O2, and living space, with no depth limitation.

3D-CMCC-PSM Uncertainty in Estimating NEE
We analyzed 300 random parameterization-sets per site to quantify model assumptions and uncertainty.The model showed different behaviors in different sites, but expectedly consistent across species behaving in a similar functional way.These results may suggest that using functional traits combinations to provide physiological parameters, instead of fixed species-specific ones, may produce still reliable and more general predictions, particularly useful in case of larger spatial/temporal simulations [76][77][78].Using a species-level parameterization, in fact, may result in a too-fine "resolution" because: (1) it would require excessive computational resources and a finelydetailed parameterization, usually inaccessible on a broad scale [79]; and (2) the model's rationale in predicting forest structure is mainly driven by competition for resources.However, there are not explicit tradeoffs, positive interactions between different tree cohorts, or intra-specific trait

3D-CMCC-PSM Uncertainty in Estimating NEE
We analyzed 300 random parameterization-sets per site to quantify model assumptions and uncertainty.The model showed different behaviors in different sites, but expectedly consistent across species behaving in a similar functional way.These results may suggest that using functional traits combinations to provide physiological parameters, instead of fixed species-specific ones, may produce still reliable and more general predictions, particularly useful in case of larger spatial/temporal simulations [76][77][78].Using a species-level parameterization, in fact, may result in a too-fine "resolution" because: (1) it would require excessive computational resources and a finely-detailed parameterization, usually inaccessible on a broad scale [79]; and (2) the model's rationale in predicting forest structure is mainly driven by competition for resources.However, there are not explicit tradeoffs, positive interactions between different tree cohorts, or intra-specific trait variability, which are fundamental to forecast forest ecosystem structure on long-run simulations [80].Having fixed species-specific parameters throughout a century would potentially result that only a very reduced amount of species would dominate the different cohorts on landscape to regional scale.
According to Figures 6-8, strong uncertainties still reside in timing for the different phenological phases.The biggest source of uncertainty in deciduous stands was driven by the amount of degree days needed to begin the vegetative period.The use of a site-specific thermal sum (GDD) to activate vegetation period is widely used, but proven to be very site sensitive, and not very effective for a regional generalization [81].On the other hand, the processes triggering bud-burst timing are still partly unknown.Moreover, those models proposing a process-oriented promoter-inhibition rationale are generally over-complex and not prone to spatial generalization [82].A possible solution in this context is to use remotely-sensed data to train a latitudinal-explicit regression, constraining GDD estimation.
3D-CMCC-PSM showed high uncertainty also in catching the beginning of the senescence phase.The new phenological scheme didn't reduce such uncertainty, since it was still using a photoperiod threshold as the senescence phase trigger [37].Another strong source of uncertainty in summer GPP may be held by the over-simplicity of soil structure and thus of the soil water routine.
As shown in [37], EC data are prone to high uncertainty.We focused on NEE fluxes to reduce the uncertainty cascade related to NEE partitioning.The next natural step will be reframing the model with a hierarchical Bayesian fashion, to quantify error propagation and parameter uncertainty from the posterior distribution [83].
Daily R eco estimation was affected by the cascade of uncertainties related to the calculation of R A and heterotrophic respiration, calculated independently.R A routine may strongly be influenced by uncertainties in R M estimation, which often resulted in R A overestimation.The R M was in fact simulated by a set of empirical relations, which involve the use of a fixed non-acclimating Q10 factor, whose generality is known to be inaccurate [84].Moreover, the rationale of Ryan's R M calculation [85] is affected by uncertainty in estimating daily increment of N pools, generally estimated by forest ecosystem models as a fixed proportion of daily C increment.

Conclusions
Soil respiration has a key role in determining NEE in a deterministic fashion [86].In general, this work showed how the inclusion of a simplistic soil carbon routine allowed prediction of trends and variability of NEE across the most diffuse European forest ecosystems.Modifications in the phenology scheme produced slight improvements in predicting GPP.However, they were still limited by the correct estimation of bud burst timing, leaf senescence starting point and duration.The use of an optimized parameter-set improved the model's performance only for those sites where the bio-geophysical processes were correctly reproduced.As a matter of fact, we showed how Mediterranean terrestrial forests, which showed lacks in representing some biological and/or physical processes, performed significantly worse than the other modeled sites, regardless of the use of optimized parameters.
In conclusion, we think that 3D-CMCC-PSM can reliably estimate NEE and R eco dynamics in a forest ecosystem, especially in scaling up daily results to monthly NEE averages.We think that 3D-CMCC-PSM is a solid basis to further explore the effects of soil structure on carbon and Water dynamics, especially in Mediterranean systems, and be used as a tool for predicting forest growth and ecosystem services, and address questions related to future scenario forecasting.

Figure 3 .
Figure 3. Soil Carbon dynamics in 3D-CMCC-PSM.The three macro pools are highlighted by red boxes (dead C pool) and blue box (live C pool, i.e., microbial).Blue-filled boxes represent the processes simulated by the soil model.

Figure 3 .
Figure 3. Soil Carbon dynamics in 3D-CMCC-PSM.The three macro pools are highlighted by red boxes (dead C pool) and blue box (live C pool, i.e., microbial).Blue-filled boxes represent the processes simulated by the soil model.
dependent factors (expressed in days).Graphically, the equation represents a skewed function of the amount of NSC allowed to be used by the trees of the specific class; the faster the leaves reach maturity, the more daily specific allocation is allowed (Figure4).

Figure 4 .
Figure 4. Graphic representation of the C tradeoff function.The axes represent respectively bud burst days (BBT), vegetative days (Veg_day), and the fraction of total NSC invested in leaf development (.The shorter BBT, the higher the maximum NSC fraction.

Figure 4 .
Figure 4. Graphic representation of the C tradeoff function.The axes represent respectively bud burst days (BB T ), vegetative days (Veg_day), and the fraction of total NSC invested in leaf development ( ∂NSC ∂t ).The shorter BB T , the higher the maximum NSC fraction.
obs represents EC daily NEE, Y sim Modeled NEE for the same day, Y obs the average EC daily NEE over the train time series.The first part of the RHS of the equation represents the square of the Pearson Correlation coefficient (R), the second the Nash-Sutcliffe Efficiency index (NSE).

Forests 2017, 8
to Eddy Covariance data on long-term daily, monthly and annual averages, over the full series of testing years (~5 years).For the validation, we used 2003-2007 for DEHai, 2004-2008 for FRPue and ITCpz, 2004-2011 for FIHyy, 2005-2012 for ITCol, 2009-2011 for ITRen.Then we evaluated how the model performed in the different seasons aggregating values for months of the same season.

Figure 5 .
Figure 5.Taylor diagrams representing 3D-CMCC-PSM performance in (a) daily, (b) monthly and (c) annual NEE estimation for the test-set.ITCol and DEHai represent DBF (red and green dots); ITRen

Figure 5 .
Figure 5.Taylor diagrams representing 3D-CMCC-PSM performance in (a) daily, (b) monthly and (c) annual NEE estimation for the test-set.ITCol and DEHai represent DBF (red and green dots); ITRen and FIHyy ENF (blue and turquoise).ITCpz and FRPue represent EBF (yellow and magenta dots).The closest a simulation lied to the "Ref" point, the better 3D-CMCC-PSM represented NEE patterns.Xand Yaxes represent NEE standard deviation (SD): the closest to 1, the better the performance.Simulations with R ≥ 0.9, and difference in SD with EC NEE less than 0.5 5 gC m −2 day −1 showed very good performance.Simulations with 0.75 ≤ R < 0.9, and difference in SD with EC NEE between 0.5 and 1 gC m −2 day −1 showed good performance.Simulations with 0.35 ≤ R < 0.75, and difference in SD with EC NEE between 1 and 1.5 gC m −2 day −1 showed sufficiently good performance.

Figure 6 .
Figure 6.Model structure-related uncertainty in estimating NEE (gC m −2 day −1 ) per DoS (Day of Simulation) by a random choice of parameter values from prior distributions.Data represent 1year simulations from randomly extracted parameterization-sets.Average daily simulations (black lines) and standard deviation (grey areas).Red dotted lines represent daily NEE simulation for the optimized parameterization set (TableS2).First column represents DBF sites, the second ENF, the third EBF.

Figure 6 .
Figure 6.Model structure-related uncertainty in estimating NEE (gC m −2 day −1 ) per DoS (Day of Simulation) by a random choice of parameter values from prior distributions.Data represent 300 1-year simulations from randomly extracted parameterization-sets.Average daily simulations (black lines) and standard deviation (grey areas).Red dotted lines represent daily NEE simulation for the optimized parameterization set (TableS2).First column represents DBF sites, the second ENF, the third EBF.

Figure 7 .
Figure 7. Inter-annual variability (IAV) for the test time-series (gC m −2 month −1 ).Observed IAV in gray boxes, simulated IAV in orange.First column represents DBF sites, the second ENF, the third EBF.

Figure 8 .
Figure 8. Patterns in daily NEE (gC m −2 day −1 ) per DoY (Day of Year) calculated from test-sets on a site level.Observed EC average patterns (black dotted line) and standard deviation (gray area).Simulated average patterns (red dotted line) and standard deviation (orange area).First column represents DBF sites, the second ENF, the third EBF.

Figure 8 .
Figure 8. Patterns in daily NEE (gC m −2 day −1 ) per DoY (Day of Year) calculated from test-sets on a site level.Observed EC average patterns (black dotted line) and standard deviation (gray area).Simulated average patterns (red dotted line) and standard deviation (orange area).First column represents DBF sites, the second ENF, the third EBF.

Figure 9 .
Figure 9. Patterns in daily soil respiration (SR) heterotrophic and autotrophic components (gC m −2 day −1 ).Patterns per DoY (Day of Year) calculated from test-sets on a site level.Microbial respiration average patterns (black-dotted lines) and standard deviation (gray areas).Autotrophic root respiration average patterns (red-dotted lines) and standard deviation (orange areas).First column represents DBF sites, the second ENF, the third EBF.

Figure 9 .
Figure 9. Patterns in daily soil respiration (SR) heterotrophic and autotrophic components (gC m −2 day −1 ).Patterns per DoY (Day of Year) calculated from test-sets on a site level.Microbial respiration average patterns (black-dotted lines) and standard deviation (gray areas).Autotrophic root respiration average patterns (red-dotted lines) and standard deviation (orange areas).First column represents DBF sites, the second ENF, the third EBF.

Table 1 .
Site descriptions and stand initialization data.Mean temperature (T) and Precipitation (prec.) are annual averages collected at site.

Table 2 .
Description of the different phenological phases for deciduous and evergreen species used in 3D-CMCC PSM.

Table 3 .
Statistics used for Model's results validation against Eddy Covariance data.

Table 4 .
Daily and Monthly Validation statistics calculated on the test-set.As stated in Table 3, R and NSE are dimensionless; RMSE and MAB are gC m −2 day −1 .

Table 5 .
Seasonal validation statistics calculated on the test-sets and aggregated by ecosystem type.As stated in Tables3 and 4, R and NSE are dimensionless; RMSE and MAB are gC m −2 day −1 .MAM stands for March, April, and May.JJA stands for June, July, and August.SON stands for September, October, and November.DJF stands for December, January, and February.DBF stands for Deciduous Broadleaf Forests, EBF for Evergreen Broadleaf Forests, ENF for Evergreen Needle leaf Forests.
, R and NSE are dimensionless; RMSE is gC m −2 day −1 .Bold values represent best performing version.

Table 7 .
Daily and Monthly Validation statistics for R eco calculated on the test-set.As stated in Table3, R and NSE are dimensionless; RMSE and MAB are gC m −2 day −1 .