Next Article in Journal
Evaluating the Legacy Effects of the Historical Predatory Seed Harvesting on the Species Composition and Structure of the Mixed Korean Pine and Broadleaf Forest from a Landscape Perspective
Previous Article in Journal
Pipe Model Can Accurately Estimate Crown Biomass of Larch (Larix olgensis) Plantation Forest in Northeast China
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Calibrating a Process-Based Model to Enhance Robustness in Carbon Sequestration Simulations: The Case of Cedrus atlantica (Endl.) Manetti ex Carrière

1
Department for Innovation in Biological, Agri-Food and Forest Systems (DIBAF), University of Tuscia, 01100 Viterbo, Italy
2
Division Impacts on Agriculture, Forests and Ecosystem Services (IAFES), Fondazione Centro Euro-Mediterraneo sui Cambiamenti Climatici, 01100 Viterbo, Italy
3
Department of Forest Development, National School of Forest Engineers, Salé 11000, Morocco
4
Forest Modelling Lab., Institute for Agriculture and Forestry Systems in the Mediterranean, National Research Council of Italy (CNR-ISAFOM), Via Madonna Alta 128, 06128 Perugia, Italy
5
Department of Agriculture and Forest Sciences (DAFNAE), University of Tuscia, 01100 Viterbo, Italy
6
National Biodiversity Future Center (NBFC), 90133 Palermo, Italy
7
BioBio Research Center, Faculty of Sciences, Mohammed V University, Rabat 10000, Morocco
*
Authors to whom correspondence should be addressed.
Forests 2023, 14(2), 401; https://doi.org/10.3390/f14020401
Submission received: 18 December 2022 / Revised: 8 February 2023 / Accepted: 12 February 2023 / Published: 16 February 2023
(This article belongs to the Section Forest Ecology and Management)

Abstract

:
To assess the degree to which it has met its commitments under the Paris Agreement, Morocco is called upon to carry out carbon assessments and transparent evaluations. Within the forestry sector, little is known about the role of Morocco’s forests in contributing to carbon uptake. With this aim, we applied for the first time in the literature the 3-PG model to Cedrus atlantica ((Endl.) Manetti ex Carrière, 1855), which represents about 131,800 ha of Morocco’s forest area (i.e., Azrou forest). Through the Differential Evolution-Markov Chains (DE-MC) we tested and assessed the sensitivity and calibrated the 3-PG model. This process-based model provided significant results regarding the carbon sequestration capacity. The results showed the following: i. Parameters related to stand properties, canopy structure, and processes, as well as biomass partitioning, are the most important or sensitive for the performance of the model; ii. The DE-MC method optimized the values of the 3-PG parameters which was confirmed by the means of the Gelman–Rubin convergence test; iii. According to the predictions of the calibrated 3-PG, the Net Primary Production in the pure Azrou forest varies between 0.35 and 8.82 tC . ha 1 . yr 1 , it is equal in average to 5.48 tC . ha 1 . yr 1 , which given the total area corresponds to 7918 tC . ha 1 .

1. Introduction

Forest ecosystems have the capacity to produce goods and services that are essential for human well-being [1,2]. These services range from providing wood products and drinking water to cultural and supporting services [3]. Forests have been affected by humans more rapidly than ever in recent years [4]. Climate change brings novel severity and timing of multiple stresses, which may significantly affect the provision of these Ecosystem Services (ES) [5,6]. Consequently, these changes introduce considerable uncertainty into the management of forest resources in order to sustain the provision of ES [7].
The Intergovernmental Panel on Climate Change (IPCC) reported a requirement for large-scale carbon dioxide (CO2) removal in most scenarios that limit global warming to 1.5 °C or well below 2 °C [8]. Net-zero emissions and global decarbonization targets recognize the important role of forest ecosystems in contributing to these ambitions [8,9]. Forests currently absorb around 30% of CO2 emissions each year and play a central role in the terrestrial ecosystem and carbon cycles [10].
Methods for carbon storage assessment can be divided into four main categories: i. inventory-based estimation which is a collection of methods generally applied to estimate carbon storage on the basis of regional forest inventory data; ii. Satellite-based estimation which makes use of active, or passive remote sensing or both in order to quantify aboveground carbon; iii. estimation through modelling which can rely on simulating the main ecophysiological processes governing forest evolution; and, iv. the combination of the three aforementioned categories [11].
The complexity of forest ecosystems and the multiple interactions between different compartments, impose major difficulties for researchers to predict accurately their systemic responses based on simple statistical relationships. Constantly updated tools to describe vegetation dynamics in an evolving context of climate change have emerged, including process-based models (PBMs) [12]. These models are built around the main eco-physiological processes underlying productivity and can include detailed representations of competition, population dynamics, and forest succession [13,14].
PBMs can perform better than the empirical growth models and showed to be more accurate for long-term projections [15,16,17,18]. Among the most widely used PBMs in forestry literature, the Physiological Principles Predicting Growth (3-PG) model [19], uses simple equations to simulate physiological processes controlling carbon balance such as solar radiation absorption, carbon sequestration and allocation, mortality, as well as water balance such as canopy interception and evapotranspiration [20]. In addition, 3-PG is parsimonious in the requirement of parameters to calibrate if compared to other models from the same class [21], and it incorporates modules that convert biological outputs, based on allometric equations, into operational data of immediate interest to forest managers such as diameter at breast height (DBH), tree height and biomass stocks which makes it appropriate for forest management purposes [13]. This model has been successfully tested in the last twenty-five years in many regions of the world and provided interesting results for a variety of tree species: Cunningamia lanceolata (Lamb.) Hook [22], Eucalyptus grandis W. Hill and Eucalyptus urophylla S.T.Blake [23], Eucalyptus grandis W.Hill and Eucalyptus camaldulensis Dehnham [20,24,25], Eucalyptus globulus Labill. [15,25], Pinus ponderosa Douglas ex C.Lawson [26,27], Pinus patula Seem. [28], Pinus sylvestris L. [29,30], Picea sitchensis (Bong.) Carrière [31] and Mediterranean maquis (e.g., Quercus ilex L.) [32] but, to our knowledge, never for Cedrus atlantica ((Endl) Manetti ex Carrière, 1855), which represents about 131,800 ha of Morocco’s forest area.
Morocco is among the most committed countries on climate and sustainable development issues. As a ratifier of the Paris agreement, this country has agreed to periodically present its Nationally Determined Contributions (NDCs) while being consistent with the enhanced transparency framework (ETF). Concerning the forestry sector, the lack of growth models for different species and growth conditions limits the tier assessment approach adopted. PBMs, if properly calibrated, could be interesting prospects and valuable tools, and this is the purpose of the present work.
Specifically, in the present work we will: i. try to define the most sensitive parameters which are going to describe the evolution of Cedrus atlantica at the level of the Azrou forest (Morocco); ii. try to optimize the 3-PG model parameters as to provide, for the first time in the literature, a calibration set for Cedrus atlantica, and iii. we will simulate the quantity of carbon sequestered in terms of Net Primary Productivity (NPP) by the cedar forest of Azrou during the period 2016–2021 using the calibrated model.

2. Materials and Methods

2.1. Study Area

The study was conducted in Azrou forest, Morocco. It is a large forest massif covering 178 km2 in the central part of the Middle Atlas (between 5.00° to 5.29° W and 33.28° to 33.52° N) (Figure 1).
The Middle Atlas is famous for its magnificent Atlas cedar and holds about 70.4% of the total area of this species in Morocco which corresponds to 93,500 ha [33]. Cedrus atlantica (Endl.) Manetti ex Carrière, endemic species of Morocco and Algeria [34], occurs at elevations between 1500 and 2400 m a.s.l [35]. Middle Atlas cedar forests contain several deciduous and evergreen tree species (e.g., Quercus rotundifolia, Quercus canariensis, Pinus pinaster, Ilex aquifolium), shrub species, short meadowlands, and other aromatic and medicinal plants.
Cedrus atlantica, which represents the main species of Azrou forest, forms pure or mixed stands with Quercus rotundifolia, Quercus canariensis, and secondary species depending on the nature of the substrate. The Cedrus atlantica stands occupy 1491.41 ha in a pure state, 7182 ha in mixture with Quercus canariensis representing 48.74% of the total area of the forest; the Quercus canariensis stands extend over an area of 4419.77 ha representing 25% of the forest [36]. In the present work, only pure cedar stands will be considered and discussed. Figure 1 illustrates the spatial extent of these stands.
Cedar ecosystem fulfils important functions and provides several services for the local population and human well-being (e.g., recreation, cultural inspiration, habitat for wildlife, place of grazing, Carbon and water regulation). However, these ecosystems suffer from recent climatic variations [37] which are amplified by anthropogenic pressures and disturbances (e.g., overgrazing, excessive clearance of woodlands, overexploitation) [38]. Besides the anthropogenic disturbance, an increase in drought occurrence and severity has been documented over a large part of the Azrou forest. This established fact was a predisposing factor in the dieback of the cedar, particularly on sites characterized by the absence of silviculture supposed to adapt stand density and structure to soil water availability [39]. In recent years, disturbance regimes in many forest ecosystems have profoundly changed with the climate being an important driver of this dynamic. Thus, given the causal relationship between climate and drought with an expected decrease in precipitations at the level of the Mediterranean region, drought events are projected to increase in frequencies, sizes, and severities within the Azrou forest [40].

2.2. Overall Approach

In order to deal with the calibration of 3-PG for C. atlantica, various data has been collected which was used to analyze model sensitivity and then to calibrate the model. Figure 2 shows the flowchart of the overall methodology used here. Data was collected from various sources: i. field data (forest resources); ii. data extracted from remotely sensed earth observations (GEE database); and, iii. data collected from other databases (the European Soil Database) [41]. Sensitivity analysis was conducted to define the most important parameters for the 3-PG model. The 20 most sensitive parameters out of 55 which represent the total of parameters in 3-PG were retained from sensitivity analysis and used for optimization according to the Differential Evolution-Markov Chains (DE-MC) approach described below. The calibrated model was then used to predict the NPP over the entire study area and compared to the observed MODIS NPP of reference for the period 2016–2021.

2.3. Used Data

The 3-PG model requires a set of input data for its initialization which includes: (1) model species parameters (species-specific eco-physiological and allometric characteristics that can be partially derived from forest inventories and literature), (2) site data (e.g., latitude, altitude, information about physics of the soil), (3) initial stand structural data (e.g., DBH, tree height, stand density), (4) climatic data (meteorological data at a monthly time step, e.g., mean daily incident solar radiation, mean monthly air temperature). Such data could be classified as in situ and ex-situ data according to their way of collection. In addition, the calibration process of this model requires observational data (5) for model validation.

2.3.1. In-Situ Data

As tree growth depends on stand density, and climate conditions, the study area was stratified according to stand density, exposure, and soil type. Due to the absence of permanent plots, the sampling design was based on temporary circular plots of 0.05 ha with their borders set by means of a Haglöf Vertex IV and its transporter including any tree whose distance from the plot center is less than or equal to 12.6 m. In total, 24 plots were randomly selected within 4 sampling strata according to an optimal allocation (considering stratification criteria: stands density and exposure). The variance of the estimate was reduced according to the Cochran formula [42] and the relative sampling error (RSE), which was found to be equal to 14.5%, was computed according to the equation presented below (Equation (1)).
R S E ( Y ¯ ) = 1 Y ¯ 1 N 2 i = 1 k ( N t ( N t n t ) s t 2 n t )
with Y ¯ representing the estimate of the population mean annual carbon increment, N and N t correspond to the total number of units within the population and within each stratum, respectively. n t and s t represent the number of sampling units and the standard deviation of the annual increase in carbomass, respectively, in each stratum from the k ensemble strata (Table 1).
Due to the large scale of the soil map, substrate representativeness was verified during the fieldwork. Indeed, 9 plots are based on a limestone bedrock while 15 plots are based on a basaltic bedrock, which implies the representativeness of samples with regard to the type of bedrock.
To derive the stand’s current and previous (2016, the start year of the simulations) characteristics, plots were inventoried. Data related to the site and tree parameters were collected. At each sample plot, the stem circumferences of all trees with a circumference greater than the pre-countable circumference of 20 cm were measured using a Tricle tape measure. Setting such a threshold for the tree measurement is justified by the fact that: i- biomass equations used here were developed on the basis of trees of medium to large dimensions, thus it is statistically inappropriate to apply them to small dimensions, ii- the sampling plots in our study have not been subject to recent disturbances which imply that the proportion of small trees in those plots is too low to impact the overall plot biomass. Furthermore, among the trees in the plot, three representative trees were measured for height (using a Haglöf Vertex IV), and their radial increments over the last five years were measured using a Haglöf Pressler corer used for core extraction and a BORLETTI calliper for measuring the thickness of the last five rings.

2.3.2. Ex-Situ Data

Other categories of data, specifically (2), (4), and (6), were derived using remotely sensed data using the cloud computing platform Google Earth Engine (GEE) [43]. Data for the available soil water variable from site parameters category (2) was extracted from the European Soil Database [41]. Below is a summary table of the data used in this study along with their metadata (Table 2).

2.4. Data Preparation and Preprocessing

2.4.1. Model Initialization

To assess the initial biomass carbon stock, the use of carbomass models developed at the level of Azrou forest [44] was relevant instead of taking into account the default parameters of the IPCC. These models with an associated error (RMSE) varying between 0.9 and 1.65 according to the tree compartment, were used to estimate the amount of carbon per compartment (i.e., aboveground part, stem, and foliage) per tree on each plot (Table 3). The conversion factors that emerged from the previous study were used to ensure the transition from biomass to carbomass in each tree compartment (Table 4). However, since root biomass was not the subject of the aforementioned study, the ratio of root biomass to aboveground biomass (0.29) proposed in [45] was used in the present work.
The average diameter increment was computed at each plot level and was used to generate the 2016 biomass stock for each compartment namely: stem biomass, root biomass, and foliage biomass.

2.4.2. Climate Data

The data preprocessing was conducted through the GEE platform using the Web programming interface to ensure work efficiency [43]. Image collections were loaded in the working environment, and then they were filtered by the area of interest and further by date to include only the images related to the study period. Then, all images were multiplied by their corresponding scale factor. Some collections did not require any additional preprocessing, solar radiation (srad) whose temporal resolution equals the simulation step of the 3-PG model (1 month), while the other collections required an additive preprocessing to bring their temporal resolution to one month.

2.5. Sensitivity Analysis

Accurate predictions and simulations depend on the accuracy of model parameter setting, climate data, and site information [46,47]. Sensitivity analysis defines the sensitivity as well as the importance of each model parameter and provides a sufficient basis for selection during model calibration. Optimizing parameters with low sensitivity increases the computation time without significantly improving the accuracy of the model [48].
The revised Morris method was chosen [49,50]. It is a one-step-at-time (OAT) method, which means that for each iteration of the algorithm, only one parameter assigns itself a new value, and the others remain unchanged [51,52]. The reason for choosing this approach is that it is well suited for models containing dozens of factors without depending on strict assumptions about the model such as additivity and monotonicity of the input-output relationship of the model. Moreover, this method is easy to understand and implement, and its results are easy to interpret [53].
The Morris algorithm starts at a randomly chosen in a k- dimensional space, where k represents the number of model parameters for which the sensitivity will be evaluated, and creates a trajectory through the k- dimensional space. Two neighbouring points differ by a distance called a “step” in only one direction. For each trajectory, a given variable can be assigned a discrete number of values called “levels” sampled evenly within the factor range of variation. This process of trajectory construction is iterated a number of times in such a way as to explore all the possible factor levels. Then, for each “repetition” and parameter, the difference in the response of the target variable considering two adjacent points is calculated, which conventionally represents the elementary effect. Based on this difference, it is possible to compute for each parameter two metrics namely μ * which refers to the mean of the absolute elementary effect, and σ which refers to the standard deviation of the elementary effect. Thereafter, considering the different level values of those metrics, model parameters could be classified into three main groups, namely parameters with a negligible effect, parameters with linear effect without interaction, and parameters with a non-linear or interaction effect, respectively.
The parameters of the Morris algorithm have been assigned the following values: (Number of levels: 20, number of repetitions: 500, step: 3). For the output variable of the Morris method, which is used for the computation of the elementary effect, it was assimilated to the fit to the observed variable (NPP observed) and expressed by its log-likelihood with normal error assumption for NPP.
The choice of the range of variation for the parameters of the present model was based on literature analysis. In the absence of a reference value specific to the genus Cedrus sp., a default value general for coniferous species was used, and the minimum and maximum values of the variation intervals were set given those considerations [17]. The sensitivity analysis was performed on all of the 55 parameters of the 3-PG model [13].

2.6. Model Calibration

The DE-MC (Differential Evolution-Markov Chains) method, which is a combination of Markov Chain Monte Carlo (MCMC) and genetic algorithms (GA), has been used for the calibration and optimization of the 3-PG model parameters. This method that combines a priori knowledge about parameters with observations is widely used in ecological research [54,55,56]. The a posteriori values of the parameters can be used as a result of the calibration, and the optimized parameter set of the model can be compared to the initial parametrization set, which seems to be extremely useful in identifying parameters that estimate the observed evolution of a given species.
The Bayesian theory combines a priori knowledge about the parameters of the model with observations of the variables that will be predicted by the model to perform a posterior estimation of these parameters [57]. The MCMC method involves the construction of a Markov chain with the parameters of the posterior distribution to obtain posterior samples of these parameters and subsequently infer the numerical characteristics of these parameters based on these samples. Bayesian theory is expressed by the following formula:
P ( θ / y ) = f ( y / θ ) g ( θ ) f ( y / θ ) g ( θ ) d ( θ )
with θ and y represent the parameters and output values simulated by the 3-PG model (e.g., NPP), respectively; P( θ /y) is the posterior probability density function of the parameters, and f(y/ θ ) refers to the observations. The conditional probability density knowing the parameters a priori is called the likelihood function, g( θ ) being the a priori distribution of the parameters, and d( θ ) corresponds to the differential of the model parameters. To solve the scaling and orientation problem of the jumping distribution, Braak [58] proposed a method combining MCMC and DE which inherits GA thereby giving DE-MC. For the DE-MC case, which is defined as the distance between the current parameters vector and the candidate vector, is calculated in DE-MC as a multiple of the difference between two parameter vectors of the current population. The selection process of DE-MC is based on the Metropolis ratio defining the probability for which a candidate could be successful [59]. Based on the results of the sensitivity analysis, the 20 most influential parameters were selected to reduce the computational time and improve the efficiency of the optimization as in Trotsiuk et al. [60]. Bayesian calibration was performed using differential evolution [58], and the MCMC algorithm from the BayesianTools library [61] with 3 chains and 4 × 10 6 iterations (see Figure 3).

3. Results

3.1. Stands Characterization

3.1.1. Adjustment of the Height-Circumference Relationship

Considering that the estimation of biomass for Cedrus altlantica is based on a two inputs model (stem circumference at 1.30 m (C in meters) and tree height (H in meters)) and that the only parameter that was measured for all the tree individuals is the circumference, (see Methods), the development of Height-Circumference model was necessary. For this specific purpose, a second-degree polynomial model was fitted using the ordinary least square method (OLS). As a result of the fitting, we obtained the following Equation ( R 2 = 0.80, p 0.0001 ).
H = 1.052 + 17.548 C 2.606 C 2

3.1.2. Results of Stands’ Characterization

Based on the aggregation of data collected at the sampling plot level, it is possible to characterize the forest by homogeneous strata. Table 5, presents by stratum level, the average of all dendrometric descriptors collected in the field and biomass in different tree compartments presented as tonnes of dry matter per hectare (tDM. ha 1 ).
Results show that low-density stands had the highest increment in circumference, 1.96 and 1.13 cm . yr 1 , respectively, for warm and cool exposures, while the lowest values were observed in highly dense stands, 1.1 and 0.93 cm . yr 1 , respectively, for cool and warm exposures. The results obtained, concerning the density of the stand and the biomass of different tree compartments, were used to initialize the model.

3.2. Sensitivity Results

The results from the Morris sensitivity analysis highlight that the parameters related to stand properties, canopy structure and physiological processes, as well as biomass partitioning are particularly the most important or sensitive for the model’s performance More explicitly, the parameters with an important overall influence on the NPP, having the highest value of μ *, were in order: alphaCx (Canopy quantum efficiency), fN0 (Value of fN fertility ratio (FR) equals to 0), rAge (relative age to give an age growth-modifier value (fAge) of 0.5), MaxCond (Maximum canopy conductance), fNn (power in the fertility equation), Topt (Optimum temperature for growth), SLA1 (Specific leaf area for mature leaves), pFS20 (Foliage-Stem partitioning ratio at B = 20 cm). It was also found that those same parameters tend to produce larger σ , which could indicate non-linearities or interactions with other parameters (Figure 4).

3.3. Calibration

Based on the Morris analysis, the 20 most influential parameters were used for calibration purposes using the Differential Evolution Markov Chain Monte Carlo algorithm (DEzs MCMC) [62]. The calibration of the model converged to parameter values close to the values used to generate the reference data. This was verified by means of the Gelman–Rubin convergence test in which the potential scale reduction factor (psrf) is equal to 1.01. The value of psrf is lower than 1.1, which confirms the convergence of the calibration [63,64]. It also emerges from this study that the posterior distribution is well-defined. In the bayesian calibration, the posterior probability for the model parameters should have a prominent peak without being bimodal. This is generally achieved with a high number of iterations (e.g., 10 6 in the case of our study) to narrow the parameter space. The normal distribution or pseudo-normal (skewed) validates all those criteria (Figure 5).
Furthermore, to evaluate the performance of the model, the predictive posterior distribution was calculated by running the model with 500 random samples from the parameters’ global posterior distribution. Then, the NPP was simulated over time for each combination of parameters. The best parameter set was then assimilated to the maximum posterior probability which corresponds to the mode of the posterior probability distribution. To verify the superiority of our parametrization dataset, the comparison of the results provided by the three parametrization sets, namely the default, the median, and the calibrated set of parameters was done using the mean squared error metric (MSE) with reference to the MODIS observations on NPP ( N P P M O D I S ) with n being the number of observation on NPP (Equation (7)). The results show that the proposed parameterization dataset is indeed more accurate than the default one and that MSE was reduced by about 18% (Table 6).
M S E = 1 n i = 1 n ( N P P p r e d i c t e d N P P M O D I S ) 2
Given these results, the obtained parametrization dataset will be used for future estimates of carbon sequestration by the Azrou cedar forest. Below is an explicit description of the 3-PG model parameters optimized in the present study (Table 7).

3.4. NPP Simulation

In order to provide an overall estimate of carbon sequestration capacity by the pure cedar forest of Azrou (we ran the 3-PG for the period 2016–2021), the estimation of sequestration by each sampling unit was performed using the calibrated 3-PG model. First, a simulation of NPP was carried out by sample plot, then these results were aggregated by stratum and their respective averages were calculated. Subsequently, in order to provide an overview of the carbon sequestered in the Azrou forest, the average annual amount of sequestrated carbon per unit area was multiplied by the area occupied by each stratum. Later, a conversion factor of 0.56 was used to convert biomass (dry weight) to carbon equivalent (C) (tonnes). These results are represented in the table below (Table 8).
Results show that dense stands have an important carbon uptake dynamic in comparison to stands with lower densities. It also emerges that for stands belonging to the same density category, cool exposures always sequestrate annually more carbon than warm exposures. In fact, for the high-density category, stands with cool exposure annually sequestrate 7.24 tC . ha 1 . yr 1 while those with warm exposure sequestrate 5.92 tC . ha 1 . yr 1 . Regarding low-density stands, those with cool exposure annually sequestrate 5.20 tC . ha 1 . yr 1 , while those with warm exposure annually sequestrate 5.08 tC . ha 1 . yr 1 .
Moreover, although it has been demonstrated in the present work that stands with cool exposure hold the highest sequestrating capacity, warm exposed stands within the forest of Azrou contribute to 67% of the annual carbon uptake while cool exposed stand contributes only to 33% of the annual carbon uptake. However, stands with warm exposure occupy 72% of the total area of the forest, while those with cool exposure occupy only 28% of the total area of the forest.

4. Discussion

4.1. Parameter Sensitivity and Optimization

In the present study, the overall sensitivity of 55 parameters of the 3-PG model in simulating NPP was analyzed using Morris’ sensitivity. As a result, it was found that NPP simulated by the 3-PG model is particularly sensitive to some parameters related to canopy structure and processes (e.g., alphaCx, SLA1, MaxCond), and conductance modifiers (e.g., fN0, rAge, fNn, Topt), as well as carbon partitioning (e.g., pFS20). While there are three possible approaches to assigning values to parameters in a model namely: i. direct measurement, ii analogy with other species and iii. parameter estimation [13]. The present work exclusively used the last two approaches given the unavailability of direct measurement for this species. It is in that sense that the results of the sensitivity analysis were leveraged and allowed to assign generic values derived from a benchmarking that included conifer species to parameters with a low ranking in sensitivity, while 20 parameters with a high ranking in sensitivity were fitted to the observations on NPP. To illustrate the importance of sensitivity analysis in understanding the behaviour of the model, some examples will be discussed in the following.

4.1.1. NPP and Canopy Processes and Structure

The maximum quantum efficiency (alphaCx; i.e., the maximum attainable efficiency when no environmental or structural modifiers limit the maximum potential photosynthesis) is directly used for the calculation of NPP. NPP results from the multiplication of ‘alphaCx’ by the absorbed photosynthetically active radiation (APAR), canopy cover, and a series of environmental and structural modifiers for the computation of the Gross Primary productivity (GPP) and then reduced through the NPP/GPP ratio (Y; i.e., NPP = GPP ∗ Y) [72,73] which represents one of the twenty most influential parameters. Similarly, other models (e.g., 3D-CMCC-FEM and CLM4.5-FATES) which simulate NPP mechanistically as the net result from GPP less Autotrophic Respiration [52,74] have shown that the number of live biomass controls mostly NPP [75]. In addition, the non-trivial role of Y is something that has been discussed by Collalti et al. [76], and represents a large source of uncertainty for models who apply the “NPP = GPP ∗ Y” given that it has been found to range from 0.22 to 0.79 [72]. GPP (and then NPP) was also found to be sensitive to ‘SLA1’, this parameter refers to the specific leaf area for mature trees which varies between stands of different ages and is necessary for the computation of Leaf Area Index (LAI) and GPP, subsequently. GPP was also found to be sensitive to the maximum canopy conductance (MaxCond), given that canopy conductance is central to calculations of canopy photosynthesis and thus NPP [77]. Similarly to our results, Almeida et al. [20] found that biomass-related outputs for Eucalyptus grandis W.Hill (root, stem, and foliage biomass) were sensitive in 3-PG to alphaCx and Maxcond. In Keryn I. et al. [78], NPP has been found to be sensitive for Eucalyptus grandis W.Hill and Pinus radiata D.Don to alphaCx and Maxcond. Those findings agree with our results showing that independently from species considered, the NPP simulated by the 3-PG model is sensitive mostly to these parameters. Moreover, results obtained from different studies using various PBMs were in good agreement with our findings. Typically, in a study conducted by Zaehle et al. [79], using the LPJ-DGVM model, both the intrinsic quantum efficiency and maximum canopy conductance were among the four most important parameters controlling NPP. In Pappas et al. [80], intrinsic quantum efficiency was of utmost importance in LPJ-GUESS parameterization, explaining most of the variability in carbon fluxes. In a different study led by Tatarinov et al. [81], using the BIOME-BGC model, the effect of ’SLA’ on NPP was strong for both Fagus sylvatica L. and Picea abies L.

4.1.2. NPP and Conductance Modifiers

In addition to parameters of the canopy structure and processes, ‘fN0’ and ‘fNn’, representing the coefficients of the relationship between the fertility index and the fertility growth modifier were found to be important for the model, which is consistent with other studies that demonstrated the relevance of the relation between the fertility indexes and the stand productivity [82]. Moreover, ‘rage’ was also found to be important for the performance of the model, this parameter intervenes in the definition of the age-related growth modifier which in turn modifies canopy conductance and quantum efficiency which are accounted for in the calculations of GPP and NPP [20]. The optimum temperature for growth (Topt) had also a high sensitivity ranking for GPP and NPP, this could be attributed to its contribution to the calculation of the effective quantum canopy efficiency and thus to the NPP [19]. In Trotsiuk et al. [60], which applied 3-PG to examine the growth of Pinus sylvestris L. and Fagus sylvatica L. mixtures along site and climatic gradients, biomass-related outputs (Wr, Wf, and Ws) were found to be sensitive to fN0, fNn, and rAge, which is consistent with our findings. In addition, our results were in good agreement with Xenakis et al. [30], who applied 3-PG on commercial plantations of Pinus sylvestris L., found that fN0 was important for foliage and root biomass while both fN0 and Topt were important for stem biomass.

4.1.3. NPP and Carbon Partitioning

In addition to parameters belonging to the canopy structure and processes as well as the conductance modifiers classes, ‘pFS20’, which refers to the partitioning ratio between foliage and stem for a representative tree of DBH = 20 cm, has been found to be important for the performance of the model. The high sensitivity of NPP to this parameter could be explained by the fact that pFS20 is a determining factor in the updating of the foliage biomass pool and thus of the Leaf Area Index (LAI) which intervenes in the calculation of the NPP through the absorbed active radiations. In a study conducted by Ulrich et al. [83], who applied the 3-PG model on Pinus ponderosa L., GPP was found to be sensitive to ‘PFS20’. In fact, it was shown that an increase of about 40% in the value of ‘PFS20’ induces an increase of about 22.5% in GPP, thus in NPP as well given the fixed proportionality between GPP and NPP. Results from other studies, e.g., the one by Massoud et al. [74] that made use of the CLM4.5(FATES) model, show the importance of leaf and stem allometry parameters controlling dynamic carbon allocation and thus the general vegetative state and size structure of the forest.

4.2. Implications, Future Perspectives, and Limitations

Results show that the calibrated model decreased the error by about 18% compared to the default parametrization dataset. The 3-PG model was calibrated using time series, and different sites since NPP is inherently dynamic in contrast to slowly varying state variables (e.g., DBH). To our knowledge, the present study is the only one in the literature to have calibrated the 3-PG model for the Atlas cedar species. The present calibration set for 3-PG will enable a better assessment of the carbon sequestration by a cedar forest instead of using the default increment as suggested by the IPCC. The calibrated model also offers the possibility of predicting the impact of climate change on forest productivity under different management options, which is of paramount importance in the context of forest disturbances being more pronounced in recent years and likely to be of the same magnitude in the coming decades [84]. Specifically, in the context of our study area, drought which is defined as a sporadic disturbance of the water cycle, is becoming more frequent in recent years. As a response, plants deploy water-use strategies to avoid excessive water consumption which are mainly translated into physiological and structural changes. Physiological responses include a decrease in stomatal conductance through stomatal closure to control the water flux. This may also induce a reduction in CO2 diffusion and thus a short-term reduction in GPP. From a structural perspective, for example, a reduction in leaf area due to early senescence could be responsible for the reduction of photosynthetic activity (GPP). Besides the carbon fluxes change, drought could also induce vegetation mortality, e.g., due to cavitation when the demand for water exceeds the supply, and the air is aspirated into the xylem or carbon starvation when autotrophic respiration exceeds photosynthesis. However, from a technical perspective, 3-PG is only able to represent the effect of drought on GPP, while density-independent mortality in general and drought-induced mortality in specific is not conceived as part of the 3-PG scheme.
The model calibrated in this study simulates the evolution of NPP over time in a pure forest of Cedrus atlantica. Forest stand dynamics are not only dependent on the ecological and physiological characteristics of the species but also on the human interventions that may prevail by station [57]. This last element has not been considered due to the unavailability of permanent plots with recorded monitoring, however, this lack has been compensated by a choice of temporary plots that were recently intact. Ideally, a network of monitoring plots should be set up in the Azrou forest to properly follow the evolution of its stands. The data thus obtained could help calibrate the 3-PG model considering its various components. Additionally, given the importance of the quality of observations for the performance of the model, using the method of Eddy Covariance could be a prominent tool to meet the quality standards enabling the monitoring of GPP with relatively a high certainty [85] instead of using NPP of MODIS having a spatial resolution of 500 m, generally exceeding the plot size, we can lead to spatial error in the case of very heterogeneous forests. It should also be mentioned that in the Azrou forest, the Atlas cedar is found either in a pure state, mixed with the holm oak or to a less degree with the Mirbeck’s oak. Owing to this species mixture, the evaluation of carbon sequestration in the whole forest of Azrou using 3-PG requires the calibration of the holm oak model for which the literature is quite rich [32,86].
From another perspective, it should be noted that the amount of carbon in soil represents a substantial portion of the carbon found in terrestrial ecosystems. In a recent study conducted at the level of the Azrou forest, soil carbon has been found to constitute between 46.4% and 59.1% of the total carbon stock [87]. The 3-PG model used for the frame of our study is only able to follow the evolution of carbon in the biomass pools, however, an upgraded version of this same model (i.e., 3PGN) couples 3PG with ICBM/2N which is a soil model and provide a way to account for the soil carbon balance [30]. Moreover, in the context of the present study, the root-to-shoot ratio provided by the IPCC has been used given the non-availability of accurate vegetation-specific values. For the purpose of better accuracy in the estimation of biomass, it will be advantageous to calibrate/develop specific allometric equations in the future.
The results obtained in this work show that the annual carbon update in the pure cedar forest of Azrou varies from 0.35 and 8.82 tC . ha 1 . yr 1 (inter-plot variability). In compliance with our results, previous studies have demonstrated the ability of the Atlas cedar to achieve high levels of the annual increment in volume ranging from 8 to 12 m 3 . ha 1 . yr 1 [88,89]. However, the comparison of these findings with the generic increment suggested by the IPCC, and adopted by the National Forest Inventory services (NFIs) which is about 0.47 tC . ha 1 . yr 1 highly underestimates Cedar forest carbon sequestration capacity and should raise questions about the relevance of the methodology used.

5. Conclusions

In this study, the sensitivity of 55 parameters of the 3-PG to the NPP was analyzed during the period 2016–2021. The 20 parameters with a high influence on NPP have been selected and optimized according to the DE-MC method. The conclusions of this study are as follows:
(1)
Parameters related to stand properties, canopy structure, and processes, as well as biomass partitioning, are the most important or sensitive for the performance of the model.
(2)
DE-MC method optimized the values of the 3-PG parameters which was confirmed by the means of Gelman–Rubin convergence test.
(3)
According to the predictions of 3-PG, the annual carbon sequestration in the pure Azrou forest varies between 0.35 and 8.82 tC . ha 1 . yr 1 , it is equal in average to 5.48 tC . ha 1 . yr 1 , which given the total area corresponds to 7918 tC . yr 1 .

Author Contributions

Conceptualization, S.L. and I.B.; methodology, S.L. and I.B.; resources, S.L. and R.V.; writing—original draft preparation, S.L., I.B., S.M. and A.C. All authors contributed to the writing and reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

A zip file containing the data and scripts used in the frame of the present study can be requested from [email protected] or [email protected].

Acknowledgments

We thank Anass Legdou, and through him the Moroccan Department of Waters and Forests (Regional Directorate of the Middle Atlas) and the National School of Forest Engineers (Salé, Morocco) for providing logistical and intellectual support. We would like also to thank the reviewers for taking the necessary time and effort to review the manuscript. We sincerely appreciate all your valuable comments and suggestions, which helped us in improving the quality of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Wang, Z.; Zhang, B.; Zhang, S.; Li, X.; Liu, D.; Song, K.; Li, J.; Li, F.; Duan, H. Changes of Land Use and of Ecosystem Service Values in Sanjiang Plain, Northeast China. Environ. Monit. Assess. 2006, 112, 69–91. [Google Scholar] [CrossRef] [PubMed]
  2. Costanza, R.; d’Arge, R.; de Groot, R.; Farber, S.; Grasso, M.; Hannon, B.; Limburg, K.; Naeem, S.; O’Neill, R.V.; Paruelo, J.; et al. The Value of the World’s Ecosystem Services and Natural Capital. Nature 1997, 387, 253–260. [Google Scholar] [CrossRef]
  3. Fabrika, M.; Valent, P.; Merganičová, K. Forest Modelling and Visualisation—State of the Art and Perspectives. Cent. Eur. For. J. 2019, 65, 147–165. [Google Scholar] [CrossRef]
  4. d’Annunzio, R.; Sandker, M.; Finegold, Y.; Min, Z. Projecting Global Forest Area towards 2030. For. Ecol. Manag. 2015, 352, 124–133. [Google Scholar] [CrossRef] [Green Version]
  5. Cudlín, P.; Seják, J.; Pokorný, J.; Albrechtová, J.; Bastian, O.; Marek, M. Forest Ecosystem Services Under Climate Change and Air Pollution. In Developments in Environmental Science; Matyssek, R., Clarke, N., Cudlin, P., Mikkelsen, T.N., Tuovinen, J.-P., Wieser, G., Paoletti, E., Eds.; Climate Change, Air Pollution and Global Challenges; Elsevier: Amsterdam, The Netherlands, 2013; Volume 13, pp. 521–546. [Google Scholar]
  6. Benabou, A.; Moukrim, S.; Laaribya, S.; Aafi, A.; Chkhichekh, A.; ElMaadidi, T.; El Aboudi, A. Mapping Ecosystem Services of Forest Stands: Case Study of Maamora, Morocco. Geogr. Environ. Sustain. 2022, 15, 141–149. [Google Scholar] [CrossRef]
  7. Reid, W.V.; Mooney, H.A.; Cropper, A.; Capistrano, D.; Carpenter, S.R.; Chopra, K.; Dasgupta, P.; Dietz, T.; Duraiappah, A.K.; Hassan, R.; et al. Ecosystems and Human Well-Being-Synthesis: A Report of the Millennium Ecosystem Assessment; Island Press: Washington, DC, USA, 2005; 137p. [Google Scholar]
  8. Shukla, P.R.; Skea, J.; Calvo Buendia, E.; Masson-Delmotte, V.; Pörtner, H.-O.; Roberts, D.C.; Zhai, P.; Slade, R.; Connors, S.; Van Diemen, R.; et al. Climate Change and Land: An IPCC Special Report on Climate Change, Desertification, Land Degradation, Sustainable Land Management, Food Security, and Greenhouse Gas Fluxes in Terrestrial Ecosystems; Intergovernmental Panel on Climate Change (IPCC): Geneva, Switzerland, 2019. [Google Scholar]
  9. UNFCCC. Adoption of the Paris Agreement; Report No. FCCC/CP/2015/L.9/Rev.1; UNFCCC: Bonn, Germany, 2015. [Google Scholar]
  10. Harris, N.; Cook-Patton, S.; Gibbs, D.; Lister, K. Young Forests Capture Carbon Quicker than Previously Thought; WRI: Washington, DC, USA, 2020. [Google Scholar]
  11. Piao, S.; Fang, J.; Ciais, P.; Peylin, P.; Huang, Y.; Sitch, S.; Wang, T. The Carbon Balance of Terrestrial Ecosystems in China. Nature 2009, 458, 1009–1013. [Google Scholar] [CrossRef]
  12. Hartig, F.; Dyke, J.; Hickler, T.; Higgins, S.I.; O’Hara, R.B.; Scheiter, S.; Huth, A. Connecting Dynamic Vegetation Models to Data—An Inverse Perspective. J. Biogeogr. 2012, 39, 2240–2252. [Google Scholar] [CrossRef]
  13. Landsberg, J.J.; Sands, P. Physiological Ecology of Forest Production: Principles, Processes and Models, 1st ed.; Academic Press: Amsterdam, The Netherlands; Boston, MA, USA, 2016; ISBN 978-0-12-810206-0. [Google Scholar]
  14. Maréchaux, I.; Langerwish, F.; Huth, A.; Bugmann, H.; Morin, X.; Reyer, C.P.O.; Seidl, R.; Collalti, A.; de Paula, M.D.; Fischer, R.; et al. Modelling forests to address key ecological questions: Lessons learned from different modelling communities and possible future paths. Ecol. Evol. 2021, 11, 3746–3770. [Google Scholar] [CrossRef]
  15. Rodríguez-Suárez, J.A.; Soto, B.; Iglesias, M.L.; Diaz-Fierros, F. Application of the 3PG Forest Growth Model to a Eucalyptus Globulus Plantation in Northwest Spain. Eur. J. For. Res. 2010, 129, 573–583. [Google Scholar] [CrossRef]
  16. Dalmonech, D.; Marano, G.; Amthor, J.S.; Cescatti, A.; Lindner, M.; Trotta, C.; Collalti, A. Feasibility of Enhancing Carbon Sequestration and Stock Capacity in Temperate and Boreal European Forests via Changes to Management Regimes. Agric. For. Meteorol. 2022, 327, 109203. [Google Scholar] [CrossRef]
  17. Testolin, R.; Dalmonech, D.; Marano, G.; Bagnara, M.; D’Andrea, E.; Matteucci, G.; Noce, S.; Collalti, A. Simulating Diverse Forest Management Options in a Changing Climate on a Pinus Nigra Subsp. Laricio Plantation in Southern Italy. Sci. Total Environ. 2023, 857, 159361. [Google Scholar] [CrossRef] [PubMed]
  18. Mahnken, M.; Cailleret, M.; Collalti, A.; Trotta, C.; Biondo, C.; D’Andrea, E.; Dalmonech, D.; Marano, G.; Mäkelä, A.; Minunno, F.; et al. Accuracy, realism and general applicability of European forest models. Glob. Chang. Biol. 2022, 28, 6921–6943. [Google Scholar] [CrossRef] [PubMed]
  19. Landsberg, J.J.; Waring, R.H. A Generalised Model of Forest Productivity Using Simplified Concepts of Radiation-Use Efficiency, Carbon Balance and Partitioning. For. Ecol. Manag. 1997, 95, 209–228. [Google Scholar] [CrossRef]
  20. Almeida, A.C.; Landsberg, J.J.; Sands, P.J. Parameterisation of 3-PG Model for Fast-Growing Eucalyptus Grandis Plantations. For. Ecol. Manag. 2004, 193, 179–195. [Google Scholar] [CrossRef]
  21. Landsberg, J.J.; Waring, R.H.; Coops, N.C. Performance of the Forest Productivity Model 3-PG Applied to a Wide Range of Forest Types. For. Ecol. Manag. 2003, 172, 199–214. [Google Scholar] [CrossRef]
  22. Zhao, M.; Xiang, W.; Peng, C.; Tian, D. Simulating Age-Related Changes in Carbon Storage and Allocation in a Chinese Fir Plantation Growing in Southern China Using the 3-PG Model. For. Ecol. Manag. 2009, 257, 1520–1531. [Google Scholar] [CrossRef]
  23. Stape, J.L.; Ryan, M.G.; Binkley, D. Testing the Utility of the 3-PG Model for Growth of Eucalyptusgrandis×urophylla with Natural and Manipulated Supplies of Water and Nutrients. For. Ecol. Manag. 2004, 193, 219–234. [Google Scholar] [CrossRef]
  24. Dye, P.J.; Jacobs, S.; Drew, D. Verification of 3-PG Growth and Water-Use Predictions in Twelve Eucalyptus Plantation Stands in Zululand, South Africa. For. Ecol. Manag. 2004, 193, 197–218. [Google Scholar] [CrossRef]
  25. Esprey, L.; Sands, P.; Smith, C. Understanding 3PG Using a Sensitivity Analysis. For. Ecol. Manag. 2004, 193, 235–250. [Google Scholar] [CrossRef]
  26. Coops, N.C.; Waring, R.H.; Law, B.E. Assessing the Past and Future Distribution and Productivity of Ponderosa Pine in the Pacific Northwest Using a Process Model, 3-PG. Ecol. Model. 2005, 183, 107–124. [Google Scholar] [CrossRef]
  27. Wei, L.; Marshall, J.D.; Zhang, J.; Zhou, H.; Powers, R.F. 3-PG Simulations of Young Ponderosa Pine Plantations under Varied Management Intensity: Why Do They Grow so Differently? For. Ecol. Manag. 2014, 313, 69–82. [Google Scholar] [CrossRef]
  28. Dye, P. Modelling Growth and Water Use in Four Pinus Patula Stands with the 3-PG Model. S. Afr. For. J. 2001, 191, 53. [Google Scholar]
  29. Landsberg, J.J.; Mäkelä, A.; Sievänen, R.; Kukkola, M. Analysis of Biomass Accumulation and Stem Size Distributions over Long Periods in Managed Stands of Pinus Sylvestris in Finland Using the 3-PG Model. Tree Physiol. 2005, 25, 781–792. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  30. Xenakis, G.; Ray, D.; Mencuccini, M. Sensitivity and Uncertainty Analysis from a Coupled 3-PG and Soil Organic Matter Decomposition Model. Ecol. Model. 2008, 219, 1–16. [Google Scholar] [CrossRef]
  31. Minunno, F.; Xenakis, G.; Perks, M.P.; Mencuccini, M. Calibration and Validation of a Simplified Process-Based Model for the Prediction of the Carbon Balance of Scottish Sitka Spruce (Picea Sitchensis) Plantations. Can. J. For. Res. 2010, 40, 2411–2426. [Google Scholar] [CrossRef]
  32. Nolè, A.; Collalti, A.; Magnani, F.; Duce, P.; Ferrara, A.; Mancino, G.; Marras, S.; Sirca, C.; Spano, D.; Borghetti, M. Assessing Temporal Variation of Primary and Ecosystem Production in Two Mediterranean Forests Using a Modified 3-PG Model. Ann. For. Sci. 2013, 70, 729–741. [Google Scholar] [CrossRef] [Green Version]
  33. Naggar, M. La Régénération Du Cèdre Dans Le Moyen Atlas Central Au Maroc. For. Méditerranéenne 2013, 34, 25–34. [Google Scholar]
  34. Fennane, M.; Ibn Tattou, M. Statistiques et Commentaires Sur l’inventaire Actuel de La Flore Vasculaire Du Maroc. Bull. de l’Inst. Sci. 2012, 34, 1–9. [Google Scholar]
  35. Benabid, A. Biogéographie, Phytosociologie et Phytodynamique Des Cédraies de l’Atlas Cedrus atlantica (Manetti). Ann. De La Rech. For. Au Maroc 1994, 27, 33–60. [Google Scholar]
  36. HCEFLCD. Etudes D’aménagement Concerté des Forêts et Parcours Collectifs des Forêts de la Province D’Ifrane; d’Azrou, F., Ed.; Report; HCEFLCD: Rabat, Morocco, 2007; p. 7. [Google Scholar]
  37. M’hirit, O.; Benzyane, M. Le Cèdre De L’Atlas: Mémoire Du Temps; Editions Mardaga: Wavre, Belgium, 2006; ISBN 978-9981-896-72-7. [Google Scholar]
  38. Moukrim, S.; Lahssini, S.; Naggar, M.; Lahlaoi, H.; Rifai, N.; Arahou, M.; Rhazi, L. Local community involvement in forest rangeland management: Case study of compensation on forest area closed to grazing in Morocco. Rangel. J. 2019, 41, 43–53. [Google Scholar] [CrossRef]
  39. Derak, M.; M’hirit, O.; Mouflih, B.; Et-tobi, M. Influence de la densité et du type de peuplement sur le dépérissement du Cèdre à Sidi M’Guild (Moyen Atlas Marocain). Forêt Méditerranéenne 2008, 1, 23–32. [Google Scholar]
  40. Moukrim, S.; Lahssini, S.; Rifai, N.; Menzou, K.; Mharzi-Alaoui, H.; Labbaci, A.; Rhazi, M.; Wahby, I.W.; El Madihi, M.; Rhazi, L. Modélisation de la distribution potentielle de Cedrus atlantica Manetti au Maroc et impacts du changement climatique. Bois Et Forêts Des Trop. 2020, 344, 3–16. [Google Scholar] [CrossRef]
  41. Joint Research Centre, Institute for Environment and Sustainability; Hiederer, R. Mapping Soil Properties for Europe: Spatial Representation of Soil Database Attributes; Publications Office: Luxembourg, Luxembourg, 2013; Available online: https://data.europa.eu/doi/10.2788/94128 (accessed on 25 May 2021).
  42. Cochran, W.G. Sampling Techniques, 3rd ed.; Wiley series in probability and mathematical statistics; Wiley: New York, NY, USA, 1977; ISBN 978-0-471-16240-7. [Google Scholar]
  43. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sens. Environ. 2017, 202, 18–27. [Google Scholar] [CrossRef]
  44. El Mderssa, M.; Faculté polytechnique de Beni Mellal, Beni Mellal, Morocco. Personal communication, 2019.
  45. IPCC. Guidelines for National Greenhouse Gas Inventories Volume4 Agriculture, Forestry and Other Land Use. Available online: https://www.ipcc-nggip.iges.or.jp/public/2006gl/vol4.html (accessed on 10 January 2021).
  46. Lin, J.C.; Pejam, M.R.; Chan, E.; Wofsy, S.C.; Gottlieb, E.W.; Margolis, H.A.; McCaughey, J.H. Attributing uncertainties in simulated biospheric carbon fluxes to different error sources. Glob. Biogeochem. Cycles 2011, 25, 2. [Google Scholar] [CrossRef]
  47. Luo, Y.; Weng, E.; Wu, X.; Gao, C.; Zhou, X.; Zhang, L. Parameter identifiability, constraint, and equifinality in data assimilation with ecosystem models. Ecol. Appl. Publ. Ecol. Soc. Am. 2009, 19, 571–574. [Google Scholar] [CrossRef] [PubMed]
  48. Li, T.; Sun, X.; Lu, Z.; Wu, Y. A Novel Multiobjective Optimization Method Based on Sensitivity Analysis. Math. Probl. Eng. 2016, 2016, 6012805. [Google Scholar] [CrossRef] [Green Version]
  49. Menberg, K.; Heo, Y.; Augenbroe, G.; Choudhary, R. New Extension Of Morris Method For Sensitivity Analysis of Building Energy Models. 2016. Available online: https://www.researchgate.net/publication/308119619_New_extension_of_Morris_method_for_sensitivity_analysis_of_building_energy_models (accessed on 18 May 2021).
  50. Morris, M.D. Factorial sampling plans for preliminary computational experiments. Technometrics 1991, 33, 161–174. [Google Scholar] [CrossRef]
  51. Franczyk, A. Using the Morris sensitivity analysis method to assess the importance of input variables on time-reversal imaging of seismic sources. Acta Geophys. 2019, 67, 1525–1533. [Google Scholar] [CrossRef] [Green Version]
  52. Collalti, A.; Thornton, P.E.; Cescatti, A.; Rita, A.; Borghetti, M.; Nolè, A.; Trotta, C.; Ciais, P.; Matteucci, G. The Sensitivity of the Forest Carbon Budget Shifts across Processes along with Stand Development and Climate Change. Ecol. Appl. 2019, 29, e01837. [Google Scholar] [CrossRef] [Green Version]
  53. Campolongo, F.; Saltelli, A.; Cariboni, J. From Screening to Quantitative Sensitivity Analysis. A Unified Approach. Comput. Phys. Commun. 2011, 182, 978–988. [Google Scholar] [CrossRef]
  54. Ren, X.; He, H.; Moore, D.J.P.; Zhang, L.; Liu, M.; Li, F.; Yu, G.; Wang, H. Uncertainty analysis of modeled carbon and water fluxes in a subtropical coniferous plantation. J. Geophys. Res. Biogeosci. 2013, 118, 1674–1688. [Google Scholar] [CrossRef]
  55. Ricciuto, D.M.; King, A.W.; Dragoni, D.; Post, W.M. Parameter and prediction uncertainty in an optimized terrestrial carbon cycle model: Effects of constraining variables and data record length. J. Geophys. Res. Biogeosci. 2011, 116, G1. [Google Scholar] [CrossRef] [Green Version]
  56. Zobitz, J.M.; Desai, A.R.; Moore, D.J.P.; Chadwick, M.A. A primer for data assimilation with ecological models using Markov Chain Monte Carlo (MCMC). Oecologia 2011, 167, 599–611. [Google Scholar] [CrossRef] [PubMed]
  57. Liu, C.; Zheng, X.; Ren, Y. Parameter Optimization of the 3PG Model Based on Sensitivity Analysis and a Bayesian Method. Forests 2020, 11, 1369. [Google Scholar] [CrossRef]
  58. Braak, C.J.T. A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: Easy Bayesian computing for real parameter spaces. Stat. Comput. 2006, 16, 239–249. [Google Scholar] [CrossRef]
  59. van Ravenzwaaij, D.; Cassey, P.; Brown, S.D. A simple introduction to Markov Chain Monte-Carlo sampling. Psychon. Bull. Rev. 2018, 25, 143–154. [Google Scholar] [CrossRef] [Green Version]
  60. Trotsiuk, V.; Hartig, F.; Forrester, D.I. R3PG—An r Package for Simulating Forest Growth Using the 3-PG Process-Based Model. Methods Ecol. Evol. 2020, 11, 1470–1475. [Google Scholar] [CrossRef]
  61. Hartig, F.; Minunno, F.; Paul, S. BayesianTools: General-Purpose MCMC and SMC Samplers and Tools for Bayesian Statistics. R Package Version 0.1.7. 2019. Available online: https://CRAN.R-project.org/package=BayesianTools (accessed on 18 May 2021).
  62. ter Braak, C.J.F.; Vrugt, J.A. Differential Evolution Markov Chain with Snooker Updater and Fewer Chains. Stat. Comput. 2008, 18, 435–446. [Google Scholar] [CrossRef] [Green Version]
  63. Gelman, A.; Rubin, D.B. Inference from Iterative Simulation Using Multiple Sequences. Stat. Sci. 1992, 7, 457–472. [Google Scholar] [CrossRef]
  64. McElreath, R. Statistical Rethinking: A Bayesian Course with Examples in R and Stan; Chapman and Hall/CRC: Boca Raton, FL, USA, 2015; ISBN 978-1-4822-5344-3. [Google Scholar]
  65. Raison, R.J.; Myers, B.J. The Biology of Forest Growth Experiment: Linking Water and Nitrogen Availability to the Growth of Pinus Radiata. For. Ecol. Manag. 1992, 52, 279–308. [Google Scholar] [CrossRef]
  66. Lu, Y.; Coops, N.C.; Wang, T.; Wang, G. A Process-Based Approach to Estimate Chinese Fir (Cunninghamia Lanceolata) Distribution and Productivity in Southern China under Climate Change. Forests 2015, 6, 360–379. [Google Scholar] [CrossRef] [Green Version]
  67. Patenaude, G.; Milne, R.; Van Oijen, M.; Rowland, C.S.; Hill, R.A. Integrating Remote Sensing Datasets into Ecological Modelling: A Bayesian Approach. Int. J. Remote Sens. 2008, 29, 1295–1315. [Google Scholar] [CrossRef] [Green Version]
  68. Navarro-Cerrillo, R.M.; Beira, J.; Suarez, J.; Xenakis, G.; Sánchez-Salguero, R.; Hernández-Clemente, R. Growth Decline Assessment in Pinus Sylvestris L. and Pinus Nigra Arnold. Forest by Using 3-PG Model. For. Syst. 2016, 25, 3. [Google Scholar] [CrossRef] [Green Version]
  69. Pinjuv, G.L. Hybrid Forest Modelling of Pinus Radiata D. Don in Canterbury, New Zealand; University of Canterbury: Christchurch, New Zealand, 2006. [Google Scholar]
  70. Forrester, D.I.; Ammer, C.; Annighöfer, P.J.; Avdagic, A.; Barbeito, I.; Bielak, K.; Brazaitis, G.; Coll, L.; del Río, M.; Drössler, L.; et al. Predicting the Spatial and Temporal Dynamics of Species Interactions in Fagus Sylvatica and Pinus Sylvestris Forests across Europe. For. Ecol. Manag. 2017, 405, 112–133. [Google Scholar] [CrossRef] [Green Version]
  71. Forrester, D.I.; Tang, X. Analysing the Spatial and Temporal Dynamics of Species Interactions in Mixed-Species Forests and the Effects of Stand Density Using the 3-PG Model. Ecol. Model. 2016, 319, 233–254. [Google Scholar] [CrossRef]
  72. Collalti, A.; Prentice, I.C. Is NPP Proportional to GPP? Waring’s Hypothesis 20 Years on. Tree Physiol. 2019, 39, 1473–1483. [Google Scholar] [CrossRef]
  73. Waring, R.H.; Landsberg, J.J.; Williams, M. Net Primary Production of Forests: A Constant Fraction of Gross Primary Production? Tree Physiol. 1998, 18, 129–134. [Google Scholar] [CrossRef]
  74. Massoud, E.C.; Xu, C.; Fisher, R.A.; Knox, R.G.; Walker, A.P.; Serbin, S.P.; Christoffersen, B.O.; Holm, J.A.; Kueppers, L.M.; Ricciuto, D.M.; et al. Identification of Key Parameters Controlling Demographically Structured Vegetation Dynamics in a Land Surface Model: CLM4.5(FATES). Geosci. Model Dev. 2019, 12, 4133–4164. [Google Scholar] [CrossRef] [Green Version]
  75. Collalti, A.; Tjoelker, M.G.; Hoch, G.; Mäkelä, A.; Guidolotti, G.; Heskel, M.; Petit, G.; Ryan, M.G.; Battipaglia, G.; Matteucci, G.; et al. Plant Respiration: Controlled by Photosynthesis or Biomass? Glob. Chang. Biol. 2020, 26, 1739–1753. [Google Scholar] [CrossRef]
  76. Collalti, A.; Ibrom, A.; Stockmarr, A.; Cescatti, A.; Alkama, R.; Fernández-Martínez, M.; Matteucci, G.; Sitch, S.; Friedlingstein, P.; Ciais, P.; et al. Forest Production Efficiency Increases with Growth Temperature. Nat. Commun. 2020, 11, 5322. [Google Scholar] [CrossRef]
  77. Farquhar, G.D.; Sharkey, T.D. Stomatal Conductance and Photosynthesis. Annu. Rev. Plant Physiol. 1982, 33, 317–345. [Google Scholar] [CrossRef]
  78. Paul, K.I.; Polglase, P.J.; Richards, G.P. Sensitivity Analysis of Predicted Change in Soil Carbon Following Afforestation. Ecol. Model. 2003, 164, 137–152. [Google Scholar] [CrossRef]
  79. Zaehle, S.; Sitch, S.; Smith, B.; Hatterman, F. Effects of Parameter Uncertainties on the Modeling of Terrestrial Biosphere Dynamics. Glob. Biogeochem. Cycles 2005, 19. [Google Scholar] [CrossRef]
  80. Pappas, C.; Fatichi, S.; Leuzinger, S.; Wolf, A.; Burlando, P. Sensitivity Analysis of a Process-Based Ecosystem Model: Pinpointing Parameterization and Structural Issues. J. Geophys. Res. Biogeosci. 2013, 118, 505–528. [Google Scholar] [CrossRef]
  81. Tatarinov, F.A.; Cienciala, E. Application of BIOME-BGC Model to Managed Forests: 1. Sensitivity Analysis. For. Ecol. Manag. 2006, 237, 267–279. [Google Scholar] [CrossRef]
  82. Bontemps, J.-D.; Bouriaud, O. Predictive Approaches to Forest Site Productivity: Recent Trends, Challenges and Future Perspectives. For. Int. J. For. Res. 2014, 87, 109–128. [Google Scholar] [CrossRef]
  83. Ulrich, D.E.M.; Still, C.; Brooks, J.R.; Kim, Y.; Meinzer, F.C. Investigating Old-Growth Ponderosa Pine Physiology Using Tree-Rings, d13C, d18O, and a Process-Based Model. Ecology 2019, 100, e02656. [Google Scholar] [CrossRef] [Green Version]
  84. Seidl, R.; Thom, D.; Kautz, M.; Martín-Benito, D.; Peltoniemi, M.; Vacchiano, G.; Wild, J.; Ascoli, D.; Petr, M.; Honkaniemi, J.; et al. Forest Disturbances under Climate Change. Nat. Clim. Chang. 2017, 7, 395–402. [Google Scholar] [CrossRef] [Green Version]
  85. Goulden, M.L.; Munger, J.W.; Fan, S.-M.; Daube, B.C.; Wofsy, S.C. Measurements of Carbon Sequestration by Long-Term Eddy Covariance: Methods and a Critical Evaluation of Accuracy. Glob. Change Biol. 1996, 2, 169–182. [Google Scholar] [CrossRef] [Green Version]
  86. López-Serrano, F.R.; Martínez-García, E.; Dadi, T.; Rubio, E.; García-Morote, F.A.; Lucas-Borja, M.E.; Andrés-Abellán, M. Biomass Growth Simulations in a Natural Mixed Forest Stand under Different Thinning Intensities by 3-PG Process-Based Model. Eur. J. For. Res. 2015, 134, 167–185. [Google Scholar] [CrossRef]
  87. Zaher, H.; Sabir, M.; Benjelloun, H.; Paul-Igor, H. Effect of forest land use change on carbohydrates, physical soil quality and carbon stocks in Moroccan cedar area. J. Environ. Manag. 2020, 254, 109544. [Google Scholar] [CrossRef] [PubMed]
  88. M’Hirit, O. Etude Ecologique et Forestiere des Cedraies du Rif Marocain; Essai sur une Approche Multidimensionnelle de la Phytoecologie et de la Productivite du Cedre (Cedrus atlantica Manetti). Ph.D. Thesis, Université de Droit, d’Economie et des Sciences d’Aix-Marseille, Aix-en-Provence, France, 1982. [Google Scholar]
  89. Toth, J. Première approche de la production potentielle du Cèdre de l’Atlas dans le sud de la France. Rev. For. Fr. 1973, 5, 381. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Map of the administrative situation of the Azrou forest and the location of pure cedar stands in this forest.
Figure 1. Map of the administrative situation of the Azrou forest and the location of pure cedar stands in this forest.
Forests 14 00401 g001
Figure 2. Flowchart of the global methodology.
Figure 2. Flowchart of the global methodology.
Forests 14 00401 g002
Figure 3. Flowchart of the DE-MC method (Differential Evolution-Markov chains).
Figure 3. Flowchart of the DE-MC method (Differential Evolution-Markov chains).
Forests 14 00401 g003
Figure 4. Results of the Morris sensitivity analysis. The 55 parameters are listed on the x-axis. A high value of μ * indicates a large influence of the parameter on the output variable while a high value of σ indicates a nonlinear or interaction effect.
Figure 4. Results of the Morris sensitivity analysis. The 55 parameters are listed on the x-axis. A high value of μ * indicates a large influence of the parameter on the output variable while a high value of σ indicates a nonlinear or interaction effect.
Forests 14 00401 g004
Figure 5. The posterior distributions of the reduced parameters of the 3-PG model developed in the present work.
Figure 5. The posterior distributions of the reduced parameters of the 3-PG model developed in the present work.
Forests 14 00401 g005
Table 1. Characteristics of the selected sampling strata.
Table 1. Characteristics of the selected sampling strata.
StratumArea (ha)The Standard Deviation of the Annual Increase in Carbomass ( tC . ha 1 . yr 1 )Number of Selected Plots
Cool exposure/Low density3680.556
Warm exposure/Low density6750.245
Cool exposure/High density1451.055
Warm exposure/High density2581.188
Total1446 24
Table 2. Summary of data collected ex-situ.
Table 2. Summary of data collected ex-situ.
CategoryVariableCollection NameBand NameSpatial/Temporal Resolution
Climate (4)tmp_min“ECMWF/ERA5/DAILY”minimum_2m_air_temperature0.25°/1 day
tmp_max“ECMWF/ERA5/DAILY”maximum_2m_air_temperature0.25°/1 day
tmp_mean“MODIS_006_MOD11A2”LST_Day_1km1000 m/8 days
prcp“ECMWF/ERA5/MONTHLY”total_precipitation0.25°/1 month
srad“ECMWF/ERA5_LAND/MONTHLY”surface_solar_radiation_downwards0.1°/1 month
frost_days“MODIS/006/MOD11A1”LST_Day_1km1000 m/1 day
Site (2)ASWEuropean Soil Database [41]SMU_EU_S_TAWC1000 m/-
NPP (6)Npp“MODIS/006/MOD17A2H”PsnNet500 m/8 days
Table 3. Carbomass models used in this study.
Table 3. Carbomass models used in this study.
ComponentModel
Tree (Aboveground part)
S C O T ( C , H ) = 53.05 C 2.09897 H 0.4063
Stem
S C O T r ( C , H ) = 53.624 C 2.19062 H 0.36418
Foliage
S C O F ( C , H ) = 0.671045 + 0.024967 C H
C: Circumference at breast height, H: Total height, SCOT: Aboveground carbomass, SCOTr: Stem carbomass, SCOF: Foliage carbomass.
Table 4. Conversion factor for each compartment of the tree.
Table 4. Conversion factor for each compartment of the tree.
CompartmentStemFoliageBranchesMean
%Carbon57.4157.3054.3056.43
Table 5. Forest inventory results by sampling stratum.
Table 5. Forest inventory results by sampling stratum.
StratumN ( trees . ha 1 )D (cm)Age (yr)SB ± SD ( tDM . ha 1 )FB ± SD ( tDM . ha 1 )RB ± SD ( tDM . ha 1 )CAI ( cm . yr 1 )
Cool exposure/High density4664917489.25 ± 19.30.50 ± 0.1248.88 ± 15.31.1
Cool exposure/Low density1803510749.78 ± 10.60.19 ± 0.0517.62 ± 6.61.13
Warm exposure/High density79246159268.08 ± 520.92 ± 0.2093.95 ± 13.40.93
Warm exposure/Low density1407114966.68 ± 16.20.16 ± 0.1322.68 ± 9.31.96
N: Density, D: Mean diameter, Age: Mean age, SB: Mean stem biomass, FB: Mean foliage biomass, RB: Mean root biomass, CAI: Mean increment in circumference.
Table 6. Prediction error of the three candidate parametrization datasets for selection.
Table 6. Prediction error of the three candidate parametrization datasets for selection.
Default ModelMedian ModelCalibrated Model
MSE0.320680.262130.26210
Table 7. Initial values, ranges, and posterior distributions of the 20 optimized 3-PG model parameters.
Table 7. Initial values, ranges, and posterior distributions of the 20 optimized 3-PG model parameters.
ParameterUnitInitial ValueRangeModeMean ± Standard DeviationDescription
pFS20-0.6[0.05, 0.8]0.290.39 ± 0.21Foliage stem partitioning at D = 20 cm
aWS-0.05[0, 0.4]0.1170.201 ± 0.109Constant in stem mass v. diameter relationship
pRn-0.2[0.0001, 0.5]0.4660.267 ± 0.141Minimum fraction of NPP to roots
gammaF11/month0.049[0.0001, 0.04]0.01260.0198 ± 0.0117Maximum litterfall rate
TminDegree °C0[ 1 , 8] 0.82 3.25 ± 2.45Minimum temperature for growth
ToptDegree °C19.5[10, 30]23.8723.98 ± 4.06Optimum temperature for growth
TmaxDegree °C35[30, 40]38.5736.11 ± 2.59Maximum temperature for growth
fN0-0.6[0.0001, 1]0.13000.3943 ± 0.2648Value of fN when FR = 0
fNn-1[0, 1]0.870.47 ± 0.28Power of (1-FR) in fN
MaxAgeYears500[350, 550]461452 ± 57Maximum stand age used in age modifier
nAge-4[1, 4.325]2.4772.537 ± 0.970Power of relative age in function for fAge
rAge-0.95[0.0001, 1.4]0.24360.6980 ± 0.3510Relative age to give fAge = 0.5
SLA1m2/kg5.5[5, 30]22.3316.62 ± 7.15Specific leaf area for mature leaves
K-0.2921[0.4, 0.6]0.400.49 ± 0.06Extinction coefficient for absorption of PAR by the canopy
MaxIntrcptn-0.25[0.1, 0.4]0.360.25 ± 0.087Maximum proportion of rainfall evaporated from canopy
alphaCxmolC/molPAR0.04212129[0.02, 0.09]0.04930.0468 ± 0.0188Canopy quantum efficiency
Y-0.47[0.44, 0.51]0.480.47 ± 0.02Ratio NPP/GPP
MaxCondm/s0.02[0.001, 0.03]0.0220.016 ± 0.008Maximum canopy conductance
CoeffCond1/mBar0.05[0.0001, 0.07]0.00300.0367 ± 0.0196Defines stomatal response to VPD
BLcondm/s0.2[0.0001, 0.3]0.11050.1479 ± 0.0861Canopy boundary layer conductance
Sources for optimized parameters ranges: [22,29,30,65,66,67,68,69,70,71].
Table 8. Carbon increment per stratum unit.
Table 8. Carbon increment per stratum unit.
StatumArea (ha)Unit Carbon Increase ± SD ( tC . ha 1 . yr 1 )Stratum Carbon Increase ( tC . yr 1 )
Cool exposure/Low density3685.20 ± 0.341913
Warm exposure/Low density6755.08 ± 0.463429
Cool exposure/High density1457.24 ± 0.771049
Warm exposure/High density2585.92 ± 0.411527
Pure cedar forest14465.477918
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Boukhris, I.; Lahssini, S.; Collalti, A.; Moukrim, S.; Santini, M.; Chiti, T.; Valentini, R. Calibrating a Process-Based Model to Enhance Robustness in Carbon Sequestration Simulations: The Case of Cedrus atlantica (Endl.) Manetti ex Carrière. Forests 2023, 14, 401. https://doi.org/10.3390/f14020401

AMA Style

Boukhris I, Lahssini S, Collalti A, Moukrim S, Santini M, Chiti T, Valentini R. Calibrating a Process-Based Model to Enhance Robustness in Carbon Sequestration Simulations: The Case of Cedrus atlantica (Endl.) Manetti ex Carrière. Forests. 2023; 14(2):401. https://doi.org/10.3390/f14020401

Chicago/Turabian Style

Boukhris, Issam, Said Lahssini, Alessio Collalti, Said Moukrim, Monia Santini, Tommaso Chiti, and Riccardo Valentini. 2023. "Calibrating a Process-Based Model to Enhance Robustness in Carbon Sequestration Simulations: The Case of Cedrus atlantica (Endl.) Manetti ex Carrière" Forests 14, no. 2: 401. https://doi.org/10.3390/f14020401

APA Style

Boukhris, I., Lahssini, S., Collalti, A., Moukrim, S., Santini, M., Chiti, T., & Valentini, R. (2023). Calibrating a Process-Based Model to Enhance Robustness in Carbon Sequestration Simulations: The Case of Cedrus atlantica (Endl.) Manetti ex Carrière. Forests, 14(2), 401. https://doi.org/10.3390/f14020401

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop