Identification of key parameters controlling demographically structured vegetation dynamics in a land surface model: CLM4.5(FATES)

Vegetation plays an important role in regulating global carbon cycles and is a key component of the Earth system models (ESMs) that aim to project Earth's future climate. In the last decade, the vegetation component within ESMs has witnessed great progress from simple “big-leaf” approaches to demographically structured approaches, which have a better representation of plant size, canopy structure, and disturbances. These demographically structured vegetation models typically have a large number of input parameters, and sensitivity analysis is needed to quantify the impact of each parameter on the model outputs for a better understanding of model behavior. In this study, we conducted a comprehensive sensitivity analysis to diagnose the Community Land Model coupled to the Functionally Assembled Terrestrial Simulator, or CLM4.5(FATES). Specifically, we quantified the first- and second-order sensitivities of the model parameters to outputs that represent simulated growth and mortality as well as carbon fluxes and stocks for a tropical site with an extent of 1×1^∘. While the photosynthetic capacity parameter (V_c,max25) is found to be important for simulated carbon stocks and fluxes, we also show the importance of carbon storage and allometry parameters, which determine survival and growth strategies within the model. The parameter sensitivity changes with different sizes of trees and climate conditions. The results of this study highlight the importance of understanding the dynamics of the next generation of demographically enabled vegetation models within ESMs to improve model parameterization and structure for better model fidelity.

1 Introduction

Earth system models (ESMs) are abstract representations of nature used to simulate physical, chemical, and biological processes across the interacting domains of the Earth system to estimate past, present, and future climate (Claussen et al., 2002; Dunne et al., 2012; Arora et al., 2013; Hurrell et al., 2013). Land surface models (LSMs), the land component of ESMs, are capable of representing vegetation dynamics through the use of dynamic global vegetation models (Foley et al., 1996; Cox et al., 2000; Krinner et al., 2005; Friedlingstein et al., 2006; Sato et al., 2007; Arora et al., 2013). The first-generation dynamic vegetation models represent plant communities and their competition using a single area-averaged representation of plant functional types (PFTs) within each land grid cell (Cox et al., 2000; Pan et al., 2002; Hickler et al., 2004). Recently, a number of vegetation models that can represent plant demographic processes have emerged to better capture coexistence and competition driven by light competition between different sizes of trees within a vertical canopy structure at different successional stages (Moorcroft et al., 2001; Thonicke et al., 2001; Sitch et al., 2003; Hickler et al., 2004; Fisher et al., 2010; Scheiter et al., 2013; Fisher et al., 2018). These demographic models allow comparison with many more observed vegetation processes than first-generation models but also contain more degrees of freedom leading to great complexity.

LSMs typically contain a suite of different parameters to resolve the carbon, water, and energy fluxes and pools at the land–atmosphere interface (Noilhan and Planton, 1989; Bastidas et al., 1999; Gupta et al., 1999; Masson et al., 2003; Sargsyan et al., 2014). Many of these parameters can be estimated directly in the field, but others are difficult or impossible to measure due to various complications such as the abstract representation of processes, technological limitations, or spatial/temporal aggregation (Entekhabi and Eagleson, 1989; Kumar et al., 2006). Parameters that are observable in the field are also often subject to large natural variability, including changes through space and time (Wood et al., 1992; Masson et al., 2003; Fisher et al., 2015). For example, vegetation parameters can be used to describe different root profiles (Vrugt et al., 2001; Zeng, 2001; Massoud et al., 2019a) or photosynthetic capacities (Leuning, 2002; Rogers, 2014); however, model parameter values are often taken from literature publications or databases, and may not represent local variation or capture seasonal or ontogenetic changes. For parameters of critical importance, even a small difference can lead to significant divergence for multi-model ensemble projections or uncertainty in model predictions from different models (Sitch et al., 2008; Dietze et al., 2014; McDowell et al., 2016; Rogers et al., 2017). Since parameters are often defined in simulations with limited prior knowledge of their mean values and variation (O'Hagan and Leonard, 1976; Kitanidis, 1986; Geromel, 1999), model uncertainty or sensitivity analyses are typically required to adequately quantify the uncertainty in model outputs and importance of parameters to guide model calibration and improvement.

There are two types of uncertainty and sensitivity analysis studies. One type of study aims to understand the model behaviors by exploring the baseline sensitivity of model outputs to parameter changes, which is normally an equal amount of deviation from the mean values of default parameters. This is commonly referred to as model sensitivity or elasticity analysis (e.g., Benton and Grant, 1999; Pappas et al., 2013; Menberg et al., 2016; Collalti et al., 2019). Another type of study aims to quantify the amount of uncertainty in model outputs and the corresponding contributions to this uncertainty by different sources, which is commonly referred to as uncertainty quantification (e.g., Xu et al., 2010; Dietze et al., 2014). It is possible that a model output is very sensitive to a particular parameter in the sensitivity analysis study, but the parameter could contribute to a small amount of uncertainty in the model output if this parameter contains a low level of variation (Dietze et al., 2014). Both types of studies are useful for model development and applications with sensitivity analysis studies focusing on understanding the baseline of model behaviors and uncertainty quantification studies focusing on guiding field and laboratory measurements. Despite the need for such studies, systematic investigation of the parameter sensitivity and output uncertainty of LSMs is not standard practice, potentially on account of the high dimensionality involved (although see Zaehle et al., 2005; Fisher et al., 2010; Pappas et al., 2013).

Today, many uncertainty and sensitivity analysis techniques are available (Sobol', 1990; Helton, 1993; Saltelli et al., 2000; Razavi and Gupta, 2016). Some of these methods examine the response of the outputs by varying input parameters one at a time and holding other parameters at their default values (Saltelli et al., 2000). However, the sensitivity index derived by this type of assessment depends on the default values of the other parameters, and the assumption that these values are satisfactory is questionable (e.g., Da Rocha et al., 1996; Sen et al., 2001; Groenendijk et al., 2011; Schwalm et al., 2010) since the LSM predictions are strongly tied (through feedbacks of momentum, energy, mass, and biogeochemistry) to the differences in their representation of the land surface (Crossley et al., 2000; Rosolem et al., 2013). Therefore, it is desirable to use more extensive sensitivity analysis techniques that examine the response of model outputs averaged over the variation of all the parameters. These “global” methods are generally preferred when computing power is not a limiting factor, as they require a relatively large number of ensemble runs. A sensitivity analysis is considered to be global when all the input factors are varied simultaneously and the sensitivity is evaluated over their entire range of interest (McRae et al., 2001; Xu and Gertner, 2008; Zhou et al., 2008). For sensitivity analysis studies, the entire range of interest could be a certain percentage deviation of default values of parameters. For uncertainty quantification studies, the entire range of interest could be the distributions of parameters estimated from laboratory measurements, field observations, and expert knowledge. Campolongo et al. (2000) suggested classifying local and global sensitivity analysis based largely on the extent of the input variable range that the technique assesses; however, this classification is ambiguous, as it depends on whether the range is sufficiently large to be perceived as global (Song et al., 2015).

The goal of this study is to conduct a comprehensive sensitivity analysis for a land surface model (Community Land Model) coupled to a demographic vegetation model (Functionally Assembled Terrestrial Simulator), or CLM4.5(FATES), at a tropical site with an extent of 1×1^∘ to (1) understand the baseline model behaviors of vegetation carbon stocks and fluxes and vegetation demography in relation to different model parameters, and (2) provide directions for improved model parameterization toward a better model fit to observations. Specifically, we aim to answer the following question: what are the main parameter controls on vegetation processes such as growth and mortality and on the resulting dynamics of carbon fluxes and stocks? Based on our understanding of simulated processes in CLM4.5(FATES), we propose to test three hypotheses. Our first hypothesis is related to photosynthetic capacity. The carbon input for vegetation growth is through photosynthesis and in most LSMs, it is simulated based on the Farquhar model (Farquhar, 1989) with the photosynthetic capacity represented by the maximum carboxylation rate at 25 ^∘C (V_c,max25) and maximum electron transport rate at 25 ^∘C (J_max25). J_max25 is simulated in proportion to V_c,max25 in many models, and previous sensitivity analysis studies (Pappas et al., 2013; Sargsyan et al., 2014; Dietze et al., 2014) have shown that V_c,max25 is generally an important parameter that affects simulated carbon fluxes. Therefore, we hypothesize that the photosynthetic capacity parameter, V_c,max25, is a key control on simulated carbon fluxes in CLM4.5(FATES) (H1). Second, for demographic models, the allometry of trees determines the amount of carbon input to different tissues (e.g., leaf, root, and stem). If more carbon is allocated to leaf compared to stem, the tree will have a higher productivity but this can also lead to lower stem growth and thus less height growth for light competition. Thus, we hypothesize that allometry parameters are important for vegetation growth and long-term carton stocks (H2), as they will determine plant's growth strategies. Finally, the carbon stock for vegetation is affected not only by the input of carbon through photosynthesis but also by the loss of carbon through mortality. We hypothesize that the parameters determining mortality are important drivers of the long-term vegetation carbon stocks (H3), as they will control carbon turnover time.

2 Materials and methods

2.1 CLM4.5(FATES)

CLM4.5(FATES) is an open-source land surface model coupled with a demographically structured dynamic vegetation model designed to predict climate–vegetation interactions. The land surface Community Land Model (CLM) represents surface heterogeneity and simulates land biogeophysics, the hydrologic cycle, biogeochemistry, human dimensions, and ecosystem dynamics (Oleson et al., 2013). It is used within various Earth system modeling frameworks, including the Community Earth System Model (CESM) and the Norwegian Earth System Model (NorESM) (Lawrence et al., 2011; Bonan et al., 2011). It is also the baseline of the land component within the Exascale Energy Earth System Model (E3SM) (Golaz et al., 2019). The vegetation model FATES is developed from the ecosystem demography (ED) model, which scales up the behavior of forest ecosystems by aggregating individual trees into representative “cohorts” based on their size and PFT, and by aggregating groups of cohorts into representative “patches” (conceptually similar to a forest plot) that explicitly tracks the time between disturbances (Moorcroft et al., 2001). The main property of the ED concept that differs from most commonly used “big-leaf” models is the capacity to predict distribution, structure, and composition of vegetation directly from their given physiological traits described by the model parameterization (Fisher et al., 2015). This is achieved via the means of trait filtering, whereby plant traits affect plant growth and survival, growth in turn affects the acquisition of light resources, and feeds back onto growth, survival, and reproduction. Differences in growth, survival, and reproduction rates thus directly control the relative distributions of vegetation types and their traits as well as the overall carbon stocks. See supplementary model description in Fisher et al. (2015) for details on specific components of the model structure. FATES represents vegetation using size-structured groups of plants (cohorts) which coexist on various successional trajectory-based land units. FATES simulates growth by integrating photosynthesis across different leaf layers for each cohort. The model allocates photosynthetic carbon to different tissues such as leaf, root, and stem based on the allometry of different vegetation types. Mortality is an important driver for the simulated forest dynamics in FATES. FATES includes five modes of mortality: (1) fixed background mortality, (2) hydraulic failure based on a threshold of very low soil moisture; (3) carbon starvation resulting from the depletion of carbon storage in plants (see Appendix C for details); (4) impact mortality resulting from the falling of big trees; and (5) fire (Fisher et al., 2015). CLM and FATES are coupled to exchange carbon and water between vegetation and soil through a common interface. FATES is designed to be modular and currently can be turned on within two land surface models: the CLM and E3SM land model (ELM). Depending on the purpose of different studies, CLM4.5(FATES) can be simulated with different modes including point mode for individual sites, regional mode for watershed or regional scales, and global mode for the global scale.

In this original version of CLM4.5(FATES), there are two challenges for the model to simulate tropical forests. First, it is difficult for the model to represent the coexistence of PFTs due to the dominance of growth and reproductive feedbacks and potentially the absence of additional stabilizing mechanisms (Fisher et al., 2010, 2018); therefore, in this initial analysis, we focus only on a single broadleaf evergreen tree PFT, which is a typical vegetation type for the study region (the Amazon). We want to point out that, because of the high species diversity in the tropics, it is always a challenge for models to capture diverse traits with a limited number of PFTs in typical ESMs. By limiting our sensitivity analysis to one PFT in this study, our sensitivity analysis will help us understand the main control on demographic rates of growth and mortality that will essentially affect the outcome of competition for multiple PFTs. Thus, we expect that our sensitivity analysis can be used to guide the selection of traits for the representation of trait trade-off for diverse tropical forests and to improve the simulation of PFT coexistence for model calibration and improvement. Second, the model generally underestimates leaf area index (LAI). We expect that our sensitivity analysis will be used as a guidance to adjust identified key model parameters in order to better fit model predictions to the observations.

The CLM4.5(FATES) tracks different size classes of plants (generally >10 size classes) through time. To facilitate our analysis, we aggregate cohorts into three size categories: small (<10 cm), medium (10–50 cm), and large trees (>50 cm). For sensitivity analysis of each size category (small, medium, and large trees), we choose to average the outputs over 30-year intervals. This is done with the view that a large amount of variability in model outputs could be caused by the transient and abrupt changes across size classes, and thus only a small amount of variability is affected by parametric variations. Our analysis shows that the identified key parameters and the corresponding magnitude of sensitivity are similar with averaging over different numbers of years from 20 to 40 years (Fig. D1).

2.2 Sensitivity analysis: the FAST method

Global sensitivity analysis aims at quantifying the contributions of input variables to the variability of the outputs of a physical model by simultaneously sampling values of parameters from their corresponding statistical distributions. There are many methods for global sensitivity analysis. Two popular variance-based approaches are the Sobol method (Sobol', 1990) and the Fourier amplitude sensitivity test (FAST) (Cukier et al., 1973). The Sobol method has received much attention since it provides a clear description of the importance index of model parameters based on variance decomposition. However, the full description requires the evaluation of 2ⁿ Monte Carlo integrals (Sudret, 2008), which is not practically feasible unless n is low (n here represents the dimensionality of the model or the number of active parameters). Compared to Sobol's method, FAST is more computationally efficient. It can be used effectively for nonlinear and nonmonotonic models (Sudret, 2008; Xu and Gertner, 2011a). FAST uses a periodic sampling approach to draw samples from the parameter space defined by probability distributions with a characteristic periodic signal for each parameter. These samples will then be fed into the model for ensemble simulations. Finally, a Fourier transformation is applied to decompose the variance of a model output into partial variances contributed by different model parameters based on the characteristic periodic signal assigned for each parameter. Only the first-order sensitivity indices referring to the “main effect” of parameters were calculated in the original method. In the 1990s, an extended FAST method able to calculate sensitivity indices referring to “total effect” was developed (Sobol', 1990; Archer et al., 1997; Saltelli et al., 1999). This “total effect” of a parameter's sensitivity refers to the sum of a parameter's individual contribution (first-order sensitivity) and the contribution from its interaction with other parameters (higher-order sensitivity) on the overall variance of the model output; that is, the total effect includes all the higher-order interactions. Xu and Gertner (2011a) further derived equations within the FAST framework to calculate specific higher-order interactions for different sampling approaches. The FAST method has found widespread use in many different fields of study including sensitivity analysis of the parameters of models that represent the land surface (Collins and Avissar, 1994), chemical reaction (Haaker and Verheijen, 2004), nuclear waste disposal (Lu and Mohanty, 2001), erosion (Wang et al., 2001), hydrologic systems (Francos et al., 2003), atmospheric systems (Kioutsioukis et al., 2004), crop growth (Wang et al., 2013), or matrix population and forest landscape models (Xu and Gertner, 2009; Xu et al., 2009).

In this study, we quantify both the first- and second-order sensitivities of the model parameters using FAST. It is possible to identify higher-order interactions with FAST; however, because of the sample size limitations for a larger trivariate parameter space, the FAST-based estimation of third-order sensitivity indices would be less reliable (Xu and Gertner, 2011a). Specifically, the first-order sensitivity is used to measure the importance of the variations in one parameter to the model outputs. If one parameter x_i is important to a model output y at time t (i.e., y(t)), we expect that the mean value of y(t) will change substantially with different values of x_i. Statistically, we expect to see a large variance of the expected value of y(t) given x_i (i.e., large V(E(y(t)|x_i)), where E(⋅) is the expected value of the output, V(⋅) is the variance calculated in the parameter space). Similarly, if the combined impact of x_i and x_j is important, we expect to see a large variance of the expected value of y(t) given x_i and x_j (i.e., large V(E(y(t)|x_i,x_j))). Therefore, we calculate the first- and second-order sensitivities, ${\mathit{\alpha }}_{x_i}$ and ${\mathit{\alpha }}_{x_i{x}_j}$, respectively, of the model parameters for each output of interest and at each time step as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\alpha }}_{x_i}\left(t\right)={\displaystyle \frac{V\left(E\right(y\left(t\right)\mathrm{|}{x}_i\left)\right)}{V\left(y\right(t\left)\right)}}\text{(2)}& {\displaystyle }\begin{array}{rl}& {\mathit{\alpha }}_{x_i,x_j}\left(t\right)=\\ & {\displaystyle \frac{V\left(E\right(y\left(t\right)\mathrm{|}{x}_i,x_j\left)\right)-V\left(E\right(y\left(t\right)\mathrm{|}{x}_i\left)\right)-V\left(E\right(y\left(t\right)\mathrm{|}{x}_j\left)\right)}{V\left(y\right(t\left)\right)}},\\ & i\ne j,\end{array}\end{array}$$

where V(y(t)) is the total variance of model output y(t). In FAST, the variances are estimated through periodic samples in the θ space between 0 and 2π, which are linked to the samples in the parameter space through a search function. Further details on the FAST toolbox used for this study can be found in Xu and Gertner (2007, 2009, 2011a) or Xu et al. (2009). We are aware that FAST can provide robust estimates of the sensitivity coefficients in high dimensions (Wang et al., 2013), especially since the CPU demands of CLM4.5(FATES) mandates application of a method like FAST due to its ability to derive sensitivity values with sparse sampling. Due to the random parametric sampling, there will be errors in the estimates of sensitivity indices. In this study, we estimate the standard error of FAST-based sensitivity index derived by Xu and Gertner (2011b), with lower errors for a larger sample size.

To better understand how parameters affect specific CLM4.5(FATES) output variables (i.e., the relationship between model parameters and outputs), we also fitted cubic splines to the scatterplots between samples of parameters identified as important by FAST and the corresponding output variable of interest using the R SemiPar package (Ruppert et al., 2003).

Figure 1 Recycled climate drivers for the study area including annual mean precipitation, relative humidity, and air temperature for the years 1948–1972. The annual radiation and air pressure are not plotted as they are quite stable across years.

2.3 Parameter selection

In total, there are more than 200 parameters for all land surface processes including surface energy exchange, hydrology, biogeochemistry, plant physiology, and demographic processes within CLM4.5(FATES). In this study, we focus solely on vegetation components and select 87 parameters that are relevant to vegetation processes, including parameters for photosynthetic processes, temperature response, allometry description, radiative transfer, recruitment, turnover, and mortality. See Tables D1–D4 in the Appendix for a complete list of the parameters used in this study, with corresponding description, units, default values, and applied ranges. Refer to Appendix A for the allometry equations, Appendix B for the temperature response curve (photosynthesis) equations, and Appendix C for the carbon storage equations used in CLM4.5(FATES).

To conduct the sensitivity analysis, we have extracted many parameters in the model that were “hard-wired” in CLM4.5(FATES). The FAST algorithm requires valid ranges to be chosen for each parameter, which creates the possible parameter space to sample from. In theory, each parameter has a corresponding observational distribution that produces the ideal space for sampling (LeBauer et al., 2013). However, in this study, there are both a large parameter set and a scarcity of appropriate data sources for Amazonian forests for many of the relevant quantities; therefore, obtaining a robust data-supported distribution for each parameter was difficult. Because we only aim to understand the baseline model structure, the parameter ranges in this study were generated by applying a uniform distribution over a range that spans ±15 % of the default parameter values of CLM4.5(FATES) (i.e., default parameter values for tropical evergreen trees). We choose a rather conservative range of ±15 % of the default CLM4.5(FATES) values so that global sensitivity indices can be estimated in the reasonable vicinity of the default parameters. We suggest that a more robust uncertainty analysis based on realistic parameter ranges is needed for guidance on additional field measurements.

Figure 2 Simulated temporal dynamics in diameter at breast height (dDBH; cm yr⁻¹ tree⁻¹) and the corresponding first-order parametric sensitivity indices. The left panels show the simulated ranges of dDBH for (a) all, (b) small (diameter <10 cm), (c) medium (10 cm < diameter <50 cm), and (d) large trees (diameter >50 cm). Shown is the mean simulation (black line) with 95 % spread of the simulation ensemble. Right panels show the sensitivity for the top six most important parameters for (e) density of all trees, (f) small tree density, (g) medium tree density, and (h) large tree density, in order of importance based on the mean parametric sensitivity across years (red is the most important and blue is the least important). The jumps seen in 10, 40, 70, and 100 years for small, medium, and large trees are due to the temporal averaging mentioned in the materials and methods section. The figures also show sensitivities of the remaining parameters in light grey (first-order sensitivity index for all other parameters) as well as the sensitivity of parameter interactions in dark grey (higher-order sensitivity index for all parameters).

Using FAST, 5000 parameter combinations are sampled from the parameter space. The sample size was determined using the heuristic method of Xu and Gertner (2011a), where it is appropriate to use 100 times the number of effective (important) parameters. The 5000 model runs cost about 32 CPU hours for each simulation, and thus we ran our simulations for a total of 160 000 CPU hours on the Los Alamos National Laboratory (LANL) Conejo supercomputer.

In this analysis, we assume the majority of CLM4.5(FATES) parameters to be non-correlated with uniform probability because our study is focused on the model parametric sensitivity for model behaviors and there is a limitation of data for estimating covariance among the 80+ parameters. However, we do need to take care of the correlation among parameters in the temperature response functions (Appendix B) in order to generate realistic temperature response curves. These parameters are tested for correlation using a published dataset (Leuning, 2002), which showed that the photosynthetic parameters for activation energy (e.g., $V_{\mathrm{c},\max ,\mathrm{ha}}$) are not necessarily correlated with the other photosynthetic parameters. However, the parameters for deactivation energy (e.g., $V_{\mathrm{c},\max ,\mathrm{hd}}$) and those related to entropy terms (e.g., $V_{\mathrm{c},\max ,\mathrm{se}}$) are highly correlated, as expected (correlation of 0.99+). Thus, each of these parameters' samples are generated from the same location in their relative parameter spaces, which maintains their correlation.

Figure 3 Simulated temporal dynamics in tree density (NPLANT; N ha⁻¹) and the corresponding first-order parametric sensitivity indices. The left panels show the simulated ranges of tree density for all trees (a) and the corresponding fraction of (b) small (diameter <10 cm), (c) medium (10 cm < diameter <50 cm), and (d) large trees (diameter >50 cm). Right panels show the sensitivity for the top six most important parameters for (e) all, (f) small, (g) medium, and (h) large trees, in order of importance. See Fig. 2 for details on legends.

2.4 Data and model setup

In this study, the CLM4.5(FATES) model simulations are set up for a 1^∘ by 1^∘ grid in a moist-tropical forest in the state of Pará, the Amazon, Brazil (7^∘ S, 55^∘ W), which is a default tropical setup for CLM. The climate conditions for this site are from Qian et al. (2006), representative of data from 1948 to 1972 and recycled for the 130-year simulations (Fig. 1). The CO₂ concentration is set as 284.7 ppm. No nitrogen deposition is simulated, as FATES currently does not have the nutrient limitation yet. We initialized the runs with a near-bare ground or a state with no vegetation but available seeds, and simulated the forest dynamics for 130 years, which we determined was enough time for the ecosystem to reach equilibrium because simulated outputs and corresponding sensitivity values for biomass, basal area, and various carbon fluxes had stabilized by this time. By choosing to start from near-bare ground and running the model until it reaches a quasi-steady-state size distribution, rather than by examining short runs initialized from observed initial forest size distributions (e.g., Dietze et al., 2014), we are deliberately allowing the ecosystem demographic structure itself to be an outcome of the parametric variance rather than a separate, possibly non-self-consistent, initial condition variance. The fire component is turned off in view that the study site has limited fire disturbances.

3 Results

In this section, we highlight the outputs of CLM4.5(FATES) from the 5000 simulations obtained for the FAST analysis and then show the important parameters that control variance in the outputs. We first investigate the forest demographic dynamics, diagnosing the growth and mortality processes simulated in CLM4.5(FATES), i.e., outputs representing the change in diameter at breast height (dDBH), the mortality rate, and the resulting basal area (BA). Then, we analyze the carbon fluxes and stocks in the model simulations including gross primary production (GPP), net primary production (NPP), LAI, and total forest biomass.

3.1 Forest demographic dynamics: growth and mortality

One of the key properties of CLM4.5(FATES) is that vegetation is represented as cohorts of varying sizes for more realistic simulation of light competition in the canopy. To illustrate how different parameters impact different size classes of trees, we group various cohorts of trees into three size classes for analysis purposes: small, medium, and large trees. Since the model runs are initialized from a near-bare-ground state, all simulated plants are considered “small” with an initial density of half-centimeter diameter saplings.

For the stem growth (dDBH averaged per tree; Fig. 2a–d), small trees have lower rates of growth compared to medium or large trees as they are mostly in the understory and thus have a lack of light for photosynthesis. However, the fraction of overall stem diameter growth is dominated by the small trees (Fig. D2) due to their high densities (Fig. 3b). The trees grow faster at the beginning of the simulation when the canopy has not reached full closure. Correspondingly, the parametric sensitivities tends to vary at the beginning of the simulation and then become stable through time. The most sensitive parameters for tree growth are the target storage carbon and stem allometry parameters; however, the importance magnitude varies with time and sizes of trees (Fig. 2e–h). The stem allometry is the most sensitive parameter at the beginning of simulation (<20 years), but the target carbon storage parameter becomes dominant after simulation year 70 (Fig. 2a). We observe that the stem allometry coefficient c is the dominant parameter that controls dDBH for medium and large trees, and the target carbon storage is the most important parameter for small trees. A higher value of stem allometry coefficient c, or a higher allocation of carbon to stem, will lead to a faster growth of diameter at breast height (DBH) in the initial life stage of small trees (Fig. D3a). However, for medium and large trees, a higher allocation of carbon to stem can lead to lower proportion of carbon allocated to leaves for productivity and thus a slower DBH growth (Fig. D3b, c). This outcome supports hypothesis H2, which states the importance of allometric parameters. The target carbon storage determines the target amount of carbon for the plant to store relative to the leaf biomass (see Appendix C for details). Smaller trees have less stem biomass and are less impacted by the stem allometry coefficient c parameter. Furthermore, small trees are vulnerable to changes in the amount of target carbon storage which affects carbon allocation to growth (see Eq. C3 in Appendix C). Our sensitivity analysis also shows specifically important parameters for different sizes of trees. For example, leaf allometry is important for small trees, V_c,max25 for medium trees, and seed allocation for large trees.

Figure 4 Simulated temporal dynamics in tree mortality rates (fraction yr⁻¹) and the corresponding first-order parametric sensitivity indices. The left panels show the mortality rate for (a) all, (b) small (diameter <10 cm), (c) medium (10 cm < diameter <50 cm), and (d) large trees (diameter >50 cm). Right panels show the sensitivity for the top six most important parameters for (e) all, (f) small, (g) medium, and (h) large trees, in order of importance. See Fig. 2 for details on legends.

In this analysis, carbon starvation emerged as the main driver for tree mortality (Fig. D4). The carbon-starvation-based mortality uses a threshold of carbon storage to trigger mortality (see Appendix C). Under shaded conditions, lower carbon stores caused by the balance of NPP, respiration, and tissue growth/maintenance should lead to a higher mortality rate. As expected, the smaller tree size classes have much higher mortality rates (Fig. 4b–d). The first-order sensitivity analysis of predicted mortality rate (percentage of mortality per year) shows that the dominant parameter for predicting mortality of large trees is the target carbon storage (Fig. 4h); however, for small and medium trees, other parameters such as allometric and photosynthetic parameters that could potentially determine their height growth and competitive advantages in the canopy are also important (Fig. 4f, g). Specifically, for medium-sized trees, the mortality rate is affected by both the stem allometry coefficient c and targeted carbon storage (Fig. 4g). For the small trees, important parameters include the photosynthetic capacity parameter (V_c,max25), stem allometry coefficient c, mortality rate under stress, and maintenance respiration, with the target carbon storage having high sensitivity for small trees in the early years (Fig. 4f).

The simulated basal area (BA) of the forest, which is the total stem cross-sectional area per ground surface area, results from the combination of both DBH growth and mortality. The BA reaches equilibrium for different sizes of trees around year 70 (Fig. 5a). Our FAST analysis shows that a key parameter that controls BA in different tree size classes is the stem allometry coefficient c (Fig. 5e–h), which is a major parameter that determines the DBH growth (Fig. 2). We also found that the target carbon storage parameter that dominantly controls mortality is an important parameter for the simulated BA (Fig. 5e–h). Different from parameters important for DBH growth at the individual tree level and mortality rate, a new parameter that becomes important for BA of small and medium trees is the minimum crown spread, which determines the ratio of crown radius to DBH. A larger crown spread can lead to a smaller number of trees in the canopy and thus a lower BA (Fig. D5). The identified important parameters for the simulated tree density and fraction of trees are very similar to those identified for the simulated BA, except that the leaf allometry coefficient b becomes very important for simulated small tree densities (Fig. 3e–h) and minimum height for fraction of trees (Fig. D6).

For the second-order sensitivity analysis, parametric interactions between stem allometry coefficient c and the proportion of carbon for seed allocation, and target carbon storage are found to be important for the prediction of total BA (Fig. D7). For trees of different sizes, parametric interactions between stem allometry coefficient c and minimum crown spread, target carbon storage, and maximum DBH are important for small, medium, and large trees, respectively. For the prediction of dDBH and mortality, the contributions of most parametric interactions are relatively small except for large trees (Fig. D8). The interactions between stem allometry coefficient c and the proportion of carbon for seed allocation, maximum DBH, and stem allometry coefficient b are important for the prediction of dDBH for large trees. With respect to large tree mortality, the interaction between stem allometry coefficient c and target carbon storage is found to be important (Fig. D9).

3.2 Forest carbon cycles: carbon fluxes and stocks

To investigate the key parametric control on carbon fluxes and stocks, we specifically investigate parameter sensitivities for GPP, NPP, LAI, and total forest biomass. Our results show that GPP and NPP increased consistently for the first 10 years of the simulations, which is expected for a forest growing from bare ground (Fig. 6). However, within a fairly short period of 5–10 years, GPP, NPP, and LAI and their variance reached a quasi-stable rate. This amount of time to reach equilibrium is much shorter compared to the basal area (Fig. 6a) and the total biomass accumulations (Fig. 6d).

The first-order sensitivity analysis based on FAST shows that, for carbon fluxes of GPP and NPP, the photosynthetic capacity parameter (V_c,max25) is the most sensitive parameter (Fig. 6e, f), which supports hypothesis H1. Furthermore, specifically for NPP, the respiration parameters such as the growth respiration fraction and leaf maintenance respiration rate show high sensitivity (Fig. 6f). For LAI, the leaf allometry coefficient b is the most important, as it determines carbon allocation for leaves (Fig. 6g). The stem allometry coefficient c is the most important for total biomass (Fig. 6h), as it determines carbon allocation to the stem, which supports hypothesis H2. See Fig. D10 for an easier comparison of parametric sensitivities for different model outputs. A common sensitive parameter is the target carbon storage, which is important for GPP, NPP, LAI, and total biomass. This results from the fact that the target carbon storage is a key driver for mortality especially for medium and large trees in the simulations (Fig. 4e–h), which account for a large proportion of total biomass (Fig. D11) and GPP (Fig. D12). This result supports hypothesis H3. For the second-order sensitivity, the contributions of most parametric interactions are relatively small (Fig. D13), as the first-order sensitivity accounts for a majority of the total variance in model outputs (Fig. 6e–h).

To understand how climate will impact sensitivity results, we also calculated the Spearman rank correlation coefficients between the first-order sensitivity index and the corresponding climate drivers. Our results show that the sensitivity of target carbon storage and maintenance respiration rate is negatively correlated with annual mean precipitation and relative humidity but is positively correlated with annual mean air temperature. This suggests that they are more important during the period of stressed conditions comprised of low precipitation, low humidity, and high temperature (Fig. 7). Sensitivity to the leaf allometry coefficient b is positively correlated with annual mean precipitation and relative humidity. This suggests the leaf carbon allocation is more important under favorable environmental conditions for growth. In general, our results suggests the climate has a larger impact on the parametric sensitivities for short-term carbon fluxes (GPP and NPP) and vegetation status (LAI) but has a smaller impact on parametric sensitivities for long-term vegetation carbon stocks.

Our bivariate spline analysis (Wahba, 1990) shows that, for V_c,max25 and target storage carbon, an increase in either of these parameters will cause an increase in the output of GPP, NPP, LAI, and biomass (Fig. 8). For the parameters related to leaf and stem allometry, however, the relations may differ depending on the output and the year of interest. At year 130, the higher leaf allocation normally leads to higher fluxes (NPP and GPP) but less biomass. Meanwhile, higher stem allocation leads to higher biomass but smaller fluxes (NPP and GPP). This suggests that the trade-offs between carbon allocation to stem vs. leaf tissues leads to a corresponding trade-off between carbon stocks and productivity in the model predictions.

Figure 5 Simulated temporal dynamics in basal area (BA, m² ha⁻¹) and the corresponding first-order parametric sensitivity indices. The left panels show the simulated ranges of BA for (a) all, (b) small (diameter <10 cm), (c) medium (10 cm < diameter <50 cm), and (d) large trees (diameter >50 cm). Right panels show the sensitivity for the top six most important parameters for (e) all, (f) small, (g) medium, and (h) large trees, in order of importance. See Fig. 2 for details on legends.

4 Discussion

4.1 Comparing parameter sensitivities to other models

While second-generation vegetation demographic models such as CLM4.5(FATES) provide new opportunities to predict the global carbon cycle, the larger number of parameters also creates challenges for identifying key processes for further investigation. In this study, we apply a global sensitivity analysis to determine the influential parameters over a specified region of the parameter space. So far, several uncertainty and sensitivity analyses have been conducted for size-structured land surface models (Pappas et al., 2013; LeBauer et al., 2013; Wang et al., 2013; Dietze et al., 2014; Collalti et al., 2019). In comparison with previous sensitivity analyses of size-structured models, our study considers a much larger number of parameters, i.e., >80 compared with ∼20–35 parameters (Pappas et al., 2013; LeBauer et al., 2013; Wang et al., 2013; Dietze et al., 2014), the difference in parametric sensitivity for different tree sizes, and the interactions among the key parameters. In general, our analysis shows similar results to sensitivity analysis on first-generation “big-leaf” vegetation models (e.g., Sargsyan et al., 2014), which show the importance of photosynthetic capacity, V_c,max25, for predicting GPP and NPP. However, we do show important parameters that are unique to LSMs with second-generation vegetation demography. Specifically, results shown here indicate the importance of leaf and stem allometry parameters, which control dynamic carbon allocation strategies based on size, and thus control the general vegetative state and size structure of the forest (Waring et al., 1998; Waring and Running, 2010). Importantly, a significant amount of variability in allometry is reported for different species and regions of the world (Feldpausch et al., 2011; Dietze et al., 2008), and thus it is critical to achieve an accurate parameterization of allometry for second-generation vegetation demographic models within LSMs. The importance of allometric parameters could result from the fact that the relationship between allometric coefficients and carbon allocation is highly nonlinear based on a power function (see Appendix A for details). This result is in agreement with the results from a recent sensitivity analysis study using a size-structured vegetation model (Collalti et al., 2019). Our sensitivity analysis also shows the importance of carbon storage for the prediction of mortality rate and thus the total biomass. This is in agreement with other sensitivity analyses of CLM which show the plant mortality rate as a key parameter for the prediction of total biomass (Sargsyan et al., 2014). However, we do want to point out that there could be a potential bias as carbon starvation is the main mechanism that kills trees in the simulations in our study; however, in reality, there could be many other causes of mortality such as wind, insects, and fire (see McDowell et al., 2018 for a review). The current implementation of hydraulic failure in CLM4.5(FATES) is only based on very low soil moisture thresholds and more mechanistic representation of plant hydrodynamics (e.g., Christoffersen et al., 2016) could result in the importance of hydraulic traits for mortality and vegetation dynamics. By exploring parametric sensitivity for small, medium, and large trees, we show that the ranking of parameter importance changes with the size of plants (e.g., Fig. 2). This result is in agreement with a recent study that showed the influence of certain functional traits varied with size (Falster et al., 2018).

In our analysis, we observed a number of key similarities in model response to parameter variations in photosynthetic capacity, mortality, and respiration parameters (Pappas et al., 2013; LeBauer et al., 2013; Wang et al., 2013; Dietze et al., 2014); however, there are differences in the order of parameter importance. For example, Dietze et al. (2014) showed that growth respiration fraction was the most important parameter for the simulation of NPP, and V_c,max25 only ranked as the seventh most important parameter. For our analysis, V_c,max25 and growth respiration fraction are the first and second most important parameters. This difference in parameter sensitivity rank may result from the fact that Dietze et al. (2014) used variable parameter ranges based on data (i.e., an uncertainty quantification study), while our sensitivity analysis uses equal percentage variations (see details in the discussion “limitation of methods” subsection). We also found that some parameters that are identified as important in other studies are not found to be important in our analysis. For example, Dietze et al. (2014) showed that water conductance that determines the upper boundary of transpiration is the second most important parameter for simulated NPP, but a similar parameter (smpso; Table D2) that defines soil water potential for opening stomata is not important in our analysis. This could be related to the fact that our site is much wetter than the temperate forests simulated by Dietze et al. (2014). Pappas et al. (2013) showed that the root distribution parameter that determines the fraction of fine roots in the upper soil layer is one of the top five parameters for the simulations of vegetation carbon fluxes and stocks; however, in our sensitivity analysis, the two root distribution parameters (root_a and root_b; Table D2) are not important for both vegetation carbon fluxes and stocks. This difference could also result from a wider range of variations ($\sim \pm \mathrm{30}$ %) in the study of Pappas et al. (2013) compared to our 15 % variations of the default parameters. Finally, our analysis shows the importance of allometry parameters, which are not considered in many previous studies (Pappas et al., 2013; LeBauer et al., 2013; Wang et al., 2013; Dietze et al., 2014).

Figure 6 Simulated temporal dynamics in carbon fluxes and stocks, and the corresponding first-order parametric sensitivity indices. The left panels show the simulated ranges for (a) GPP (kg C m⁻² yr⁻¹), (b) NPP (kg C m⁻² yr⁻¹), (c) LAI (m² m⁻²), and (d) total biomass (kg C m⁻²). Right panels show the sensitivity for the top six most important parameters for (e) all, (f) small, (g) medium, and (h) large trees, in order of importance. See Fig. 2 for details on legends.

4.2 Comparing simulations with observations

The goal of our study is not to reproduce the observations but instead to identify important parameters that can be better estimated for the model to fit observations. Thus, we lay out potential parameter estimation improvements to achieve this goal. We do want to highlight three caveats. First, improved estimation of the most sensitive parameters may not be most efficient if they have relatively small uncertainty or variability across different species and locations. Second, even if the estimates for most sensitive parameters are perfect, we may still not be able to fit model predictions to observations if there is deficiency in the representation of key processes in the model. Third, the recycled climate drivers from 1948 to 1972 may not match the observational periods. Given observation data limitations for our site, we conduct a qualitative comparison of our model simulations to ranges reported in the literature for the tropics. Not surprisingly, our model results show a variation of model–data mismatch for key vegetation states. For LAI (Fig. 6c), our simulated range is between ∼1.9 and 6.0 m² m⁻², which is lower than the observed range of ∼3.0–6.9 m² m⁻² based on LAI estimated from MODIS (Knyazikhin et al., 1999) during 2000–2016 within a 0.5^∘ window around our site. Our sensitivity analysis showed that leaf allometry coefficient b and target carbon storage are two key parameters for simulated LAI (Fig. 6g), and we expect that a better estimation of these parameters with data could potentially improve the model simulations. For GPP (Fig. 6a), the simulated range is between ∼1.0 and 3.0 kg C m⁻² yr⁻¹, which is also lower than the observed range of ∼2.4–3.7 kg C m⁻² yr⁻¹ based on extrapolation from eddy flux tower observations and climate during 1981–2010 (Jung et al., 2009). Our analysis suggests that photosynthetic capacity, as represented by V_c,max25, target carbon storage, and top-of-canopy specific leaf area, is an important parameter (Fig. 6a), and an improved estimation of them could help improve model simulations of GPP. We are not able to access on-site data for other model outputs. Therefore, we compare our model outputs with ranges from multiple tropical sites to evaluate their validity. For biomass (Fig. 6d), the simulated range of ∼2.5–12.5 kg C m⁻² is lower than the observed range of ∼7.3–21.3 kg C m⁻² from 21 transects within three tropical sites (Hunter et al., 2013). For BA (Fig. 5a), the simulated range of ∼5.0–30.0 m² ha⁻¹ is also lower than the observed range of ∼17.1–35.2 m² ha⁻¹ from five tropical sites (Hunter et al., 2013). Our results show that stem allometry coefficient c is the most important control on BA and biomass, and an improved parameterization on stem allometry could help improve the model simulations. For the DBH growth, there are large variances in the observed values across different sites with the range of 0–3 cm yr⁻¹ (Lieberman et al., 1985; Worbes, 1999; Adams et al., 2014). The simulated average DBH growth is between 0 and 0.4 cm yr⁻¹ but could be as high as 4 cm yr⁻¹ for medium and large trees (Fig. 2). Based on our sensitivity analysis (Fig. 2), we expected an improved parameterization of both allometry coefficient c and target carbon storage could help fit the model predictions to data.

Figure 7 Spearman's correlation coefficients between climate drivers and six most important parameters identified for (a) GPP, (b) NPP, (c) LAI, and (d) biomass.

We compare our mortality simulations with an extensive dataset of observed mortality of 1781 species from 14 pan-tropical large-area ForestGEO forest dynamics plots (Johnson et al., 2018). For this study, the forest plots ranged from 2 to 52 ha each, with 371 ha in total, in which all recorded stems are ≥1 cm diameter at breast height. Our comparison shows that the CLM4.5(FATES) simulations of medium and large tree mortality (Fig. 4c, d) are close to the 95 % confidence interval of observed values, which is about ∼0.5 %–5.7 % per year. However, for small and medium trees (Fig. 4b, c), the simulated mortality rate of ∼15 %–30 % and ∼1 %–10 % is high when compared to the observed 95 % confidence interval of mortality rate of ∼0.6 %–11.3 % and ∼0.8 %–3.0 % for small and medium trees, respectively. The high predicted mortality rate of small trees could result from the fact that the model predicts a very high mortality rate for very small trees (<1 cm), as they cannot survive after establishment due to low light conditions in the simulations. Since the small trees have such a large fraction of the population in our simulations (Fig. 3b), the overall mortality rate (Fig. 4a) of ∼15 %–30 % is also high when compared to observations (∼0.6 %–11.3 %); however, if we separate the mortality rate of very small trees from the calculation of the overall mortality, then the simulated mortality rates of 1 %–10 % (Fig. D14b) are in the range of observations. The very high mortality rate range of smaller trees (∼10 %–30 %; Fig. D14a) spans the reported seedling/sapling mortality rate, e.g., ∼15 %–21 % per year from 1- to 20-year-old tropical forest stands in Costa Rica (Dupuy and Chazdon, 2006). However, there is potential for improvement for site-level simulations as the current recruitment algorithm within CLM4.5(FATES) depends only on the availability of seed bank but not on the density, light, and water availability. The relatively high mortality rate of small and medium trees could also be linked to the fact that CLM4.5(FATES) uses the perfect plasticity approximation (PPA) to simulate the canopy light availability for understory trees (Fisher et al., 2018), which may create canopy closure too fast for the small- and medium-sized trees to survive under low light conditions. We expect that future improvements in recruitment and representation of the light environment within the PPA could be helpful for a better prediction of tree mortality for small- and medium-sized trees. It is also possible that the observed mortality rate for small trees could potentially be underestimated if all the trees in a certain size classes died at a shorter time frame than the census intervals. Our sensitivity analysis indicates that key model parameters that can be better estimated for improved mortality predictions include stem allometry parameters, V_c,max25, target carbon storage, and mortality rate under stress (Fig. 4).

Another reason for the data–model discrepancy could result from the limited representation of diverse tropical species or traits with the simulation of a single PFT. This is a limitation of many LSMs, as they typically only have two to three PFTs for tropical forests (e.g., only evergreen and deciduous for tropical trees within CLM). CLM4.5(FATES) has the potential to better represent the trait diversity through trait filtering under different environmental conditions (Fisher et al., 2010, 2015). One critical component to incorporate traits into the model is to represent the trait trade-off and coordination for different PFTs. Through our sensitivity analysis, we have identified key parameters for vegetation dynamics, which can be targeted for the representation of trait trade-off and coordination in the tropics. For example, our study shows that a higher stem carbon allocation could reduce the GPP and a higher V_c,max25 could increase GPP (Fig. 8). The potential exploration of trade-off and coordination between these two parameters could be critical to resolve different PFTs and represent the trait variations. Even though the simulated ranges of the model outputs are different than the observations, our sensitivity analysis should still be valid in view that a primary end goal of this research is to identify important parameters that can be better estimated for the model to better fit observations. For example, Holm et al. (2019) utilized results from our study to implement their tropical forest parameterization, specifically by increasing their target carbon storage parameter to obtain higher survival and lower growth.

In addition to directly comparing the model outputs to observations, we want to highlight that the sensitivity analysis will also allow us to explore the functional relationships between model parameters and outputs. Future synthesis studies that show these functional relationships using data across different sites could be very useful to evaluate the fidelity of model structure to represent the key processes that control these relationships.

Figure 8 Relations between outputs of CLM4.5(FATES), including GPP, NPP, LAI, and biomass (units shown in Fig. 6), and the most sensitive parameters, i.e., V_c,max25 (units are $\mathrm{\mu }\mathrm{mol}{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), target storage carbon (unit is the ratio of leaf biomass), leaf and stem allocation (unitless parameters) for simulation year 10 (red) and 130 (blue). Shown are the mean relations with the 95 % confidence intervals in grey envelopes. These figures show how an output will generally increase or decrease when a given parameter is changed.

4.3 Limitation of methods

Our study is the first global sensitivity analysis for CLM4.5(FATES); however, it is subjected to several limitations that could be improved for future studies. First, our study uses an arbitrary choice of parameter ranges (±15 %), which determines the variance in the model outputs and the corresponding results of the sensitivity analysis. However, we expect that our analysis can reveal the importance of parameters given equal percentage of variations, which can help us gain a better understanding of the model structure. We do acknowledge that uncertainty analysis studies that specifically consider the potential ranges of values in tropical forests based on observations could provide insights on which additional measurements are needed to explain variance in the model prediction.

Second, we only consider the correlation in pairs of parameters that determine temperature responses for deactivation energy and entropy in photosynthesis. We do want to point out that the potential correlation among other parameters (Díaz et al., 2016), such as the trade-off between mortality and growth parameters and the correlation among coefficients for allometric equations, could affect the simulated model output ranges and the sensitivity results. However, our exploration of parameter sensitivity assuming their independence could still help us understand the baseline parameter control on model behaviors (Xu and Gertner, 2009). The exploration of trade-offs and coordination among different parameters requires data analysis for multiple traits of the same species. The Predictive Ecosystem Analyzer (PEcAn) framework (LeBauer et al., 2013) could be a useful tool to synthesize plant trait data to estimate model parameter distributions. The challenge is that, even though there are great efforts in the research community to compile plant trait data across the globe (Kattge et al., 2011a, b), there are still limited datasets with observations of multiple traits for the same species. Future uncertainty analysis studies that explicitly consider the prior distributions and correlations for all the parameters can build on this analysis and gain further insights on where the uncertainty in the model predictions comes from.

Third, it is possible that the parameter sensitivity could be different if we use different model inputs, different sites, and different structures of subcomponents within the model. For example, using site-level climate drivers, instead of the reanalysis meteorological drivers used in this study (Qian et al., 2006), could lead to different sensitivity values since our analysis showed that parameter importance is quite sensitive to different climate conditions. Furthermore, there are ongoing development activities to improve different components of the models. For example, there are current efforts to incorporate different representations of tree allometry within CLM4.5(FATES), which have different formulations between size and biomass, e.g., Chave et al. (2014), or the current formulation of the photosynthetic process in the CLM4.5(FATES) can be replaced with a model that more accurately represents the allocation of nitrogen and thus the photosynthetic process (see Xu et al., 2013; Ali et al., 2016). Therefore, model improvements such as these can affect corresponding sensitivity analysis results. To understand the impact of site-level variations on model dynamics, similar sensitivity analysis across different sites can be conducted to understand how climate variability will affect the sensitivity analysis results.

Finally, although the FAST method presented in this study can provide a comprehensive analysis of the parameters that control vegetation demography, it is mostly built on statistical relationships (e.g., Fig. 8). A complementary approach is to use more tractable theoretical ecology models (e.g., Farrior et al., 2016 and Falster et al., 2018) to approximate the underlying model input–output relationships, which can provide more mechanistic understanding of model behaviors.

5 Conclusions

LSMs have many parameters that could potentially affect the outcome of their simulations. In this study, we use the FAST analysis to conduct a high-dimensional global sensitivity analysis on CLM4.5(FATES). We use an intermediate complexity of simulation: runs are sufficiently long to permit short-term physiological variance to propagate into the long-term forest demographic structure. Even though we do not explore competitive dynamics between different PFTs, our sensitivity analysis will guide us on the selection of key plant traits for the consideration of trait trade-off and coordination in order to improve PFT coexistence within CLM4.5(FATES).

Our analyses show that the target carbon storage and stem allometry parameters are important for the simulation of DBH growth for individual trees and tree mortality. The photosynthetic parameter,V_c,max25, is the most important for the simulation of carbon fluxes including GPP and NPP. The combination of stem allometry, target carbon storage, and V_c,max25 dominantly control the simulation of total BA and long-term carbon stocks. These identified growth and survival parameters will help us better understand the key control of fast and slow carbon and vegetation dynamics within the next generation of demographically enabled LSMs toward improved model parameterization and model structure.

The results of the sensitivity analysis presented here can be utilized to construct the parameter-output response surface for the CLM4.5(FATES) model, which can assist future efforts for model calibration or diagnosis. These findings may help us better understand the overall model structure and guide the estimation of key model parameters with significant control over vegetative processes in these models for better model fitting to data. The FAST analysis provides a promising means of analyzing complex LSM components and can be a powerful tool in understanding the necessarily high-dimensional representation of living systems within Earth system models.

Appendix A: Allometry equations

The following equations are cohort-based calculations for allometry in CLM4.5(FATES). Interested readers are referred to Fisher et al. (2015) for more information. The parameters used for the allometry equations include dbh2h_m, dbh2h_c, dbh2bd_a, dbh2bd_b, dbh2bd_c, and dbh2bd_d (all are unitless variables). Specifically, the dead wood biomass (BD; kg C) is calculated as a function of diameter (DBH; cm), height (h; meter), and wood density (g cm⁻³):

$$\begin{array}{}\text{(A1)}& \mathrm{BD}=\left({\text{dbh2bd}}_a\right)\left(h^{{\text{dbh2bd}}_b}\right)\left({\text{DBH}}^{{\text{dbh2bd}}_c}\right)\left({\text{density}}_{\mathrm{wood}}^{{\text{dbh2bd}}_d}\right).\end{array}$$

The height (m) is calculated based on DBH (cm) as follows:

$$\begin{array}{}\text{(A2)}& H={\mathrm{10}}^{\mathrm{dbh}\mathrm{2}{\mathrm{h}}_c}\left({\mathrm{DBH}}^{\mathrm{dbh}\mathrm{2}{\mathrm{bd}}_m}\right).\end{array}$$

Appendix B: Temperature response curve

The parameters used for the temperature response curve equations include the equation to calculate the maximum carboxylation rate, V_c,max25, the maximum electron transport rate, J_max, and the triose phosphate use (TPU) limited carboxylation rate, TPU (also all parameters here are unitless) (Fisher et al., 2015). The temperature response equations for $V_{\mathrm{c},\max ,z}$, J_max,z, and TPU_z are

$$\begin{array}{}\text{(B1)}& {\displaystyle }\begin{array}{rl}{V}_{\mathrm{c},\max ,z}=& V_{\mathrm{c},\max ,\mathrm{25}}\left(e^{\frac{\mathrm{vcmaxha}}{\left(\mathrm{0.001}\mathrm{rgas}\right)(t_{\mathrm{frz}}+\mathrm{25})}}\right)(\mathrm{1}-{\displaystyle \frac{t_{\mathrm{frz}}+\mathrm{25}}{t_{\mathrm{veg}}}})\\ & \left({\displaystyle \frac{\mathrm{vcmaxc}}{\mathrm{1}+e^{-\mathrm{vcmaxhd}+\left(\mathrm{vcmaxse}\right)\left(t_{\mathrm{veg}}\right)}}}\right)\end{array}\text{(B2)}& {\displaystyle }\begin{array}{rl}{J}_{\max ,z}=& J_{\max ,\mathrm{25}}\left(e^{\frac{\mathrm{jmaxha}}{\left(\mathrm{0.001}\mathrm{rgas}\right)(t_{\mathrm{frz}}+\mathrm{25})}}\right)(\mathrm{1}-{\displaystyle \frac{t_{\mathrm{frz}}+\mathrm{25}}{t_{\mathrm{veg}}}})\\ & \left({\displaystyle \frac{\mathrm{jmaxc}}{\mathrm{1}+e^{-\mathrm{jmaxhd}+\left(\mathrm{jmaxse}\right)\left(t_{\mathrm{veg}}\right)}}}\right)\end{array}\text{(B3)}& {\displaystyle }\begin{array}{rl}{\mathrm{TPU}}_z=& \mathrm{tpu}\mathrm{25}\left(e^{\frac{\mathrm{tpuha}}{\left(\mathrm{0.001}\mathrm{rgas}\right)(t_{\mathrm{frz}}+\mathrm{25})}}\right)(\mathrm{1}-{\displaystyle \frac{t_{\mathrm{frz}}+\mathrm{25}}{t_{\mathrm{veg}}}})\\ & \left({\displaystyle \frac{\mathrm{tpuc}}{\mathrm{1}+e^{-\mathrm{tpuhd}+\left(\mathrm{tpuse}\right)\left(t_{\mathrm{veg}}\right)}}}\right),\end{array}\end{array}$$

where t_frz is the freezing point of water in Kelvin (273.15 K).

Appendix C: Carbon storage in CLM4.5(FATES)

The target carbon storage is the cushion parameter shown in Table D3. Specifically, a higher value of this parameter will lead to a higher allocation of carbon to storage and thus a lower allocation to growth at the specific time step. Also, carbon storage plays an important role for the simulated mortality through the parameter that controls the mortality rate under stress, stress_mort in Table D3. The tree will be under stress when it has low carbon storage (< leaf biomass). Therefore, the target carbon storage parameter and the mortality rate under stress parameter play a large role in determining the level of mortality that occurs in the simulations.

Carbon storage, b_store (in kg C/cohort), plays a very important role in both growth and mortality (Fisher et al., 2015). Specifically, CLM4.5(FATES) assumes a target carbon storage determined by the multiplication of leaf biomass (b_leaf) and the target carbon storage parameter (i.e., the target amount of carbon plants store relative to leaf biomass; S_cushion, variable cushion in Table D3). At the specific time, the carbon balance for growth and storage is calculated as follows:

$$\begin{array}{}\text{(C1)}& C=\mathrm{NPP}-T_{\mathrm{md}}{f}_{\mathrm{md},\min },\end{array}$$

where T_md is the maintenance respiration and f_md,min is the minimum fraction of the maintenance demand (storage priority parameter in Table D1) that the plant must meet each time step, which represents a life-history-strategy decision concerning whether leaves should remain on in the case of low carbon uptake (a risky strategy) or not be replaced (a conservative strategy).

The fraction of the carbon balance for each cohort allocated to the carbon storage pool (f_store) will be determined by the following equations:

$$\begin{array}{}\text{(C2)}& f_{\mathrm{store}}=e^{(-f_{\mathrm{tstore}}{)}^{\mathrm{4}}},\end{array}$$

where

$$\begin{array}{}\text{(C3)}& f_{\mathrm{tstore}}=\max (\mathrm{0},{\displaystyle \frac{b_{\mathrm{store}}}{S_{\mathrm{cushion}}{b}_{\mathrm{leaf}}}}).\end{array}$$

Thus, the target carbon storage parameter, S_cushion, can affect carbon allocations. Specifically, a higher value of S_cushion will lead to a higher allocation of carbon to storage and thus lower allocation to growth at the specific time step.

Carbon storage also plays an important role for the mortality. Specifically, carbon starvation mortality (M_cs) is calculated as follows:

$$\begin{array}{}\text{(C4)}& M_{\mathrm{cs}}=S_{\mathrm{m}}\max (\mathrm{0},\mathrm{1}-{\displaystyle \frac{b_{\mathrm{store}}}{b_{\mathrm{leaf}}}}),\end{array}$$

where S_m is the stress mortality factor (i.e., stress_mort in Table D3).

Appendix D: Appendix figures and tables

Figure D1 Comparison of first-order parametric sensitivity for medium (10 cm < diameter <50 cm) tree density averaged over 20, 30, and 40 years.

Figure D2 Simulated total change in diameter at breast height (dDBH) from CLM4.5(FATES) for all trees and its fractional distribution for small (diameter <10 cm), medium (10 cm < diameter <50 cm), and large trees (diameter >50 cm). Shown is the mean simulation (black line) with 95 % spread of the simulation ensemble.

Figure D3 Impacts of stem allometry on the change in diameter at breast height (dDBH) averaged over the simulation years 100–130 for trees of different sizes. The shaded area shows the 95 % confidence interval of these relations.

Figure D4 Mortality outputs from CLM4.5(FATES), including the mechanisms of M1 – background mortality, M2 – hydraulic failure, M3 – carbon starvation, and M4 – impact mortality. The bottom panel shows the total mortality, which is the sum of M1–M4. Shown is the 95 % (light grey) spread of the simulation ensemble, along with the mean simulation (black lines).

Figure D5 Impacts of minimum crown spread on the basal area (BA) averaged over the simulation years 100–130 for trees of different sizes. The shaded area shows the 95 % confidence interval of these relations.

Figure D6 First-order parametric sensitivity indices for tree density of all trees (a) and the corresponding fraction of (b) small (diameter <10 cm), (c) medium (10 cm < diameter <50 cm), and (d) large trees (diameter >50 cm). See Fig. 2 for details on legends.

Figure D7 Second-order sensitivity index of the model parameters for the basal area (BA) outputs from CLM4.5(FATES) for (a) all trees, (b) small trees, (c) medium trees, and (d) large trees. Shown are the top eight most important parameter interactions in order of importance based on the mean parametric sensitivity across years (red is the most important and blue is the least important)

Figure D8 Second-order sensitivity index of the model parameters for the change in diameter at breast height (dDBH) outputs from CLM4.5(FATES) for (a) all trees, (b) small trees, (c) medium trees, and (d) large trees. Shown are the top eight most important parameter interactions in order of importance based on the mean parametric sensitivity across years (red is the most important and blue is the least important).

Figure D9 Second-order sensitivity index of the model parameters for the mortality outputs from CLM4.5(FATES) for (a) all trees, (b) small trees, (c) medium trees, and (d) large trees. Shown are the top eight most important parameter interactions in order of importance based on the mean parametric sensitivity across years (red is the most important and blue is the least important).

Figure D10 Comparisons of parametric sensitivities and the corresponding 95 % confidence intervals for different model outputs at year 130. Shown are the identified four most important parameters for (a) GPP, (b) NPP, (c) LAI, and (d) total biomass.

Figure D11 Fraction of total biomass for trees of different sizes, including small (diameter <10 cm), medium (10 cm < diameter <50 cm), and large trees (diameter >50 cm).

Figure D12 Fraction of total GPP for trees of different sizes, including small (diameter <10 cm), medium (10 cm < diameter <50 cm), and large trees (diameter >50 cm).

Figure D13 Second-order sensitivity index of the model parameters for the GPP, NPP, LAI, and biomass outputs from CLM4.5(FATES). Shown are the top eight most important parameter interactions in order of importance based on the mean parametric sensitivity across years (red is the most important and blue is the least important).

Figure D14 Mortality outputs from CLM4.5(FATES) for trees with DBH smaller than 5 cm (a) and all trees with DBH between 5 and 10 cm (b). Shown is the 95 % (light grey) spread of the simulation ensemble, along with the mean simulation (black lines).

Table D1 Parameter sets used in this study – part 1.

Table D2 Parameter sets used in this study – part 2.

Table D3 Parameter sets used in this study – part 3.

Table D4 Parameter sets used in this study – part 4.