Parametric Sensitivity of Vegetation Dynamics in the TRIFFID Model and the Associated Uncertainty in Projected Climate Change Impacts on Western U.S. Forests

Changing climate conditions impact ecosystem dynamics and have local to global impacts on water and carbon cycles. Many processes in dynamic vegetation models (DVMs) are parameterized, and the unknown/unknowable parameter values introduce uncertainty that has rarely been quantified in projections of forced changes. In this study, we identify processes and parameters that introduce the largest uncertainties in the vegetation state simulated by the DVM Top‐down Representation of Interactive Foliage and Flora Including Dynamics (TRIFFID) coupled to a regional climate model. We adjust parameters simultaneously in an ensemble of equilibrium vegetation simulations and use statistical emulation to explore sensitivities to, and interactions among, parameters. We find that vegetation distribution is most sensitive to parameters related to carbon allocation and competition. Using a suite of statistical emulators, we identify regions of parameter space that reduce the error in modeled forest cover by 31±9%. We then generate large initial atmospheric condition ensembles with 10 improved DVM parameterizations under preindustrial, contemporary, and future climate conditions to assess uncertainty in the forced response due to parameterization. We find that while most parameterizations agree on the direction of future vegetation transitions in the western United States, the magnitude varies considerably: for example, in the northwest coast the expansion of broadleaf trees and corresponding decline of needleleaf trees ranges from 4 to 28% across 10 DVM parameterizations under projected future climate conditions. We demonstrate that model parameterization contributes to uncertainty in vegetation transition and carbon cycle feedback under nonstationary climate conditions, which has important implications for carbon stocks, ecosystem services, and climate feedback.

ecosystem. As plant types better suited for the new climate conditions begin to thrive, they can out-compete the historical vegetation. This transformation can have impacts on the local water cycle and, if changes are large enough, on the global carbon cycle. Dynamic vegetation models (DVMs) are used to analyze the impacts of climate on vegetation and project future changes. DVMs are used across scales ranging from watersheds, informing management decisions, to the globe as an important component of the global carbon cycle. The mathematical representations of processes in DVMs are imperfect and introduce uncertainty in model projections. In this study we identify processes that introduce the largest uncertainties in the vegetation distribution modeled by a DVM. We then identify settings of model parameters that improve the modeled historical vegetation when compared to observations. While unique model parameterizations simulated the historical vegetation well, each parameterization simulated different future vegetation transitions under future climate conditions. Our study shows that model parameterization introduces uncertainty in the projected future changes in vegetation. Quantifying and appreciating this uncertainty are important when considering management decisions based on projected vegetation changes at the local to global scale.

Introduction
Global climate change is increasing the vulnerability of some forested ecosystems to tree mortality (Allen et al., 2010;Anderegg et al., 2013). The western United States experiences seasonal droughts that, in the future, are projected to occur at higher temperatures (Ault et al., 2014;Diffenbaugh et al., 2015;Overpeck, 2013). Hotter droughts increase forest vulnerability to tree mortality Breshears et al., 2013;Williams et al., 2013), insect attack, and disease (Meddens et al., 2012). Regional-scale forest dieback has already been observed (Breshears et al., 2005;van Mantgem et al., 2009), and tree mortality is projected to continue under changing climate conditions (Adams et al., 2009;Allen et al., 2010).
Vegetation models are commonly used to understand how climate change may affect forest health and vegetation transitions (e.g., Coops & Waring, 2011;Halofsky et al., 2013;Lenihan et al., 2008). Dynamic vegetation models (DVMs) in particular simulate establishment, growth, competition, and mortality, making them a useful tool for understanding how vegetation species distributions and associated carbon stocks respond to climate forcings. DVMs, being process-based, can be applied at global (Friedlingstein et al., 2006;Sitch et al., 2008) to local scales  and have been used to project changes in vegetation distributions and associated carbon consequences in the western United States Jiang et al., 2012;Kerns et al., 2018;Rehfeldt et al., 2012;Shafer et al., 2015).
Many processes in DVMs must be parameterized resulting in numerous uncertain parameters. Parameter settings are determined based on observational measurements whenever possible (e.g., Bonan & Levis, 2006;Hickler et al., 2012;Krinner et al., 2005;Sitch et al., 2003), but parameters often do not explicitly represent physical, measurable processes. DVMs are commonly tuned for the region of interest by adjusting parameters one-at-a-time until (quasi-) steady-state simulations (often spanning thousands of model years) agree with observations from flux tower, forest inventory, or remote sensing data . Better agreement with observations can be achieved through compensating errors; therefore, many different model parameterizations may appear to improve performance while having very different values for individual parameters, highlighting the importance of systematic characterization of the sensitivity of model behavior to parameter settings. This "parametric uncertainty" is one major source of uncertainty in earth system model predictions, the others being uncertainties due to internal variability (i.e., natural fluctuations), imperfect model structure (i.e., imperfect mathematical representation or omission of earth system processes), and forcing scenario (Hawkins & Sutton, 2009). Many sensitivity analyses have been performed to quantify parametric uncertainty in global climate models (e.g., Stainforth et al., 2005;Murphy et al., 2004;Murphy et al., 2007;Collins et al., 2011;Sanderson et al., 2008;Williamson et al., 2013;Yamazaki et al., 2013;Qian et al., 2015;Li et al., 2019). Sensitivity studies have been performed that explored parametric uncertainty in DVMs (e.g., Fisher et al., 2010;Pappas et al., 2013;Ricciuto et al., 2018) but are commonly done offline (i.e., not coupled to the atmosphere). Parameters often represent processes related to nonlinear interactions between the land-surface and the atmosphere, such as photosynthesis, but the offline approach does not capture complex interactions and feedback between the land surface and the atmosphere. For example, if a meadow is overtaken by a forest, the local air temperature and humidity will be affected by the change in albedo and energy partitioning between sensible and latent heat; however, the driving climate data is not affected by the change in land surface cover leading to an unphysical mismatch in the conservation of mass (water/carbon), energy, and momentum. Coupling highly nonlinear systems that were parameterized and tested uncoupled has been shown to produce challenges (Bonan & Levis, 2006). Accounting for land-atmosphere interactions in sensitivity analyses captures the amplifying or dampening effects of parameter adjustments on land-atmosphere feedback, which is particularly important when the model is intended for making future projections. Parametric uncertainty may contribute considerably to estimates of future vegetation transitions and the evolution of the carbon cycle (Booth et al., 2012;Huntingford et al., 2009).
One approach for assessing the parametric uncertainty of complex systems with dynamical feedback is through the use of perturbed parameter ensembles (PPEs) where multiple parameters are adjusted simultaneously and the high-dimensional parameter space is sampled efficiently. Parametric uncertainty lies within the structural limitations of the model but has been shown to generate a wide range of model sensitivities in response to forcings (Stainforth et al., 2005). In addition to characterizing parametric uncertainty, such sensitivity analyses can inform model developers about structural uncertainty by identifying key processes influencing model behavior (Pappas et al., 2013).
The objectives of this study were to (1) identify parameters and processes that strongly influence simulated terrestrial ecosystem dynamics in a fully coupled DVM, (2) demonstrate an effective approach for assessing parametric uncertainty and improving model performance through parameter adjustments in DVMs, and (3) assess uncertainty in vegetation state changes under nonstationary climate conditions among plausible model parameterizations.
We assessed the sensitivity of vegetation state variables to model parameter settings in regional climate model simulations coupled with dynamic vegetation. We took an iterative approach starting with a oneat-a-time parameter sensitivity experiment to identify influential parameters in a DVM. We then simultaneously varied the influential parameters in a large PPE of model simulations. We used statistical emulators as a computationally cost-effective proxy of the DVM to estimate responses for the unsampled parameter space. We used a suite of emulators of vegetation properties (e.g., fractional coverage and biomass) to explore the sensitivity of vegetation dynamics to parameter values and parameter interactions. Through this process, we identified regions of parameter space that improved model performance in our chosen evaluation criteria. We selected an ensemble of model parameterizations that satisfied evaluation criteria and ran simulations under preindustrial, contemporary, and future climate to assess the uncertainty in forest coverage attributable to differences in model parameterization.

Model Description
Simulations of western U.S. climate and vegetation were generated on the Weather@home distributed volunteer climate modeling system (Guillod et al., 2018;Massey et al., 2015;Mote et al., 2015). The vast computational capacity of Weather@home enables large ensembles of short simulations to be generated rapidly lending to its popularity as an event attribution tool (e.g., Mera et al., 2015;Mote et al., 2016), but it has also been used to make statistically robust projections of regional climate change . Weather@home uses the UK Met Office's regional atmospheric model HadRM3P (0.22°longitude by 0.22°latitude horizontal resolution,~25km; 5-min time step) nested within the global atmospheric model HadAM3P (1.875°longitude by 1.25°latitude; 19 vertical levels; 15-min time step) with one-way daily coupling. Both global and regional models use the Met Office Surface Exchange Scheme, version 2 (MOSES2; Essery et al., 2001), which simulates the exchange of heat, momentum, water, and carbon between the land-surface and the atmosphere. Soil information was updated in HadRM3P to be consistent with the North American Land Data Assimilation Phase 2 (NLDAS-2; Xia et al., 2012; see the supporting information for details). For this application we used the western U.S. region (27-54°N and 218-254°E) with an alternative parameterization of HadAM3P/RM3P (Table S1 in the supporting information) selected from Li et al. (2019) that markedly reduced summer warm and dry biases over the western US.
Vegetation dynamics were simulated by the Top-down Representation of Interactive Foliage and Flora Including Dynamics model (TRIFFID; Cox 2001). TRIFFID uses CO 2 fluxes at the land-atmosphere interface to update plant distributions and soil carbon. Environmental conditions in the atmospheric and land surface models are translated into net primary productivity (NPP) rates (carbon assimilation). TRIFFID allocates NPP aggregated over the coupling period into the growth and expansion of five plant functional types (PFTs): broadleaf (BL) trees, needleleaf (NL) trees, C3 grass, C4 grass, and shrubs (SHs). Competition among PFTs occurs based on a tree-SH-grass hierarchy, with competition occurring between the two tree and grass functional types based on canopy height. Canopy height, roughness length, albedo, leaf area index (LAI), surface resistance, canopy heat capacity, canopy catchment capacity, and root density are calculated based on the new vegetation distribution and fed back to the land surface model. This approach allows climateinduced changes in biophysical properties of the land surface to feed back on to the atmosphere, altering surface-atmosphere fluxes of heat, water, momentum, and CO 2 .
TRIFFID couples a photosynthesis-stomatal conductance model (Cox et al., 1998) to an ecological population model based on the Lotka-Volterra competition equations (Cox 2001). The vegetated state is described by three primary diagnostics: fractional coverage, LAI, and canopy height of PFTs. TRIFFID updates the vegetation state using a carbon balance approach allocating accumulated carbon (NPP) from MOSES into growth (biomass) and expansion (fractional coverage). Allocation and competition are highly parameterized and dependent on the primary variables calculated by TRIFFID. Allocation is related to the phenological status of LAI, which is updated by temperature. Competition acts to reduce growth through density-dependent litter production and limits expansion based on relative canopy heights between competing species. Competition for light, water, and nutrients is not explicitly represented; resources are initially accessible to all PFTs. The canopy height of the PFT best suited for the local climatic envelope increases and that PFT begins to out-compete the surrounding vegetation. This approach approximates the shading out of surrounding vegetation, limiting growth by increasing litterfall. The updated canopy height and LAI are used to calculate the biophysical variables required by MOSES. Canopy height determines the aerodynamic roughness lengths used in the calculation of surface-atmosphere fluxes of heat, water, momentum and CO 2 . Water freely available for evaporation is dependent on the canopy catchment capacity, which varies linearly with LAI. Albedo is calculated from PFT-dependent parameters. TRIFFID is limited by a onepool soil decomposition model, which has been shown to enhance soil respiration and amplify the carbon-climate feedback (Jones et al., 2005). Mortality is represented by a constant background rate that varies by PFT and occurs evenly across space and time.
TRIFFID was built with an "equilibrium mode" in which an implicit time-stepping scheme is employed to bring the model into equilibrium in relatively few simulation years. In equilibrium mode, the climate model runs with stationary vegetation for a user-defined number of years, over which NPP is accumulated. TRIFFID then runs uncoupled from the climate model for 10,000 iterations, recycling the accumulated NPP as input for each iteration, before feeding state variables back to the climate model. TRIFFID has been evaluated and applied across a range of spatial scales (e.g., Cox et al., 2000;Sitch et al., 2008) and across ecosystems (e.g., Cox et al., 2004;Piao et al., 2009). TRIFFID was used in the Coupled Climate-Carbon Cycle Model Intercomparison Project (C4MIP; Friedlingstein et al., 2006) and is incorporated into the Hadley Centre Global Environmental Model version 2 (HadGEM2; Martin et al., 2011) and HadGEM3 (Walters et al., 2011) family of earth system models used in the fifth and sixth phases of the Coupled Model Intercomparison Project, respectively. TRIFFID and MOSES form the foundation of the Joint UK Land Environment Simulator (JULES; Clark et al., 2011) used in the most recent assessment of the global carbon budget (Le Quéré et al., 2016).

Model Simulations
Atmospheric forcings, including greenhouse gas concentrations, volcanic aerosol optical depth, ozone concentrations, sulfur dioxide emissions, and solar irradiance, followed those in the Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor et al., 2011) historical experiment for simulations including years before 2006 and the RCP 8.5 experiment for years after 2005. HadAM3P/HadRM3P were run here with prescribed sea surface temperatures (SSTs) and sea ice concentrations (SICE). Simulations before 1985 used the HadISST observational SST and SICE dataset version 2.1.0.0 (Titchner & Rayner, 2014) as in Guillod et al. (2018). The Operational Sea Surface Temperature and Sea Ice Analysis observational product (OSTIA; Donlon et al., 2012) was used in simulations between 1986 and 2015. Future SSTs were created by applying a delta calculated from CESM1-CAM5 to the OSTIA observational SSTs, similarly to . The delta was calculated by taking the CMIP5 CESM1-CAM5 average of three ensemble members (r1i1p1-r3i1p1), concatenating the historical experiment with the RCP 8.5 experiment, Gaussian smoothing with an 11 year window (standard deviation of 5.7 years), and taking the monthly difference of the smoothed SSTs between January 2046 and December 2075 and January 1986 and December 2015 (e.g., January 2046 minus January 1986). Future SICE was created with a similar methodology, detailed in the Supplementary Information.
We initialized simulations from an existing spinup period  and performed an additional four years of spinup (1900)(1901)(1902)(1903) to initialize the preindustrial simulations, in order to bring the deep soil temperature and moisture into equilibrium with the preindustrial climate state. Each simulation was initialized with the same state variables but with slight perturbations to the atmosphere's potential temperature field to introduce atmospheric initial condition (IC) variability .
The assessment of parameter uncertainty in DVMs is challenged by the large computational costs of bringing multiple parameterizations of a DVM into equilibrium with the climate. The long residence time of carbon in soils and vegetation can require thousands of model years for carbon pools to reach equilibrium. A primary strength of Weather@home2 is the ability to rapidly generate many short simulations; however, it is ill-suited to perform very long simulations; thus, using Weather@home2 for assessing parameter uncertainty in DVMs requires nonstandard techniques. We utilized TRIFFID's equilibrium mode (Cox 2001) to bring the vegetation into quasi-equilibrium with the climate in ensembles of short simulations. For the parameter sensitivity experiments (Table 1; Steps 1-3) we ran 4-year simulations with two rounds of TRIFFID equilibrium iterations taking place after the second and fourth years. NPP was accumulated over the first 2 years and recycled by TRIFFID for 10,000 iterations, the new vegetation state variables were fed back to the climate model, and another 2 years of NPP was accumulated and again recycled for 10,000 iterations in a second equilibrium round. We introduced variability in the climate conditions, and thus the accumulated NPP, through initial atmospheric condition ensembles and replicated this process for many different model parameterizations. For the production runs (Table 1) we ran 5-year simulations with one round of TRIFFID equilibrium iterations taking place after the fifth year.

Evaluation Data
Evaluating DVMs against satellite observations presents challenges where disturbances (natural and anthropogenic) have recently occurred, because the vegetation may be at various stages along its successional pathway. For this reason DVMs are commonly evaluated by comparing the simulated vegetation distribution after spinup to the potential climax vegetation (Bachelet, Ferschweiler, et al., 2015;Shafer et al., 2015). Küchler (1965) mapped and defined the potential natural vegetation that would occur without human interference at the climax of succession. We reclassified potential climax vegetation types from Küchler (1965) into NL trees, BL trees, grasses, and SHs (Table S2) in accordance with the Vegetation-Ecosystem Modeling and Analysis Project (VEMAP; VEMAP Members, 1995) and determined the dominant vegetation type in each grid cell. Grid cells dominated by SHs or grasses were grouped together as recommended by Kelley et al. (2013).
Evaluation of dominant vegetation types emphasizes that the correct PFT dominates in the correct climatic zones of the domain. Other DVM studies use a similar metric for model evaluation (Bachelet, Ferschweiler, et al., 2015;Shafer et al., 2015). In addition, we evaluated the contemporary forest fractional cover against an observational data set derived from satellite and Forest Inventory Analysis (FIA) data (Ruefenacht et al., 2008). We evaluated the forest biomass against a forest inventory based data set (Wilson et al., 2013) and the canopy height against a satellite derived data set for the year 2005 (Simard et al., 2011), regridded to the HadRM3P model grid using inverse distance weighting. Grid cells where the landsurface type occupying the largest area within a grid cell was pastureland, cropland, or urbanland, as defined by Meiyappan and Jain (2012), were masked out. The domain was divided into level I ecoregions (Omernik & Griffith, 2014). We placed additional emphasis on the northwestern coast and the forested mountain ecoregions during evaluation. We then restricted our uncertainty analysis to these two ecoregions where more thorough model evaluation was performed.

Methods
We present a three-step methodology for constraining plausible parameter space and selecting DVM parameterizations that satisfy a range of evaluation criterion. In step 1, we identified parameters with influence on the modeled vegetation state. In step 2, we explored the sensitivity of the modeled vegetation state to the parameters selected in step 1, identifying regions of parameter space that improve the agreement of the simulated vegetation with observations. In step 3, we selected multiple model parameterizations that satisfied a range of evaluation criteria. We then ran an ensemble of early and mid-21st century simulations and compared the projected change in vegetation distributions among parameterizations. Table 1 describes the experimental design for each step and summarizes the methodology.

Step 1: Identification of Influential Parameters
DVMs simulate ecosystem processes related to plant growth, competition, and mortality, including photosynthesis, evapotranspiration, biochemical, and biophysical processes. Standard representations of these processes are commonly used across DVMs (Sitch et al., 2008). Many parameters are comparable among DVMs and have been evaluated by previous studies (Booth et al., 2012;Li et al., 2016;Pappas et al., 2013). We identified 38 parameters that either exhibited sensitivity in previous experiments or were of interest due to their role in modeled processes (Table A1). We set plausible ranges (Table A2) for each parameter based on literature review, observational constraints, and expert elicitation (personal communication). If no information was available, the parameter range was set to half and doubles the default setting. We acknowledge that the parameter range influences the resulting model sensitivity; thus, we error on the side of initially setting implausibly wide ranges and including parameters with little influence on model response. Throughout the methodology, the parameters and plausible ranges are systematically constrained and implausible parameter space is ruled out. One-at-a-time parameter perturbations were performed with each parameter set to seven values with equal geometrical spacing across the defined range. For PFT-specific parameters, perturbations were made across all PFTs preserving the plant trait relationships among PFTs. At each parameter setting, we performed four-year simulations under contemporary forcings (January 2004 to December 2007) with TRIFFID equilibrium rounds taking place after the second and the fourth year. Simulations were initiated with three unique ICs, resulting in a 798-member ensemble (38 parameters × 7 values per parameter × 3ICs).
We were interested in the sensitivity of three primary response variables describing the vegetation state and competition dynamics in TRIFFID: fractional coverage, LAI, and canopy height. In this study, we focused on woody ecosystems and analyzed results for BL trees, NL trees, and SHs in the continental U.S. west of 106°W. Individual parameter influence was quantified as the range of the response variable across parameter settings, after averaging over ICs. To assess the magnitude of influence of a parameter, we compared the range in the response variable to uncertainty due to ICs alone. To estimate the reference uncertainty attributable to ICs alone, we performed simulations with the default parameterization and nine unique ICs. The range in each response variable across ICs was set as a threshold; if the range in the response variable across parameter settings was larger than this threshold range, we concluded that the parameter was influential. Parameters influencing any of the three vegetation state response diagnostics of BL, NL, or SHs were retained for subsequent sensitivity analysis in step 2. All other parameters were held at the default value for the remainder of the experiment. Step 1 Identification of influential parameters a. Model simulations with one-at-a-time parameter adjustments (Table A1) 38 parameters 7 parameter settings 3 initial atmospheric conditions 4-year simulations 2 TRIFFID equilibrium rounds (1 round consists of 10,000 iterations) b. Determine influential parameters (Figures 2, S1, and S2) 18 parameters selected Step 2 Sensitivity analysis a. Model simulations with simultaneous parameter adjustments (Table A1) 18 parameters selected in step 1 250 unique parameterizations 5 initial atmospheric conditions 4-year simulations 2 TRIFFID equilibrium rounds b. Use statistical emulation to explore parameter influences and interactions (Figures 3, 4, and S3-S9) Train emulators to predict response variables Explore one-at-a-time sensitivity Explore parameter interactions c. Constrain plausible parameter space ( Figure 5) Determine influential parameters Define plausible parameter ranges Step 3 Parameter set selection a. Model simulations within plausible parameter space (Table 2) 10 parameters/ranges selected in step 2 140 unique parameterizations 5 initial atmospheric conditions 4 year simulations 2 TRIFFID equilibrium rounds b. Evaluate model performance ( Figure 6) Eliminate parameterizations that do not meet evaluation criteria c. Select model parameterizations ( Figure S11 and

Step 2: Parameter Sensitivity Analysis
In step 2 we explored interactions among parameters and identified regions of parameter space that improve model performance. DVMs have complex model structure with most processes being influenced by multiple parameters within the same equations. One-at-a-time sensitivity analyses can miss important interactions (Saltelli & Annoni, 2010) leading to over/under estimation of sensitivities. However, a comprehensive assessment of sensitivity would require simulations to be performed at all possible parameter combinations, requiring unrealistic computational resources. We combined PPE modeling with statistical emulation as a computationally cost-effective solution to thoroughly explore model parametric sensitivity and parameter interactions.
A PPE of simulations was generated with the 18 parameters selected in step 1 (Table A1: bold). To efficiently sample the 18-dimensional parameter space, a Latin hypercube sampling design (McKay et al., 1979) was used to generate an ensemble of unique parameterizations (see Table A2 for parameter ranges). We ran the model with each unique parameterization for 4 years from January 1904 to December 1907, with TRIFFID running in equilibrium mode after the second and the fourth year of the simulation. The design allows vegetation transitions in the first equilibrium round to feedback onto the climate prior to the final equilibrium round.
The PPE, consisting of 250 unique parameterizations averaged over five ICs per parameterization, was used in the sensitivity analysis. Sensitivity was evaluated as the change in the vegetation state from the reference state after two equilibrium time steps. Three response metrics were used to evaluate the change in the vegetation state: (1) The absolute change in grid-cell fractional coverage averaged over the domain; (2) the spatial correlation between the resulting fractional coverage and the initial fractional coverage; and (3) the absolute change in grid cell above ground biomass averaged over the domain. We evaluated the sensitivity of these response metrics for BL trees, NL trees, and SHs.
To explore sensitivity, we used a Gaussian process emulator as a proxy for the model. In an ideal sensitivity analysis we could evaluate the model response at all possible combinations of parameter settings. Since this is computationally infeasible, we ran a strategically designed subset of model parameterizations, ⃑ x, and record the model simulated response, ⃑ y, then use an emulator to estimate the function The emulator interpolates between the simulated parameterizations to estimate the mode response at all other possible combinations of parameters. For the emulator, we used a kriging function to produce a predictive distribution for the model output given any input parameterization. The kriging model is fit assuming a linear trend a priori and adjusting with a zero mean stationary Gaussian process with unknown hyperparameters estimated by maximum likelihood. We used the kriging function km from the DiceKriging package (Roustant et al., 2012) in R (R Core Team, 2016) with the default covariance kernel and no "nugget" term, which assumes that there is no uncertainty at the model simulated points used to condition the emulator.
An emulator was trained to predict each of the three response metrics of interest for each of the three woody PFTs (9 emulators in total) as a function of DVM parameterization. Emulator validation can be found in Figure S3 in the supporting information. Once trained and validated, we used the emulators to estimate the individual parameter sensitivity and parameter interactions. The individual sensitivity was explored by holding all but one parameter fixed at the default value and using the emulator to predict the response when the parameter of interest is set to uniform, incremental values across the perturbation range ( Figures S4-S6). The emulated response for individual parameter adjustments was used to identify regions of parameter space where model behavior is improved, and whether model behavior at the edges of parameter ranges without prior information is plausible. We defined a sensitivity index as the variance in the emulated response metric across individual parameter ranges in order to compare the sensitivity of different PFTs across response metrics. Parameters were culled between steps 2 and 3 based on qualitative examination of Figures S4-S6.
We did not set a quantitative threshold for the sensitivity index below which parameters would be discarded and above which parameters would be retained specifically because of the dependence of the sensitivity index on the parameter range.
To quantify the contribution of each parameter to the total variance, we performed a Fourier amplitude sensitivity test (FAST; Saltelli & Bolado, 1998). The FAST parses the contribution to the total variance into main effects and interaction terms for each parameter. The main effects are the ratio of an individual parameter's contribution to the total variance and the interaction terms represent the proportion of variance contributed by interactions with all other parameters. Though the FAST analysis does not parse the interactions between each pair of parameters, we used it to help identify the influential parameters to be perturbed in step 3. Noninfluential parameters were held at the default value for the remainder of the experiment.
The default TRIFFID parameterization simulates the spatial distribution of forest coverage well when compared to an observational data set derived from satellite and FIA data (spatial correlation of 0.82±0.03); however, the fraction of forest cover is low across the entire domain (−45.9 ±3.8%). To identify regions of parameter space that the simulated fractional forest cover more accurately than the default parameterization, we calculated the grid cell absolute difference between the observation-based forest cover (Ruefenacht et al., 2008) and model-simulated fractional forest cover of each parameterization and averaged over forested ecoregions (northwestern coast and forested mountains level 1 ecoregions). We trained an emulator to predict the absolute error in fractional forest cover from model parameterizations. The Figure 1. Pairwise settings of parameters influential on the fractional coverage of needleleaf trees refined throughout the methodology. Parameter values in the step 2 perturbed parameter ensemble (PPE; light grey) exemplify the initial plausible parameter space. In step 3 the parameter space was constrained to improve model performance (dark grey). A final set of 10 parameterizations (black) were selected as improved model variants (Table S3). The default setting is shown in orange.
emulator was used to estimate the absolute error (%) at pairwise combinations of parameter settings while randomizing all other parameter settings. From the emulated response we identified interactions between pairs of parameters and regions of parameter space that reduced the absolute error. Results guided the selection parameter ranges for step 3.
In step 1 and 2 parameter adjustments were applied to all PFTs uniformly, maintaining the plant trait relationship between PFTs. These plant trait relationships were originally defined for global-scale applications and may not be appropriate for our region of interest. In step 3 we varied parameters independently for individual PFTs, allowing for the relationship between PFTs to vary.

Step 3: Parameter Set Selection
The objective of step 3 was to select parameterizations of the DVM that meet evaluation criteria to capture a range of plausible model behaviors. The parameter space determined in step 2 was resampled with a Latin hypercube design to generate 140 unique parameterizations. Four-year simulations were performed with the same protocol as in step 2, and results were averaged over five ICs for each parameterization.
Changes to the DVM were made through adjusted parameter settings; no structural changes were made. We therefore use the default model parameterization as a benchmark and judge alternative model parameterizations based on whether they improve agreement between the simulated vegetation state and observations. For reference we ran a 40-member IC ensemble with the default model parameterizations referred to hereafter as the default ensemble. Four-year simulations were performed for December 1903 through November 1907 with TRIFFID running in equilibrium mode after the second and the fourth year. Using a coupled landatmosphere model allows the vegetation transitions after the first equilibrium round to feed back onto the atmosphere during the third and fourth year, which, in turn, impact the vegetation transitions in the second equilibrium round. Simulations were initialized from the preindustrial spinup described in section 2.2. The initial vegetation state was identical for all simulations. equilibrium time steps is shown for seven unique settings per parameter. Each point represents a different parameter setting averaged over three initial conditions. Parameters are grouped by the modeled process they influence (A-G). A 9-member initial condition ensemble with default parameter settings (black open circles) was used to estimate the spread in model results due to initial conditions alone. Adjustments to parameters shown in color introduce uncertainty resulting in a larger spread than initial conditions in fractional coverage, leaf area index ( Figure S1), or canopy height ( Figure S2). Parameters shown in color for needleleaf (green), broadleaf (purple), and shrubs (brown) were selected for perturbation in step 2.
We evaluated the simulated dominant vegetation type at each grid cell over the domain and calculated the percent agreement with observation-based potential dominant vegetation type (described in section 2.3). The default ensemble was used to establish an evaluation threshold: the gridcell agreement of potential dominant PFT (NL, BL, or SH/grass) between the default ensemble and the observational data set ranged from 72.7 to 74.5% of grid cells in the domain; the median percent agreement (73.6%) was set as a threshold. Parameterizations that did not improve the agreement with observational data were discarded. The remaining parameterizations were subjected to further evaluation.
The dominant PFT distribution metric evaluates how the model simulates the relative proportions of PFTs in specific locations. We further evaluated the model-simulated magnitude of the fractional forest coverage against an observational data set derived from satellite and FIA data. The biomass and canopy height of forested ecoregions were evaluated against observational data from Wilson et al. (2013) and Simard et al. (2011), respectively. Areas of large historical land use and land cover changes were masked out since observational data sets are for the contemporary period and simulations are performed without socioeconomic activities and under preindustrial atmospheric conditions. The observational data sets and land use land cover mask are described in section 2.3. We selected 10 model parameterizations that simulated all three evaluation criteria well (forest coverage, biomass, and canopy height). We consider the selected 10 parameterizations to be improved model variants of TRIFFID for the western United States.
To visually depict the methodology, the evolution of model parameter space after step 1 is shown in Figure 1 (step 1 eliminated parameters but did not otherwise constrain parameter ranges). Settings of parameters most influential on NL forests (see Figure S7 for additional parameters) are shown in pairwise parameter space with the default setting highlighted in orange. Parameter settings in the 250-member PPE in step 2 are shown in light grey. Statistical emulators were built to interpolate between simulated points, then the emulators were used to identify influential parameters, explore parameter interactions, and identify regions of parameter space that improved model performance. The newly defined parameter space was resampled to select 140 parameter settings for the PPE of simulations in step 3 (dark grey). Model simulations were evaluated against observational data and 10 parameterizations that performed well in evaluation criteria were selected (black). The 10 model parameterizations were considered unique model variants and used to assess uncertainty in model response to climate forcing.

Uncertainty in Response to Climate Forcing
Multiple model parameterizations yielded similar equilibrium vegetation distributions that agreed well with observations. Model variants with alternate parameterizations may respond differently to changes in forcing (e.g., climate conditions, temperature, water limitations). To assess the uncertainty in the forced response of alternate model parameterizations, we evaluated the projected change in equilibrium vegetation state under a preindustrial climate compared to both the contemporary climate and a projected mid-21st century climate under Representative Concentration Pathway 8.5 (Meinshausen et al., 2011).
Preindustrial simulations were initiated from the equilibrium vegetation state obtained for each of the 10 parameterizations selected in step 3. Additional 6-year, 8-month simulations (January 1900 to August 1906) with two equilibrium rounds occurring after the third and sixth year were performed to ensure that the modeled vegetation state had reached an equilibrium distribution under preindustrial atmospheric conditions. For each parameterization, the vegetation state at the end of the additional spinup period was used to initialize 5-year simulations for 1 September 2004 to 30 September 2009 with one TRIFFID equilibrium round occurring at the end of the simulation. The same procedure was followed for the period 1 September 2054 to 30 September 2059. Atmospheric forcings for the preindustrial, historical, and future time periods (described in section 2.2) applied anthropogenic radiative forcings of 0.2, 1.9, and 5.1 W/m 2 , respectively. The resulting change in global mean temperature was 1.42°C between the historical period and the preindustrial period and 3.59°C between the future period and the preindustrial period. We introduced uncertainty in the modeled climatic conditions over the 5-year periods by using a 30-member IC ensemble for each parameterization, resulting in two 300-member ensembles of 5-year simulations. We hereby refer to these two ensembles as the historical and future experiments, respectively.
We focused our assessment of uncertainty due to model parameterization on the conversion of NL forest in response to climate forcings. We restricted our analysis of changes in the forest fraction within the northwestern coast and forested mountain ecoregions since model evaluation of fractional coverage, biomass, and canopy height was performed only in these ecoregions.

Step 1: Identification of Influential Parameters
An assessment of the sensitivity of a response to parameters in unavoidably influenced by the chosen parameter perturbation ranges. Achieving measures of sensitivity standardized across parameters is challenging, if even possible. We acknowledge, therefore, that the intercomparison of parameter sensitivities is dependent on the experimental design.
Parameters with the largest influence on vegetation response were related to the allocation of carbon to growth and expansion (Figure 2, group A). The sensitivity is driven by the strong dependence of modeled PFT competition on the representation of plant structure (i.e. height, biomass, and LAI) and thus carbon allocation. This result highlights the importance of improved representation of plant competition for light, water, and nutrients among PFTs.
The carbon allocation parameters A_WL, A_WS, B_WL, and ETA_SL relate biomass to canopy height, h, as where W is the total stem carbon (Cox 2001). A strong response to perturbations of all four parameters was observed in both canopy height ( Figure S2, group A) and fractional coverage (Figure 2, group A). B_WL had a disproportionate response because it is in the exponent (1) and, without any prior information on parameter sensitivity, the perturbation range was set to half and double the default value. Because of the functional form of allocation (1), only A_WL and B_WL are required to adjust the shape and scale of carbon allocation functions, and thus, A_WS and ETA_SL were held at their default value for the remainder of the experiment.
Competition among PFTs in TRIFFID is primarily driven by LAI and canopy height. Once a minimum LAI (LAI_MIN) is attained, plants begin to compete with surrounding vegetation for light and space. When a PFT's actual LAI is near the minimum LAI (LAI_MIN), more carbon is allocated to growth, but as the actual LAI approaches the maximum LAI (LAI_MAX), more carbon is allocated to expanding the areal fractional coverage of the PFT. NL trees, BL trees, and SHs were all sensitive to perturbations in LAI_MIN and LAI_MAX in all three response metrics (Figures 2, S1, and S2, group A), highlighting the importance of carbon allocation and competition to the vegetation state and distribution.
Fractional coverage and LAI were sensitive to some parameters related to growth and turnover (Figures 2 and S1, group B). BL trees were particularly sensitive to the leaf growth rate (G_GROW) and the temperature threshold for leaf drop (TLEAF). Litterfall is a linear function of the leaf, root, and stem wood turnover rates in TRIFFID. Stem wood turnover rates are not well constrained in the literature, whereas leaf and root turnover rates are better known. For this reason, only the rate of stem woody biomass turnover (G_WOOD) was selected for additional analysis in step 2.
Mortality is simply represented by a background disturbance rate (G_AREA) that varies by PFT; episodic disturbances such as fire, windthrow, or insect attack are not explicitly represented. The default values for each PFT were chosen to yield realistic lifecycles (Cox 2001), but there is evidence that these rates may be increasing under anthropogenic climate change (van Mantgem et al., 2009). As expected, the fractional coverage of all three PFTs was sensitive to changes in G_AREA.
Parameters related to photosynthesis, respiration, and transpiration influence NPP, which is then allocated into growth and expansion of PFTs. The minimum and maximum temperatures at which plants can photosynthesize (TLOW and TUPP, respectively) vary by species. While well constrained in the literature for specific species, these parameters are aggregate representations of all species of, for example, NL trees. Wide ranges were set for both TLOW and TUPP, and sensitivity was observed in both LAI ( Figure S1) and fractional coverage (Figure 2). There is a high degree of uncertainty in observational estimates of the growth respiration fraction (R_GROW; the energy cost for constructing organic compounds fixed by photosynthesis). The MODIS NPP satellite data product assumes the growth respiration fraction to be 25% of NPP (Zhao et al., 2005). Land surface models used by the Land Data Assimilation System (LDAS) also assume a growth respiration fraction of 25%. Here, R_GROW was varied between 12% and 50% resulting in detectable changes in LAI and fractional coverage. Constraining the plausible range of R_GROW, and improving the model representation of the temperature dependence of respiration would reduce modeled uncertainties. Notably, we did not detect model sensitivity to variations in quantum efficiency (ALPHA); however, Li et al. (2016) and Pappas et al. (2013) found perturbations in quantum efficiency to strongly influence gross primary productivity. Stomatal conductance is calculated via a simplified Leuning model (Jacobs et al., 1994) where relative humidity in the Ball-Woodrow-Berry model has been replaced with a dependence on the critical humidity deficit (DQ_CRIT; Cox et al., 1998). We did not detect sensitivity to the minimum leaf conductance (GLMIN) and CO 2 compensation point (F0) in fractional coverage, LAI, or canopy height (Figures 1, S1, and S2 respectively). The absence of model sensitivity to parameters related to the influence of atmospheric CO 2 concentrations on photosynthesis is partially due to the experimental design since simulations are performed under near constant atmospheric CO 2 concentrations. Changes to the critical humidity deficit (DQ_CRIT) moderately influenced the fractional coverage of SHs ( Figure 2, group F). Photosynthesis becomes limited when the atmospheric demand for moisture is high and the critical humidity deficit is crossed. This only occurs in specific seasons and for short intervals. The response metrics used in this study are a function of the aggregated NPP; thus, model sensitivity to the DQ_CRIT parameter was low. To better evaluate the influence of DQ_CRIT, a sensitivity experiment using response metrics that capture extreme conditions would be required.
Parameters related to radiative balance and albedo (Figure 2, group D) exhibited little to no effect on the vegetation state and were thus held at the default value for the remainder of the experiment. The limited influence of radiation parameters on model response highlights the poor representation of competition for light within the canopy. The large response to photosynthetically active radiation extinction coefficient (KPAR) perturbations was due to the unrealistically large range over which this parameter was perturbed. The KPAR parameter has been constrained in the literature (Chen et al., 1997), and little sensitivity was observed when constrained to a realistic range.
LAI and fractional coverage were sensitive to parameters related to nutrient composition (group E: SIGL, NL0, and NR_NL). The parameter relating stem nitrogen to leaf nitrogen (NS_NL) showed little influence and was left at its default value in subsequent simulations. Parameters related to hydrology had little effect on the modeled vegetation state, except for rooting depth (ROOTD), which controls the soil moisture accessible to the plant. ROOTD has a large influence on the fractional coverage of SHs, but little influence on tree PFTs. The lack of model sensitivity to large adjustments in the rooting depth parameter for tree species indicates that the model may not adequately represent water limitations on growth.  Clark et al., 2011) using gross primary productivity and latent heat of evaporation as response metrics and found similar parameters to be most influential. Additionally, Li et al. (2016) found that parameters related to the canopy catchment capacity (DCATCH_DLAI) and the change in roughness length with height (DZ0V_DH) strongly influence modeled gross primary productivity and latent heat of evaporation. We chose to vary these parameters in subsequent steps because of the importance of gross primary productivity and latent heat of evaporation to land-atmosphere feedback that are inherently captured in our study design but not explicitly evaluated in step 1.

Step 2: Parameter Sensitivity Analysis
In step 2 we constrained the parameter space by identifying the most influential parameters and processes, estimating interactions among parameters, and identifying regions of parameter space that improve the simulated fractional forest cover. We describe how the sensitivity analysis informed the selection of parameters and perturbation ranges used in step 3.
We evaluated the change in vegetation state after equilibrium simulations with three response metrics: the regional mean of the absolute grid-cell change in fractional coverage, the spatial correlation of fractional coverage, and the regional mean absolute grid-cell change in biomass. The emulated one-at-a-time sensitivity of each response variable to individual parameters is shown in Figures S4-S6 and summarized by the sensitivity index (the variance of the emulated one-at-a-time sensitivity; section 3.2) in Figure 3 for NL trees, BL (BL) trees, and SHs. The sensitivity index was standardized for each response metric (i.e., rowwise standardization) to identify the parameters most influential on each response variable. The sensitivity index is unavoidably dependent on the perturbation range set for each parameter. Parameters without observational constraints were set to one half and double the default setting, and thus, the sensitivity index for parameters not constrained by observations should not be directly compared to parameters with observationally constrained ranges. For each parameter, the sensitivity index shown in Figure 3 can be compared across response metrics.
The vegetation state variables of all three PFTs were sensitive to parameters related to carbon allocation and competition ( Figure 3, group A). Nearly all response metrics were sensitive to A_WL and B_WL. B_WL was the single most influential parameter on the biomass of all three PFTs. In step 3, we varied B_WL uniformly for all PFTs, as in step 2, but allowed the relationship of A_WL between PFTs to vary by changing A_WL of each PFT independently of one another. We took the same approach for other parameters that exhibit influence on specific PFTs. LAI_MIN influenced the vegetation state of all three PFTs but only SHs showed sensitivity to LAI_MAX. In step 3 we varied LAI_MIN independently on all three PFTs, but LAI_MAX was only varied for SHs.
Different growth and turnover parameters (Figure 3, group B) showed influence on specific PFTs. NL trees and SHs were sensitive to the turnover rate for woody biomass (G_WOOD), and BL trees were highly sensitive to the temperature threshold for leaf drop (TLEAF). In step 3 growth and turnover parameters were varied only for sensitive PFTs and the range for TLEAF was constrained according to Woodward & Williams (1987; Table 2).
The fractional coverage and spatial correlation of all PFTs were sensitive to the background mortality rate (G_AREA). Lower disturbance rates reduced loss of BL and NL trees, whereas higher background disturbance rates prevented SHs from expanding into grassland ecosystems. In step 3 we set the background disturbance rate to 0.3% and 0.25% per year for BL and NL trees, respectively (Table 2) but held G_AREA at the default setting for SHs and grasses.
The minimum and maximum temperatures for photosynthesis (TLOW and TUPP) are fairly well constrained in the literature (Collatz et al., 1991) with thermal inhibition of C3 plants occurring above 35°C (Berry & Raison, 1981). Increasing TUPP increased the fractional coverage of NL forest, but HadAM3P/RM3P is known to have a warm summer bias in inland regions even with the bias-reducing alternative climate model parameterization we used from Li et al. (2019). We are cautious not to overly tune the vegetation model to a biased climate so we left TLOW and TUPP fixed at the default value but allowed TUPP to vary between 31 and 36°C for NL trees.
The vegetation equilibrium state was moderately sensitive to nutrient-related parameters. The top leaf nitrogen concentration (NL0) and the specific density of leaf carbon (SIGL), together define the nitrogen content per one-sided leaf area. The default values were based on Schulze et al. (1994), who suggests that the nitrogen content in BL trees and grasses should be 0.00157 gN/m 2 . The nitrogen content in NL trees is higher than BL trees and grasses. Observations (Schulze et al., 1994) suggest that NL0xSIGL should span a range of 0.00167 to 0.0059 gN/m 2 for NL trees. Although issues related to using one-sided leaf area for NL trees introduces some uncertainty, we held SIGL fixed at the default value and varied NL0 such that NL0xSIGL spans the observational range. The MOSES2-TRIFFID model does not explicitly represent nutrient limitation on growth. Improving the model representation of nutrient limitations on photosynthesis would introduce physiological constraints, but also additional parametric uncertainty that must be adequately quantified and accounted for in projections.
The rooting depth (ROOTD) had little influence on NL and BL trees since roots are always able to access soil moisture at deeper layers in the soil column. While transpiring vegetation extracts water from all soil layers as a function of the root depth and density, the canopy conductance becomes limited by soil moisture availability when the water content in a soil layer lies below a critically limiting soil water content. The critically limiting soil water content strongly controls plant response to moisture stress; however, it has little physical basis (Kala et al., 2016;Powell et al., 2013); thus, it is challenging to calibrate against observational data and by plant species. The critical point for BL and NL trees is set to a low value (i.e., trees are free to photosynthesize even at low soil moisture levels), which affects model sensitivity to adjustments of the rooting depth parameter. SHs have a shallower rooting depth and limited access to soil moisture. Increasing the ROOTD allowed SHs to access water in deeper soil layers and better compete with trees, resulting in notable changes in vegetation composition in dry, forest-dominated ecosystems. In step 3 we varied the rooting depth for only SHs to change the relative rooting depths between PFTs and influence competition.
Several parameters showed little to no influence on any of the response metrics and were held at the default value for the remainder of the experiment (DCATCH_DLAI, DQCRIT, DZ0_DH, NR_NL, R_GROW, and G_GROW). The low sensitivity of model response to adjustments in these parameters may be due to the response variables chosen in this study; these parameters may be influential on other model processes.
Additionally, the low sensitivity may be attributable to a model structural deficiency, missing processes, or poorly represented processes.
We quantify parameter interactions with the FAST analysis and show results for the fractional coverage of NL, BL, and SHs in Figure 4 (results for spatial correlation and biomass are shown in Figures S8 and S9, respectively). The proportion of the overall variance contributed by individual parameters is shown in darker shades and interactions with all other parameters are shown in lighter shades. Parameters related to carbon allocation and competition (group A) have the most influence on fractional coverage of all PFTs, and the interactions among parameters contribute substantially to the overall variance. B_WL is highly influential on the fractional coverage of SHs (brown), and over 50% of the contribution to the overall variance is attributable to interactions with other parameters. For BL trees (purple), the temperature below which trees drop their leaves (TLEAF) is the most influential parameter on the fractional coverage since this parameter controls the climate conditions in which BL trees can establish, grow, and compete. While B_WL, LAI_MAX, and DQCRIT all contribute at least 10% of the overall variance in BL tree fractional coverage, the majority of the contribution to the variance is due to interactions with other parameters. NL tree (blue) fractional coverage is sensitive to several processes and parameters: carbon allocation and competition (A_WL, B_WL, LAI_MIN), turnover (G_AREA), and nutrient composition (SIGL, NL0). For NL trees, interactions among parameters were smaller than the main effects of each parameter.
Parameters were selected for perturbation in step 3 based on the emulated one-at-a-time sensitivity analysis, which identified the most influential parameters, and the FAST, which parsed the main effects from the interactions. To identify regions of parameter space that improved the simulated fractional forest cover, we trained an emulator to predict the absolute grid-cell error in fractional forest cover. We used the emulator to estimate the error at all combinations of pairs of selected parameters, while randomizing the settings of all other parameters. Results shown in Figure 5 identify regions of parameter space with the lowest error in  fractional forest cover. Setting A_WL, B_WL, and LAI_MIN to the lowest end of the respective perturbation ranges results in the lowest error in simulated fractional forest coverage. A_WL, B_WL, and LAI_MIN are parameters that arise from the functional form of the model and cannot be directly observed or estimated. We selected parameter ranges for step 3 ( Table 2) by identifying that the lowest values of A_WL, B_WL, and LAI_MIN result in the least error, and inferring that the error may be further reduced by lowering the perturbation range. Without prior information on these parameters the lower bound was set as one half the default value in steps 1 and 2. In step 3 we further reduced the lower bound of the perturbation range with the expectation of further reducing the error in simulated fractional forest cover. We validated the bounds set for the allocation parameters A_WL and B_WL by comparing the modeled biomass and canopy height (equation (1)) against forest inventory analysis data (Berner & Law, 2016) to verify that the selected ranges simulate a realistic relationship between biomass and canopy height ( Figure S10).
Lower settings for G_WOOD reduce the error in fractional forest cover ( Figure 5) so we constrained the plausible range for step 3 between the lower bound in step 2 and the default value. G_WOOD was varied on NL trees and SHs independently. TLEAF shows little influence on the emulated error in forest coverage ( Figure 5) although it contributes largely to the variance of BL trees (Figure 4). The western United States is heavily dominated by NL forest, and since TLEAF is a BL specific parameter, it does not have a notable influence on the overall forest cover. LAI_MAX and ROOTD are excluded from Figure 5 since they were only perturbed for SHs in step 3 (Table 2).
Step 2 further eliminated several parameters that exhibited little influence on the equilibrium vegetation state. The range of the remaining parameters was chosen to increase the simulated fractional forest cover, reducing the error when compared to observations. The selected parameters and ranges used in step 3 are shown in Table 2. Perturbations to B_WL were applied to all PFTs uniformly, while for all other parameters, PFTs were treated as independent from one another resulting in 14 free parameters.

Step 3: Parameter Set Selection
In step 3 we identified parameterizations of the DVM that meet evaluation criteria. The vegetation distribution after two equilibrium time steps of 140 parameterizations selected from the refined parameter space in step 2 was evaluated against the potential natural vegetation distribution from Küchler 1965. Parameterizations not meeting the dominant PFT evaluation criteria (section 3.3) were discarded, resulting in 53% of parameterizations being retained for further evaluation. The fractional forest coverage (NL + BL) of the retained parameter sets was evaluated against observations within forested ecoregions including the northwestern coast and forested mountains (section 2.3; Figure 6 map inset). The retained parameterizations reduced the error in forest fractional coverage when compared to the default ensemble ( Figure 6). Parameterizations with the least error in fractional forest cover (points in the lower right corner of Figure 6) obtained the increase in fractional cover by reducing the canopy height and biomass of the forest. This relationship is attributable to the functional form of the model: in TRIFFID, carbon (NPP) is either allocated to growth (biomass) or expansion (fractional coverage). Thus, parameterizations that allocate more carbon to expansion improve the simulated fractional forest coverage, but the resulting forest is less dense with a lower canopy height. This relationship is depicted in Figure 6 for the PPE (filled circles) and the default ensemble (open circles). Parameterizations simulating fractional forest coverage near observed coverage (diamond) have simulated biomass below half of the observed biomass. Parameterizations with biomass near that observed simulate low fractional forest cover. Parameterizations and observations in Figure 6 are shaded by canopy height, and the default ensemble is outlined by canopy height (m). The canopy height generally covaries with the biomass, but some parameters shift the relationship between biomass and canopy height. There are several parameterizations that yield unrealistically tall forests with high biomass but very low fractional coverage (points in the upper left corner of Figure 6).
While the functional form of the model governs the shape of the relationship among fractional coverage, biomass, and canopy height, the magnitude of the fractional coverage and biomass is directly proportional to the total amount of carbon available (NPP) for allocation to growth and expansion. NPP simulated by HadRM3p-MOSES2 is generally biased low when compared to MODIS satellite observations (not shown), which limited improvements in DVM performance through parameter adjustments. Results in Figure 6 are sized by the 4year average monthly NPP for each simulation; each IC member is shown for the default ensemble (40 ICs), but results from the PPE are averaged over the 5-member IC ensembles for each parameterization. In the default ensemble, increasing the NPP increases both the biomass and the fractional coverage. This relationship is slightly confounded in the PPE since several parameters (NL0, ROOTD, TUPP) influence NPP.
A subset of 10 parameterizations (Table S3) were selected that improve the simulated fractional forest cover when compared to the default ensemble, but not at the expense of the biomass. On average the 10 selected parameterizations reduced the error in fractional coverage by 31±9%, while maintaining the above ground biomass to within ±30% of observations. The selected parameterizations have dominant vegetation distributions that agree well with potential natural dominant vegetation distributions (74.4% to 77.3% agreement; Figure S11) and show little to no drift in the equilibrium vegetation state during additional preindustrial spinup simulations (not shown). The 10 selected model parameterizations perform similarly in evaluation metrics, and yet have dissimilar parameter settings (Figures 1 and S11 and Table S3). This exemplifies the importance of parameter compensations and interactions, which our methodology for selecting model parameterizations implicitly takes into account. The evolution of the plausible parameter space of parameters adjusted in the final model parameterizations are depicted in Figure 1 (NL parameters) and S7 (all other parameters). The 10 Observation based estimates of fractional forest coverage (Ruefenacht et al., 2008), biomass (Wilson et al., 2013), and canopy height (Simard et al., 2011) are shown (diamond) after masking land-use land-cover change (Meiyappan & Jain, 2012).

10.1029/2018MS001577
Journal of Advances in Modeling Earth Systems selected parameterizations (black) show a notable spread in parameter settings but do not deviate greatly from the default setting (orange). Each parameterization is considered a unique model variant. We quantified the uncertainty among parameterizations to assess the influence of parameter compensations and interactions on model response to climate forcing.

Uncertainty in Response to Climate Forcing
We evaluated uncertainty in the response to climate forcing due to model parameterization using the 10 parameterizations selected in step 3. We assessed the change in equilibrium vegetation state under the contemporary climate (historical) and a projected mid-21st century climate (future) compared to a preindustrial climate (section 3.4). The change in NL forest fraction in the future ensemble relative to the preindustrial equilibrium vegetation state is shown for each parameterization (Figure 7) averaged over 30 IC ensemble members. We restrict our analysis to the northwestern coast and the forested mountains ecoregions since our evaluation and selection of model parameterizations was performed with additional observational data in these regions. Parameterizations generally agree on a net increase of NL fraction across the inland mountains and a decline in NL fraction in the northwestern coast. However, the spatial distribution, magnitude, and in some cases, direction of changes vary widely among parameterizations (Figure 7).
The ecoregion mean change in the historical and future ensembles are shown in Figure 8 for NL versus BL (a and c) and NL versus SH (b and d). All parameterizations agree on an increase in both NL and BL trees in the forested mountain ecoregion (a) by the future period. The IC spread within an individual parameterization is on the order of tenths of a percent change in forest cover, whereas the spread across parameterizations is around 10% change in forest cover. Much of the forested mountain ecoregion lies at high elevations where growth is temperature limited during parts of the year. The growth and expansion of trees is due to an increase in NPP attributable to a temperature increase, earlier spring growing conditions, and CO 2 fertilization. In the historical ensemble, the expansion of NL trees is associated with a loss of SHs in most parameterizations (Figures 8 and 8b). By the future period growing conditions are suitable for SHs to grow and expand in some areas of the forested mountain ecoregion as well. Averaging over large ecoregions masks finer spatial transitions, but this study is not intended for projecting local vegetation transitions.
Here we simply illustrate the uncertainty in large-scale patterns of model response to climate forcing due to parameterization.
In the northwestern coastal ecoregion (Figures 8c and 8d) most parameterizations project a transition of NL to BL trees. However, the magnitude of the transition varies from a 4% increase in both NL and BL to a 28% conversion of NL to BL across parameterizations in the future ensemble. Very little change in SH fraction is observed in most parameterizations, with the exception of parameterization 6, which showed an increase in both NL and BL trees and a corresponding decline in SH fraction.
The uncertainty in vegetation state response to climate forcing among parameterizations is large in comparison to uncertainty due to ICs. Figure 8 shows the ensemble mean estimate ±3 standard errors estimated from a 30-member IC ensemble for each parameterization. For most parameterizations, ±3 standard errors are on the order of tenths of a percent change. The uncertainty due to parameterizations is much larger than the Figure 8. Change in fractional coverage of (a and c) needleleaf trees and broadleaf trees and (b and d) needleleaf trees and shrubs in the historical ensemble (hist; dark) and future ensemble (future; light) from the preindustrial equilibrium vegetation state (preind) for the 10 parameterizations selected in step 3. Fractional change is averaged over the forested mountain (a and b) and the northwestern coast (c and d) ecoregions (map insets). Parameterizations are numbered as in Figure 7 with the historical and future ensemble mean connected with grey lines. The mean change ±3 standard errors (estimated from 30 initial condition ensemble members per parameterization) are shown for each parameterization.
uncertainty due to ICs, illustrating that model parameterization is an important source of uncertainty in vegetation transitions and carbon budgets.
We tested for statistical differences in the vegetation state simulated by different model parameterizations under preindustrial, historic, and future climate conditions. We calculated the dissimilarity between the simulated NL fractional coverage of 10 different model parameterizations under preindustrial conditions using the Euclidean distance. We compared the dissimilarity matrix with the historical and with the future dissimilarity matrices using a Mantel test, which calculates the Pearson correlation between two dissimilarity matrices, then permutes the rows and columns of the matrices to evaluate significance. We found that the dissimilarity in simulated NL fractional coverage under preindustrial conditions among the 10 model parameterizations was significantly different from the dissimilarity in NL fractional coverage simulated under historical conditions (p<0.01). The dissimilarity in NL fractional coverage among the 10 parameterizations was significantly different between the preindustrial and future periods (p<0.01) and between the historical and future periods (p<0.01). We repeated this analysis for the fractional coverage of BL trees and SHs and found that the dissimilarity matrices for each time period were significantly different from one another for both PFTs (p<0.01). This analysis was performed with the dist function in R (R Core Team, 2016) and the mantel function from the vegan package (Oksanen et al., 2011).

Discussion and Conclusions
DVMs are often parameterized to match observed vegetation distributions; however, this study demonstrates that many model parameterizations are plausible given available benchmark data. The resulting uncertainty in parameter values could contribute an important amount of uncertainty in vegetation transitions and carbon cycle feedback under nonstationary climate conditions. The impacts of parametric uncertainties may be amplified or reduced by land-atmosphere feedback. To quantify parameter sensitivity while accounting for land-atmosphere feedback, we performed coupled model simulations. Our results highlight the importance of considering model sensitivity to parameterization when choosing model variants, which can contribute to considerable uncertainty in forced model response.

Process Sensitivity
Parameters related to carbon allocation and competition had the most influence on modeled equilibrium vegetation distribution. Competition in TRIFFID, and some other dynamic vegetation models (Gotangco Castillo et al., 2012;Arora & Boer, 2006), is simulated based on the Lotka-Volterra representation of competitive ecological processes. A dominance hierarchy (tree-SH-grass) determines the outcome of competition between PFTs. This model structure is designed to represent natural competition for light and is thus a good representation of competition between trees and grasses, for example. However, this assumption breaks down when stand composition is complex or when BL and NL trees compete for light. The competition coefficients for within species competition are based solely on canopy height; competition for water, light and nutrients are not explicitly represented. A more complete understanding of Note. Parameters in bold were identified as influential in step 1 and included in the parameter sensitivity analysis in step 2. interspecies and intraspecies competition for light and nutrients is required before these processes can be incorporated into DVMs. Examining the model structure illustrates why the parameters with the highest sensitivity are A_WL and B_WL, which relate stem biomass to canopy height. Neither of these parameters can be measured through observational studies, they are simply a coefficient and an exponent in the functional relationship developed for this model (Cox 2001). However, the simulated relationship between biomass and height can be evaluated against observations ( Figure S10). Currently, all models use empirical allocation schemes (Franklin et al., 2012) and improved understanding is needed to mechanistically represent carbon allocation in DVMs such that parameters could be constrained through observations. Some limiting aspects to modeling natural systems are inevitable but the uncertainty introduced can be quantified and accounted for through thorough sensitivity analyses.
TRIFFID was also highly sensitive to parameters that allocate carbon (NPP) into growth versus expansion. Plants begin to compete for space once the LAI surpasses a minimum LAI required for competition and allocate more carbon into expansion as they approach the maximum LAI. The functional form places a high degree of importance on these two parameters, which have only abstract physical meanings. The functional form constrains the relationship between fractional coverage and biomass limiting model improvements through parameter settings ( Figure 6). This methodology for performing a sensitivity analysis not only assessed the model parameter sensitivity but also allowed us to identify which parameterized processes are most influential on model performance. We highlight the dependence of the simulated vegetation state on the structural representation of growth and competition in the TRIFFID DVM.
Several parameters exhibited very little influence on the modeled vegetation distribution including the rooting depth (ROOTD), the maximum ratio of internal to external partial pressure of CO 2 (F0), and radiation parameters. Low sensitivity may be attributable to a physical insensitivity; trees that are not water limited would not exhibit sensitivity to adjustment of ROOTD. The choice of response variables also impact the detected sensitivity; for example, had this experiment set out to quantify the parametric sensitivity of NPP to a doubling of atmospheric CO 2 concentrations then the F0 parameter would likely have had a larger influence on the results. Additionally, low model sensitivity may be a result of incomplete process representation. For example, competition among tree species is partially determined by radiation dispersion through the canopy, but the radiative transfer model used in TRIFFID assumes homogenous diffusion of light throughout the canopy, which does not adequately capture the effects of canopy structure on competition. The low sensitivity of modeled plant competition to the adjustment of radiation parameters may be attributable to poor process representation. Identifying parameters with low sensitivity can, in some cases, reveal model structural limitations.

Constraining Model Parameterizations
DVM parameterization contributes a significant amount of uncertainty to equilibrium vegetation distributions and vegetation response to climate change. Combining dynamically downscaled PPEs with statistical emulation techniques is a computationally cost effective method for understanding parameter and process sensitivity and selecting model variants. The approach we took for parameterizing an "off-the-shelf" model for a specific domain does not optimize individual parameter settings but rather allows for model parameter interactions to emerge and perhaps compensate slightly for one another. Our approach constrains

10.1029/2018MS001577
Journal of Advances in Modeling Earth Systems the plausible parameter space in each step without optimizing performance for any specific diagnostic.
Here we select 10 model parameterizations that perform similarly well in evaluation metrics, and yet have very different sensitivities to nonstationary climate conditions. We illustrate how parameter interactions can compensate for one another such that many different parameterizations yield qualitatively similar equilibrium dominant vegetation distributions ( Figure S11) but respond differently to forcings. Accurately simulating the vegetation distribution over a large domain with a range of climatic zones does not imply a similar transient response from a DVM; that is, space-for-time is not adequate.
While the functional form of the TRIFFID model governs the shape of the relationship among fractional coverage, biomass, and canopy height, the magnitude of the fractional coverage and biomass is directly proportional to the total amount of carbon available (NPP) for allocation to growth and expansion (Cox 2001). By performing coupled land-atmosphere simulations we account for the feedback between vegetation transformations and atmospheric conditions capturing the resulting impacts on growth and NPP.

Uncertainty in Response to Climate Forcing
Results from all 10 parameterizations show an increase in NL forest cover in the forested mountains and a transition of NL forest to BL forest in the northwestern coast under future climate conditions ( Figure 8). The expansion in the forested mountains is attributable to increased NPP under a warmer climate since growth is primarily limited by temperature at high elevations. The transition of NL forest to BL forest along the northwestern coast varies widely among parameterizations ranging from an~4% increase in both NL and BL to añ 28% conversion of NL to BL. The uncertainty due to parameterization is large enough to have implications on carbon stocks, ecosystem services, and potential management decisions. Parametric uncertainty compounds on model structural uncertainties, which have both been shown to profoundly influence the ecosystem response to changing climate conditions (Cramer et al., 2001;Fisher et al., 2010;Sitch et al., 2008). The terrestrial carbon cycle is the least constrained component of the global carbon cycle (Le Quéré et al., 2018) and improved process representation and parameterization of carbon assimilation, growth, and carbon stocks in earth system models has the potential to further constrain the response to increasing atmospheric CO 2 (Fatichi et al., 2019).
The conversion of NL forest to BL forest in the coastal ecoregion between the preindustrial and historical periods is proportional to conversions between the historical and future periods within parameterizations ( Figure 8c). This indicates that model sensitivity to changes in climatic conditions is dependent on the parameterization for some regions. Assessing the sensitivity in the response to climate forcing to parameter settings would require a larger ensemble of parameterizations and is beyond the scope of this study. However, we demonstrate that parameter sensitivity experiments used in parameterizing DVMs to simulate historical vegetation states may not adequately assess uncertainty in the modeled response to changing climate conditions. The response to climate forcing may be more sensitive to parameters related to the interactions between climate and growing conditions (e.g., temperature limits of photosynthesis, phenological temperatures, nitrogen concentrations, and rooting depths). To adequately quantify uncertainty in projected changes in vegetation, DVMs need to be subject to thorough sensitivity analyses that account for the sensitivity to forced response. Efforts toward establishing benchmarking data (Collier et al., 2018) and metrics (Kelley et al., 2013) that take into account internal variability and trends may help constrain the uncertainty in future projections (Fisher et al., 2018). Parametric uncertainty that is not taken into consideration in future climate projections could contribute considerably to climate feedback (Koven et al., 2015;Sanderson et al., 2015) and the evolution of the carbon cycle (Booth et al., 2012;Fisher et al., 2010;Huntingford et al., 2009;Pugh et al., 2018).
We demonstrated that the parameterization of DVMs contributes a considerable (and often unaccounted for (Koven et al., 2015)) amount of uncertainty to projected change in vegetation state. We assessed the uncertainty in modeled vegetation response to nonstationary climate attributable to parameterization and found notable differences in simulated vegetation transitions among parameterizations that performed similarly under preindustrial conditions. Simulations were performed without large-scale disturbances, such as wildfire or insect attack, which would alter the magnitude and possibly direction of the projected changes (Harris et al., 2016). When combined with large-scale disturbances, the differences among parameterizations may be amplified or reduced. Interactions between disturbances and model parameterization may lead to further uncertainties in transformational vegetation changes under future climate conditions.