Parametric Controls on Vegetation Responses to Biogeochemical Forcing in the CLM5

Future projections of land carbon uptake in Earth System Models are affected by land surface model responses to both CO2 and nitrogen fertilization. The Community Land Model, Version 5 (CLM5), contains a suite of modifications to carbon and nitrogen cycle representation. Globally, the CLM5 has a larger CO2 response and smaller nitrogen response than its predecessors. To improve our understanding of the controls over the fertilization responses of the new model, we assess sensitivity to eight parameters pertinent to the cycling of carbon and nitrogen by vegetation, both under present‐day conditions and with CO2 and nitrogen fertilization. The impact of fertilization varies with both model parameters and with the balance of limiting factors (water, temperature, nutrients, and light) in the pre‐fertilization model state. The model parameters that impact the pre‐fertilization state are in general not the same as those that control fertilization responses, meaning that goodness of fit to present‐day conditions does not necessarily imply a constraint on future transient projections. Where pre‐fertilization state has low leaf area, fertilization‐induced increases in leaf production amplify the model response to the initial fertilization via further increases in photosynthesis. Model responses to CO2 and N fertilization are strongly impacted by how much plant communities can increase their rates of nitrogen fixation and also directly affected by costs of N extraction from soil and stoichiometric flexibility. Illustration of how parametric flexibility impacts model outputs should help inform the interpretation of carbon‐climate feedbacks estimated by in fully coupled Earth system model simulations.


Introduction
The response of the terrestrial biosphere to carbon dioxide enrichment and increasing nitrogen (N) deposition represents a major source of uncertainty in predictions of the evolution of the climate system (Friedlingstein et al., 2006;Arora et al., 2013;Friedlingstein et al., 2014). In the CMIP5 multimodel intercomparison, the smallest response to CO 2 forcing was predicted by the Community Land Model v4 (CLM4; Thornton et al., 2007;Lawrence et al., 2011), the only model in that comparison that represented the impacts of nitrogen availability on plant growth, leading to suggesting that the addition of nitrogen limitation might substantially reduce the capacity for vegetation to respond to CO 2 fertilization (and hence reduce the degree to which the terrestrial biosphere can mitigate increasing atmospheric CO 2 . The newly released CLM version 5 (CLM5;  includes numerous updates relative to the CLM4 and CLM4.5 (Oleson et al., 2013). The introduction of several new dynamical features, including substantial modification of the vegetation nitrogen cycling, a new plant hydrodynamics scheme (Kennedy et al., 2019), altered surface hydrology, snow physics, and dynamic crop representation, has implications for responses to environmental forcing. Wieder et al. (2019) report the global carbon-cycle implications of these changes and note that, in CLM5, carbon uptake responds more strongly to CO 2 and less strongly to N fertilization than its predecessors, with both fertilization responses being closer to experimental results than the previous versions of the model.
Large quantities of analysis and computing resources will likely be devoted to assessing the impacts of transient CO 2 and N on CMIP class models, such as CLM5, in the context of global biogeochemical feedback analysis (e.g., Arora et al., 2013;Friedlingstein et al., 2014). Land surface models, in particular the elements related to biogeochemical cycling, are nonetheless subject to very high degrees of parametric and structural uncertainty Bonan & Doney, 2018) Given this, we consider that the illustration of how parametric choices affect the responses to biogeochemical forcing in such models is an important feature of model documentation.
In this paper, to that end, we give an overview of the new components of the biogeochemical cycling in CLM5 and assess the sensitivity of the CLM5 model parameterization of ecosystem biogeochemistry. We use a range of single site simulations, combined with one-at-a-time parameter uncertainty analyses, and CO 2 and N fertilization experiments, to investigate the dynamics of the model under steady state and step-change modifications in forcing. We further highlight areas of process and parametric uncertainty that appear to require greater attention in future model development and testing cycles.
The majority of land surface model sensitivity tests, uncertainty quantification efforts, and calibration exercises are conducted under present-day conditions Fer et al., 2018;Lu et al., 2018;Ricciuto et al., 2018). Predicted global biogeochemical feedbacks, however, are primarily related to the response of the system to forcing, and not necessarily to those system properties that control the mean state. Thus, our motivating questions here are: 1. Are parameters that dominate responses to CO 2 and nitrogen fertilization also dominant over the pre-fertilization model state? 2. How do fertilization impacts and parameter sensitivities vary with site conditions, as expressed via water and nutrient limitations?
The analysis herein does not focus on the relative fit of the model to the data collected at the individual sites. The fit at individual sites of the globally adjusted model typically cannot be seen as indicative or otherwise of its skill in the absence of specific parameterization using the local vegetation, soil, and hydrological characteristics. A focus on model fit at individual sites would preclude detailed consideration of how parameter variation impacts fertilization responses, as presented here. Assessment of the skill of the CLM5 model against a broad suite of globally relevant data products and manipulation experiments is the focus of  and Wieder et al. (2019). Future studies will consider the impact of alternative parameterizations, informed by these analyses, on the global performance and skill of CLM5.

CLM5
Version 5 of the CLM was released to the research community in 2018 as part of the CESM2.0 release, and open-source releases and ongoing development are housed at the GitHub website (https://github.com/ ESCOMP/ctsm/). An overview of model developments is presented (within this special issue) by , a description of the global nitrogen cycling responses to fertilizaion by Wieder et al. (2019). The plant hydraulics model is documented by Kennedy et al. (2019). Further manuscripts will describe the crop, urban, and land use models, sensitivity to forcing and calibration of fast-timescale processes.

Updates to Biogeochemistry Cycling in the CLM5
The CLM5 biogeochemistry scheme builds on the initial implementation of carbon-nitrogen interactions in CLM4 (Thornton et al., 2007;Lawrence et al., 2011), and on modifications to the biogeochemical cycling (notably, a vertically stratified soil decomposition processes and a slowing of the soil mineral N loss pathways) included in the CLM4.5 (Bonan et al., 2012;Koven et al., 2013). The full technical description of the CLM5 is available as an appendix to  and also online (https://escomp.github.io/ ctsm-docs/doc/build/html/tech_note/index.html). We refer to this document here, in lieu of repeating large sections of the documentation an as appendix.
In this paper, we focus on the response to CO 2 and nitrogen deposition, which are primarily mediated by the biogeochemical processes within the model. The primary changes to the biogeochemical cycling in CLM5 are within the vegetation nitrogen cycle. In particular, the new model contains three new prognostic elements that are likely to impact responses to fertilization.
First, the Leaf Utilization of Nitrogen for Assimilation module (LUNA) simulates distribution of N between different leaf assimilation processes (Ali et al., 2016;Xu et al., 2012). Most land surface models predict photosynthetic parameters directly from leaf nitrogen content ( (Ainsworth & Rogers, 2007), changing temperature , irradiance (Niinemets & Kull, 1998), and soil moisture (Keenan et al., 2009). LUNA predicts the optimal partitioning of N among the maximum rate of carboxylation, V c,max , the maximum rate of electron transport, J max , and other leaf N components, for given leaf nitrogen per unit area (N leaf ) and environmental boundary conditions, assuming colimitation of photosynthesis by both electron capture and carboxylation processes. V c,max is thus primarily controlled by N area and the prevailing meteorology.
Second, the Fixation and Uptake of Nitrogen (FUN) model was incorporated into CLM5. FUN simulates the dynamics and carbon economics of N acquisition from the environment, building on (Brzostek et al., 2014;Shi et al., 2016). A major issue confronting N cycle models is how to manage "excess" carbon under conditions where carbon assimilation is faster than that which can be supported by nutrient uptake . Previous versions of CLM managed this by reducing the gross photosynthetic flux down to the level at which growth could be supported by N uptake. This process occurred after the calculation of stomatal conductance (which is linked to assimilation rates via the Ball-Berry model) and therefore created inconsistencies between C assimilation and water cycling, especially under conditions of N limitation, as discussed extensively in the literature (Medlyn et al., 2011;Walker et al., 2014). The FUN model, in contrast, employs the principle that for the alternative sources of nitrogen acquisition by plants (active uptake from soil, retranslocation from senescent tissues, and symbiotic N fixation), there is a concurrent cost in terms of carbon (Bloom et al., 1985;Rastetter et al., 2001;Jiang et al., 2017;Terrer et al., 2018), which is determined by the abundance of nitrogen in the soils and plant tissues and the temperature-dependant cost of N fixation. Within FUN, plants spend excess carbon on nitrogen uptake by preferentially using pathways that are "cheaper" to them. Where N is scarce in the environment, carbon that would otherwise have been used for plant growth can be deployed to acquire (pay for) more N. Thus, N limitation leads to both reduced growth and higher expenditure of C on N uptake. The switching of N uptake to cheaper pathways results in the preferential use of symbiotic N fixation when soil mineral N concentrations are low (Vitousek et al., 2002), whereas CLM4 and CLM4.5 predicted N fixation as a linear function of net primary productivity (Cleveland et al., 1999;Wieder et al., 2015). The FUN model is documented within the CLM technical note (at https://escomp.github.io/ctsm-docs/doc/build/html/tech_note/ FUN/CLM50_Tech_Note_FUN.html).
Third, CLM5 contains the capacity to modify tissue C:N ratios with spatial and temporal variation in the cost of N uptake. Plant trait databases (e.g., Kattge et al., 2011) typically indicate a wide variability in tissue N concentrations within a single plant functional type (Wright et al., 2005). Ghimire et al. (2016) introduced the capacity for variable tissue C:N ratios, whereby the total nitrogen supply in each time step is partitioned between tissues in proportion to their relative "demand" (ascertained from the target C:N ratio and the carbon allocated to the pool). The FUN model, however, was originally conceived in the context of a static tissue C:N ratio in that C is initially allocated between growth and N uptake in a manner that exactly tracks a target C:N ratio. To allow a variable C:N ratio in the context of FUN, it was necessary to modify how much C is expended on N uptake versus on growth from this default calculation. The degree to which tissue C:N content varies as N becomes limiting is unclear in the literature, and N cycle models typically include semi-heuristic and/or bounded representations of tissue N adjustments under N limited conditions (Zaehle & Friend, 2010;Ghimire et al., 2016). Here we include a placeholder algorithm that modifies C:N ratios as a function of (1) the environmental cost of N acquisition and (2) the degree to which the plant is already far from its target N content. When N is scarce in the environment, less C is spent on uptake, increasing C:N ratios. When plants become N limited and C:N ratios increase, however, expenditure increases. This algorithm is also documented in the CLM5 technical description (https://escomp.github.io/ ctsm-docs/doc/build/html/tech_note/index.html) and includes three tuning parameters the sensitivity of which is investigated here.

Sites
We utilized four observational flux tower sites, in order to capture a range of environmental conditions, across tropical evergreen broadleaf forest (Caxiuanã, CAX −1.   for NWR. Conducting and analyzing these parameter sensitivity experiments globally is beyond the scope of this study. Instead, our aim is to provide detailed illustrations of how model parameter choices impact a model outputs under a range of realistic forcing conditions.

Simulation Protocol
We spun up the default version of the model for 400 years at each site, in "accelerated spin-up" mode (whereby the decomposition of slower cycling carbon and nitrogen pools is increased for the duration of the spin-up, and the pool sizes are modified accordingly at the end of the spin-up-see . We used preindustrial CO 2 concentrations (274 ppm) and recycled the available site-level half-hourly meteorological inputs for the duration of the simulations. Four hundred years was sufficient to bring all soil and vegetation carbon and nitrogen pools into equilibrium. Subsequent to this, parameters were perturbed (see below), and a second spin-up conducted for 180 years for each site, parameter, and value combination, with accelerated mode turned off. Restarting from the state at the end of the perturbed parameter spin-ups, we ran "control" transient simulations for each site and parameter combination starting in 1760 and invoking historical transient CO 2 and N deposition until 2018. For each site/parameter/value combination, we also ran one elevated CO 2 and one elevated nitrogen deposition simulation. Following Wieder et al. (2019), who conducted global runs with the default parameterization of the model, we increased CO 2 to 550 ppm, starting in 1998, and for Nitrogen, we added 5 g N·m −2 ·year −1 , also starting in 1998.

Perturbed Physics Ensemble
The CLM5 has a large array of parameters that might have direct and indirect impacts on biogeochemical cycling in the model. Here we focus on a set of eight parameters, chosen for their direct mechanistic impacts on responses to CO 2 and N fertilization. These parameters are narrowed down from the much larger complete set of CLM5 parameters using a process informed primarily by the model structure, as well as the extensive iterative investigations of the CLM5 model during the development process, as well as a prior investigation of 14 candidate parmeters. We focus on biogeochemically important inputs; parameters relevant to the hydrology and surface energy balance processes (with the exception of the stomatal slope) will be assessed by studies currently underway, and parameters of the new plant hydraulics scheme are discussed in Kennedy et al. (2019).
To reduce the complexity involved in understanding and attributing changes in model behavior with parameter variation, we conducted one-at-a-time (OAAT) perturbations of the parameters around the default state, which (compared to a global parameter sensitivity analysis) allows more straightforward visualization and interpretation of the results. At each site, we conduct simulations with five values of each parameter (one of which is the default) across the ranges described below. We note that a global parameter sensitivity analysis would also be of interest, and indeed have conducted many such analyses, but concluded that OAAT analysis provides clearer insight into the functionality of the parameters with respect to biogeochemical forcing.

Default Parameter Determination
During the development of CLM5, we conducted extensive testing into the possibility of inverse calibration of the model parameter space using neural-network-derived emulators. We did not use the results of this effort for a variety of reasons, including, (1) the high dimensionality of the parameter space; (2) the degree of non-orthogonal behavior in the parameter response surfaces; (3) the tendency of the modeled ecosystems to die under lower pre-industrial CO 2 conditions when calibrated using present-day CO 2 ; (4) the dominance of high-productivity areas over the calibration of a particular Plant Functional Type (PFT) at the expense of marginal areas, resulting in marginal areas with many dead grid cells; (5) the subjectivity in weighting of the various observational data products; and (6) philosophical issues concerning whether the structural validity of the calibrated model could be independently assessed. The parameterization used in the default version of the model consequentially was derived from the mean observed trait values where possible, in combination with targeted manual tuning of parameters (including allocation parameters and nitrogen costs) that were less well constrained by data, where excessively unrealistic values for simulated states were detected. We recognize that this more subjective calibration method is less than ideal. Ongoing studies will focus on more limited-scope calibration attempts to overcome the difficulties described (e.g., Fer et al., 2018).

Parameters Selected for Perturbation Experiment
The parameters included in our sensitivity analysis are listed in Table 1. Their impact on model functionality is described below. The CLM5 shorthand version of the parameter name is given in italics. Note. For the Medlyn slope parameters, PFT-specific multipliers of the standard errors (se) from  were used, where the standard errors were 0.25 for evergreen needleleaf, 0.09 for evergreen broadleaf, 0.25 for deciduous needleleaf, and 0.36 for deciduous broadleaf trees. v1-v4 correspond to the multipliers used in the four sensitivity analysis runs (two higher and two lower than the default).
Specific leaf area, SLA, determines the area of leaves (in m 2 ) derived from 1g of leaf biomass. High SLA values correspond to thinner, "cheaper" leaves, which lead to higher leaf area index for a given leaf investment. In CLM5, high SLA leaves also have lower N per unit area (given a target mass-based C:N ratio) and thus lower photosynthestic capacity per unit area.
Fine root mass per unit leaf mass, froot_leaf , is the ratio of fine root biomass to leaf biomass. In CLM5, fine root biomass affects the capacity of vegetation to acquire both water and nutrients via the resistance of the rhizosphere to water flow and the cost of N acquisition, respectively. Lower values should in principle incur a cost to vegetation in terms of below-ground resource acquisition, while higher values imply greater construction and maintenance costs.
Stem to leaf mass ratio, stem_leaf , controls the fraction of biomass to stem tissues as a fraction of the leaf allocation. Increases in this parameter both decrease the leaf area index (LAI, m 2 /m 2 ) achievable for a given amount of (C and N) allocation to growth, and increase the equilibrium woody biomass. Woody biomass in CLM5 has no functionality, and thus, there is no cost to plants associated with low values of this parameter.
Nitrogen costs, n_costs, are a set of parameters that determine the environmental cost of nitrogen uptake from the soils. FUN uses a bulk approach to account for the various costs to plants of low environmental N availability, including carbon expended on root exudates, symbiotic relationships with soil fungi, expenditure on fine root growth, and/or increased metabolic rates to support active uptake. This abstraction means that uncertainties associated with the cost functions, in particular, are large. These six parameters were defined and calibrated by Brzostek et al. (2014). In CLM5, we maintain the ratios of the parameters as defined by the previous FUN model calibration but allow the magnitudes of the parameters to change concomitantly, and with plant type, since the definitions of soil N availability vary between CLM5 and Brzostek et al. (2014). The six n_costs parameters include the rate constants that relate the cost of ectomycorhhizal uptake to the fine root carbon pool (ekc_active) and soil N concentration (ekn_active), respectively, and their counterparts akc_active and akn_active for arbuscular mycorrhizal uptake and kc_nonmyc and kn_nonmyc for nonmycorrhizal active uptake.
The fraction of carbon that can be used for fixation, fracfixers, is the fraction of the net primary production (NPP) or gross primary productivity (GPP) minus autotrophic respiration (not including carbon spent on N uptake) that can be used for symbiotic N fixation. It is thus a proxy for the ecosystem-level fractional productivity of symbiotically N-fixing plants. The constant fracfixers parameter here acts in lieu of, for example, a prognostic model of N fixer distributions or an input map of their abundance, both of which could plausibly be added in future generations of CLM.
Leaf carbon:nitrogen target, leafcn uses the concept of a "target" leaf C:N ratio, around which flexibility is allowed, given the environmental costs of nitrogen. As with most land surface schemes, higher leaf N concentrations lead to both higher maintenance respiration costs and higher photosynthetic capacity.
Slope of the stomatal conductance model, medlyn_slope. CLM5 replaced the Ball-Berry stomatal conductance scheme with the Medlyn stomatal conductance scheme (Medlyn et al., 2011). The slope parameter determines the degree of stomatal opening for a given combination of assimilation capacity, CO 2 , and vapor pressure deficit, in the absence of plant hydraulics limitations. It is included here on account of its critical impact on stomatal responses to CO 2 .
C:N flexibility, the cn_flex_a, cn_flex_b, & cn_flex_c, parameters as described in the CLM5 technical note, determine the response of tissue C:N ratios to depletion of N in both the environment and the plant tissues themselves. Only the third parameter (c), which determines how tightly plants retain C:N ratios close to their target, had an important impact in this sensitivity analysis, and so we focus here on the impact of that parameter.
Other variables we tested in earlier versions of this analysis included rate constants controlling the growth and maintenance respiration rates, the fraction of ectomycorrhizal fungi, and parameters "a" and "b" of the flexible C:N ratio model. None of these had important impacts on the fertilization responses, and growth and respiration rates affected the model state in a linear and predictable fashion. They are omitted here for the sake of simplification, but model output related to their impacts can be accessed via the available simulations.

Parameter Ranges
SLA and leafcn are readily observable and thus among the most common traits represented in the TRY plant trait database (www.try-db.org). We derive values from the analysis by Kattge et al. (2011) for their mean and distribution. For both parameters, values are log-normally distributed. We here take the PFT average log-normalized standard deviation range from TRY and apply it across all PFTs (thus, the lower bounds are closer to the mean than the higher bounds; Table 1). The range of fracfixers was set to vary from 0% to 100% of the net productivity being available for fixation. medlyn_slope varies across the standard errors reported in de Kauwe et al. (2015). Reasonable ranges of the cn_flex parameters are unknown, since these physiological controls are poorly understood. Thus, we explore an order of magnitude in each direction from the default, to obtain an understanding of the general sensitivity to this parameterization. The same issues apply to N costs parameters.

Responses to Forcing
We calculated the impacts of the CO 2 and N deposition by assessing the ratio of the fertilized to control simulations, averaging years 15-19 after fertilization began. For some parameter values, the constraint on growth placed the control system in a state close to or at death (zero live carbon pools). Addition of N or CO 2 therefore caused either a very strong recovery of the plants, or no recovery at all, where biomass was already zero, leading to some large non-linearities in the relationship between parameters and environmental forcing. To focus on the responses of parameter combinations that give more reasonable initial vegetation conditions, parameter combinations that caused LAI to be less than 60% of the control simulation value were excluded from the fertilization plots.

Overall System State
In this section we assess first the impacts of parametric variation on the pre-fertilization system state. The parameter responses of the site-level water and nitrogen limitation (see below) in the pre-fertilization parameter ensemble are shown in Figure 1 and GPP, NPP, LAI, and N area in Figure 2. In each of these figures, the y axis corresponds to the output variable, and the x axis represents the range of parameter perturbation with −1 to +1 corresponding to the range of variation in each parameter, depicted in Table 1. Lines of differing colors illustrate responses to variation in the target parameters, with the central crossover point in each figure corresponding to the default case.

Water and N Limitation Across Sites
The four sites varied in the degree to which they are limited by water and N availability in the pre-fertlization state (Figure 1), which in turn impacts the parameteric controls over sensitivity to forcing. Water limitation in CLM5 can be expressed via a model property termed B tran , which is calculated as a function of the (prognostic) leaf water potential. B tran is used to reduce assimilation and thus stomatal conductance by scaling the photosynthetic capacity, V c,max . Leaf water potential and B tran vary diurnally, so we use the average daily minimum (B tran,mn ) here as a metric of drought stress (1 is no drought stress; 0 is total stomatal closure). For N limitation, the percentage of total NPP (where NPP is GPP-r auto ) that is used to pay for nitrogen uptake (PNU) is a direct indicator of how N stress affects the carbon budget. Thus, the sites chosen span a range of initial limiting conditions. Underlying limitations by different boundary conditions (water, nutrients, light, and temperature) affect the ability of the system to respond to fertilization and/or parameter modifications, as we describe further in the following sections.

Parametric Control Over System State
In this section, we discuss the impacts of the parameter variation on the initial model state. Many parameters have both first-order direct impacts, and potential indirect effects via modification other elements of the model. Here we attempt to describe how these mechanisms operate in CLM5 within the parameter space explored here.
Specific leaf area (slatop) has a large first-order impact on LAI , and on N leaf (Figure 2), at all sites. slatop is used directly in the calculation of both of these quantities. The impact on GPP is always positive, despite the thinner leaves having lower N, and thus lower photosynthetic capacity, per unit area of leaf. The impact of SLA on NPP can, however, sometimes be negative (CAX and NWR), since the fixed-fraction allometry model in CLM5 can generate shaded leaves that are in negative carbon balance at high LAI.
Root:leaf ratio (froot_leaf ) has numerous impacts, including (1) on LAI, via tissue allocation fractions; (2) on NPP, via changes in maintenance costs of fine root pools; (3) on the uptake of water, via root length density; and (4) on nutrient limitation, via root-mediated N uptake costs. The first two impacts represent the costs of roots, and the latter two their benefits. CAX and HVF exhibit higher LAI with lower root allocation, Figure 2. Response of control model states and fluxes to parameter perturbation across −1 to +1 range of parameter variation (Table 1) for gross primary productivity (GPP), net primary production (NPP), leaf area index (LAI), and N leaf at the CAX, NWR, HVF and MET sites.
whereas NWR has large reduction in LAI (and GPP and NPP) with lower root allocation. MET exhibits a mixed response crossing an optimal value. froot_leaf is a primary control over drought stress at both NWR and MET ( Figure 1), illustrating the impact of low root density on water uptake capacity.
Low values of stem_leaf ratio resulted in higher LAI at all sites, as well as declines in total vegetation carbon ( Figure S3 in the supporting information). The impact of increasing LAI is typically to increase light interception and hence GPP; thus, these responses are typically quite large for most sites. Lower stem allocation and higher LAI also increase N demand, reducing N leaf and increase PNU at NWR and MET (Figure 1).
Increasing N costs increases the respiration losses associated with nitrogen uptake. All sites showed a moderate decline in LAI, NPP, GPP, and vegetation carbon when costs were high, but lesser impacts when costs were reduced ( Figure 2). In contrast to the allocation parameters, N costs does not appear to be a first-order control over the pre-fertilization system state.
Increasing fracfixers affects the maximum amount of carbon that plants can pay toward symbiotic nitrogen uptake (N fixation). For all sites and variables here, it has a relatively moderate and positive impact on the control state. As is the case for N costs , fracfixers does not appear to be a first-order control over the prefertilization system state. Interestingly, Figure S1 illustrates that while fracfixers had a large impact on the N fixation rates, this did not actually affect the overall soil nitrogen content, implying that other factors were dominant over the equilibrium soil nutrient budget.
Changing leafcn over one standard deviation from the mean TRY values had substantial impacts on all the output variables at all sites. Higher leafcn values decrease N leaf , decreasing photosynthetic capacity and GPP.
In particular, at MET and HVF the higher leafcn values had N leaf and photosynthetic capacities so low that GPP was almost zero (Figure 2). GPP continues to increase almost linearly with declining leafcn (increasing N leaf ) at all sites except HVF, indicating a lack of saturation of photosynthesis with respect to photosynthetic capacity. This is consistent with the analysis by , who find that V c , max numbers in CLM5 might be low compared to observations. For NPP, respiration costs of modified N leaf and altered PNU ( Figure 1) buffered the larger changes in GPP. At CAX, for example, the lowest leafcn had a GPP 225 g C·m· −2 ·year −1 higher than the default, while the corresponding NPP increase was only 52 g C·m· −2 ·year −1 .
Changes in medlyn_slope, the stomatal slope parameter, across the range shown here, had limited impacts on the prefertilized system state, with the exception of NWR, where lower values of the slope parameter generated higher photosynthesis, consistent with this site experiencing severe drought stress and thus benefiting from increased water efficiency.
cn_flex_c affects the degree to which plants attempt to maintain tissue C:N values near to their target when environmental costs are varying. For the pre-fertilization state impacts here, the role of this parameter was limited, but for the step-change fertilization it was more important, so we discuss its role further in the section below.

CO 2 and Nitrogen Deposition Responses 3.2.1. Timescale Dependence of Responses to Fertilization
The initial changes in carbon economy induced by the fertilization experiments are, where conditions permit, amplified by increases in LAI. Figure 3 illustrates, for the default parameter set, the evolution through time of GPP, NPP, LAI, B tran,mn , PNU, N fixation and soil Nitrogen, illustrating the complexities of the response.
At MET (purple lines), for example, the initial GPP stimulation by CO 2 (26%) resulted in a slow but prolonged LAI increase, (given the lower overall carbon turnover at this colder site). At 80 years, LAI had increased by 50%, amplifying the much GPP and NPP stimulation to 53-69% and 49-71% respectively, (varying through the cycling climate forcing data). For CO 2 fertilization at CAX (blue lines), in contrast, the initial stimulation of GPP (∼1.37) lead to an almost 50% increase in LAI, but at this site, water limitations (lower B tran,mn ), prevented any further increase in GPP. At HVF (yellow), the initial stimulation in the first year was lower (25%) potentially on account of the higher carboxylation capacity at this site (HVF has a saturating response to decreasing leafcn, Figure 2). NWR (red) had moderate (25%) increases in initial GPP and LAI (∼30%). The stimulation of GPP via LAI here is typically greater than is typically reported in elevated CO 2 experiments, but King et al. (2005) report a large response to ∼530 ppm CO 2 in Aspen stands, with increases in LAI and biomass of around 60%.
At all sites, the percentage of NPP spent on N uptake (PNU) increased with CO 2 fertilization, but most dramatically at MET, where initial N expenditure was lowest initially on account of low growth rates. Further, N fixation rates increased markedly at all sites ( Figure 3) from their default values. An increase of x2-3 was found at HVF and NWR. At MET, the very low initial rates (0.01 gNm −2 y −1 ) mean that the amplification was extremely large and is thus omitted from the plot. At CAX, fixation rates initially increased and then declined through time in tandem with increases in soil Nitrogen. Soil N in general was not in equilibrium after 85 years of simulation time, but this did not appear to greatly affect the fixation rates at NWR and HVF. The soil N dynamics observed raise interesting questions about potential transient responses of CLM5, which will be covered in future analyses.
N addition most strongly affects carbon expenditure on N uptake, and thus NPP, rather than GPP directly. At NWR, an initial NPP stimulation of 1.3 was magnified by subsequent increases in LAI to nearly 1.8 times the control value, causing an overall NPP increase of 1.7, reflecting the high N limitation at this site in the pre-fertilization state. None of the other sites exhibited such a dramatic change. In particular, for CAX, which was heavily water limited in the first instance, adding nitrogen only resulted in stimulation in LAI, GPP and NPP of around 10%. At all sites, addition of N caused a radical drop in the rates of N fixation, and gradual, non-equilibrated increases in the soil N content.
The time dependence of the responses generates an interesting issue for interpretation. The first year responses are illustrative of the short-term immediate impacts of change, whereas the longer-term responses are heavily dominated by the feedback responses. Further, depiction of absolute versus relative differences can serve to highlight different aspects of the model response. Here we report values from 15-20 years after fertilization began, to be consistent with the methodology used by Wieder et al. (2019) in their analysis of global fertilization of the default model, and to capture the dynamics of the whole system feedbacks captured in Figure 3. Thus, Figures 4 and 5 illustrate the responses to 15-20 years of CO 2 and N fertilization simultaneously, for GPP and LAI. Note that each figure is a different scale to best illustrate the parametric impacts. Corresponding figures for NPP, vegetation carbon, N fixation, and N area are in the supporting information ( Figures S2-S5).

Impact of Parameters
The parameter response space of the model is complex, and parameter impacts on fertilization responses vary with changing initial conditions and with treatment. Here, we discuss the major impacts of each parameter on the responses and how they relate to expectations based on model structure.
The impact of slatop on fertilization was different at every site (Figures 4 and 5). For MET, slatop was the most important control over N responses, with higher values allowing a greater (or faster) LAI/GPP feedback in response to N fertilization. At NWR, the opposite behavior was found, where high SLA reduces the impact of N. At this site, high pre-fertilization LAI meant that the response of GPP was more saturated, and so the feedback from increasing leaf area was dampened. At HVF, higher values of slatop slightly reduced both CO 2 and N fertilization for the same reason (although the range of variation at HVF is only 5% between the highest and lowest values, as opposed to 40% at NWR).
froot_leaf had contrasting impacts at different sites. At MET, low values increased the CO 2 fertilization, probably due to the increased drought stress under those conditions (Figure 1). At HVF, low froot_leaf decreased both CO 2 and N fertilization on account of higher initial LAI, mirroring the response to slatop.
Low stem_leaf ratio caused lower responses to N and CO 2 fertilization CAX and HVF for GPP and LAI at (Figures 4 and 5), again mirroring the NPP-LAI nonlinearity issue. At NWR, N limitation restricted CO 2 responses, nullifying the impact of stem_leaf . In contrast, at MET, lower stem_leaf allowed greater responses to N plausibly from increasing the speed of the NPP/LAI feedback.
Higher N costs increase the responsiveness of the system to N at CAX, NWT, and HVF. The magnitude of the impact on N fertilization tracks the extent of N limitation in the pre-fertilization state (Figure 1). Very low values of N costs appeared to reduce the responsiveness of the system to CO 2 fertilization, particularly at CAX and HVF. These configurations potentially are limited in their ability to direct carbon to fixation pathways, since estimated soil N uptake costs remain low.
fracfixers (the maximum fraction of NPP available for N fixation) had a particularly dominant role over the model response to step changes in N and to CO 2 . For GPP and LAI, all sites displayed a trade-off, whereby high fracfixers systems responded strongly to CO 2 but not N, and low fracfixers systems responded strongly to N but not CO 2 (Figures 4 and 5). At MET, the changing fixation capacity had no impact on the response of GPP or LAI to N fertilization. Only slatop and stem_leaf had large impacts on N fertilization for this site (which had the smallest PNU in the first instance), suggesting fundamental metabolic limitations to carbon economics and growth (e.g., a short growing season).
The impact of leafcn on fertilization responses was not as dominant a control over fertilization responses as it was over the system state. Complex responses were illustrated, reflecting the differential costs and benefits of leafcn on the system. For N fertilization, at CAX and NWR, lower values of leafcn had smaller responses to N fertilization, reflecting the higher initial LAI of these configurations and the nonlinear GPP feedback discussed earlier. At CAX, the highest leafcn had smaller CO 2 responses. Figure 1 illustrates that this configuration had less water stress (both lower LAI and lower photosynthetic capacity should lead to smaller stomatal conductances and evapotranspiration) and thus had less water stress to be alleviated by the high CO 2 .
Surprisingly, given the importance of medlyn_slope on stomatal responses to CO 2 , we did not find a widespread impact of the Medlyn parameter on the fertilization response. Earlier studies using a less constrained range of medlyn_slope values indicated that it can have a potentially large impact on system state via impacts on water status (results not shown) but that its capacity to modify the fertilization response is limited. The linear response of the stomatal conductance term to CO 2 means that the direct impact of medlyn_slope is likely to be restricted to changes in the prefertilization state.
The impact of cn_flex_c was small compared to the more dominant parameters. Higher values of cn_flex_c imply less carbon spent on N uptake when tissue N is lower than the target amount, allowing N ratios to vary more. Impacts of this parameter are manifested through their impact on N leaf ( Figure S5). For N deposition, N leaf always increases, and for CO 2 it always declines. Smaller values of cn_flex_c corresponded in all cases to smaller variations in N leaf under both CO 2 and N fertilization. Lower values of cn_flex_c corresponded to larger CO 2 stimulation of GPP, as these leaves maintain higher photosynthetic capacity under more N limiting conditions. The impact on NPP and LAI was less clear, given the consequence of high N leaf for both respiration and N uptake costs. At CAX, NWR, and HVF, the very highest values of cn_flex_c also had a higher response to CO 2 than the default. The reasons for this discontinuity are unclear; however, the impact of cn_flex_c in the CLM5 code is subject to at least two "caps," (whereby plants cannot spend less than 0.5 or more than 1.0 of the carbon that would retain the target LAI), potentially leading to complex feedbacks once these caps are reached.

Discussion
The responses of the land surface components of Earth system models to CO 2 fertilization are the subject of intensive investigations, and vast quantities of computing time will be expended on the release version of each model, under the auspices of the various CMIP6 model intercomparison projects and associated activities (Meehl et al., 2014;Eyring et al., 2016). Throughout the CMIP activities, the carbon-cycle feedbacks reported by each modeling group will be from a single default instance of the parameter space. Interpretation of these simulations should be informed by solid comprehension of which model properties dominate the size of the feedback strengths, both within and between models and, ideally, a sense of how well constrained these properties are. We intend these illustrative simulations to provide initial guidance of this process. At the time of writing, we do not yet know the carbon-climate feedback strength of the default version of CLM5/CESM2. In off-line simulations forced with a data atmosphere, CLM5 shows a stronger response to elevated CO 2 and a weaker response to N fertilization than previous versions of the model (Wieder et al., 2019). Further, the impact of both types of fertilization was found, in all versions of the model, to be highly variable with PFT (from 0% to 110% for N and −10% to +55% for CO 2 ). In this study we set out to illustrate how variations in parametric space give rise to responses to forcing. While we set out primarily to to determine responses to model parameters, our findings concerning the impact of pre-fertilization model state, water, and N limitation also have important implications for the interpretation of simulations with the default model.
Given the dominant impact of N fixation (fracfixers) on the model fertilization responses, it seems likely that the new capacity, introduced in CLM5, for plants to use carbon to pay for the fixation of addition nitrogen under CO 2 fertilization is at least partly responsible for the change in the global responses and also that any changes in this parameter would be likely to have large impacts on this result.
Aside from the straightforward direct impacts of N fixation rates on model predictions, interpretation of other configurations is less straightforward. In some cases the sign of the impact of individual parameters changes between sites, while other parameters can demonstrate nonlinear and counterintuitive impacts as the result of numerous interacting feedbacks. To summarize our findings, we note a variety of key outcomes of our analysis on the factors controlling responses to CO 2 and N fertilization in the CLM5: 1. The initial model state has a dominant effect on responses to CO 2 and N fertilization. The initial LAI, and whether the system is on a linear or saturated part of the LAI versus GPP relationship, has a large impact on any fertilization response. Parameter impacts in more productive systems are likely to generate smaller LAI/GPP feedbacks and lower parametric variation as a result. Many of the direct impacts of the parameters assessed here were masked by their large impacts on the initial model state.
The clearest examples of this were the slatop and stem_leaf impacts, but this impact was illustrated for many model configurations.
2. The degree of water or nutrient limitation in the initial state is also a dominant factor over the model responses to CO 2 fertilization and N deposition, and parameteric control thereof. Examples of this include how low froot_leaf increases water limitation at NWR, increasing responsiveness to CO 2 , and how high leafcn reduces water stress and thus decreases the CO 2 response at CAX.
3. Other parameters with a direct impact on fertilization responses included those related to the costs of acquisition of nitrogen (N_costs) and the stoichiometric flexibility of the system (cn_flex_c) illustrating the sensitivity of the model to relatively unconstrained processes within the nitrogen cycle.
4. Many of the parameters that are dominant over the initial state (slatop, stem_leaf leafcn) either have limited impacts on fertilization responses or have impacts that are largely traceable to their impacts on the initial model state. In contrast, numerous parameters that more directly and uniformly impact fertilization responses (fracfixers, N_costs) are less important for the initial model state. We thus suggest that model calibration efforts focused on parameters that affect the control model state are unlikely to constrain the major axes of responses to CO 2 and N fertilization.
5. The results can in general be interpreted as being in line with the expected dynamics of the model, given sufficient output data and consideration of the dynamics of the individual sites. Given the complexity of the model, this helps generate confidence that the model outcomes are derived from the intended process mechanisms and are not subject to significant deviations from the intent of the model structure.
6. The large degree of variability in the responses between sites and output variables will likely hinder attempts to generate land surface models with globally optimized parameters. This analysis should help inform and guide the design of such efforts. Sensitivity tests across broader geographical ranges would help to eludicate more clearly how changes in baseline water and N limitation impact parameter responses and spatial patterns in responses to fertilization.

Model Improvements Highlighted by this Analysis
While there are many potential improvements that could be made to any land surface model, this particular assessment of the parameter space of the CLM5 revealed several dynamics within the representation of biogeochemistry which might be improved in future versions.
First, the response (e.g., at CAX) of NPP to increasing slatop was often negative, despite increasing GPP. The fixed-fraction allocation scheme in the model will, under some circumstances, generate leaves that respire more than they assimilate under very high-productivity conditions. This highlights the need for either more realistic tissue allocation, or for a scheme that optimizes leaf allocation to prevent leaf layers in negative carbon balance (e.g., . Second, in other cases, the heuristic C:N stoichiometry algorithm can increase tissue N content in response to environmental forcing in a manner that, again, can act to reduce net assimilation by modifying the balance between assimilation and respiration. In principle, there exists an optimal leaf C:N ratio that maximizes leaf C export under a given set of conditions. Numerous models of this type have been proposed in the plant optimal functioning literature (Anten & During, 2011;Franklin et al., 2012;McMurtrie & Dewar, 2011Thomas & Williams, 2014;Van Wijk et al., 2003), and future model development efforts could revisit this topic.
Third, one issue requiring focus is the dominance of the fraction of fixers over fertilization responses. fracfixers is a proxy for the community level composition of nitrogen fixing plants. Technically, because plant density is not represented in the CLM5, this parameter represents the fraction of net assimilated carbon (gpp -r auto ) that can be used for symbiotic N fixation. Ideally, frac_fixers would be prognostic, and the fraction of N-fixing plants would respond to changes in environmental conditions and the relative competitive abilities of N-fixing plants. In general, nitrogen fixers are known to primarily exist in early successional parts of ecosystems (Vitousek et al., 1989), and therefore, representation of their ecological dynamics would likely include a representation of vegetation demographics and succession (Fisher et al., 2018;Trugman et al., 2016).
Lastly, we recognize the importance of other limiting nutrients, particularly phosphorus, in controlling the influence of CO 2 fertilization (Reed et al., 2015). Phosphorus dynamics were not integrated into this version of the CLM, but alternative configurations of the CLM exist (e.g., Yang et al., 2014;Zhu et al., 2016) that will provide contrasting biogeochemical predictions in future studies.

Applicability of Results to Other Land Surface Models
Our finding that Nitrogen fixation plays a dominant role over responses to CO 2 fertilization concurs with earlier studies. Notably, Meyerholt and Zaehle (2015) investigate six alternative models of N fixation implemented within the O-CN model (Zaehle & Friend, 2010). They conclude that their optimizing algorithm, which, as in CLM5, switches N sources depending on their relative costs, represents a promising means by which N fixation might be simulated. Interestingly, this model also includes a maximum rate of fixation per unit (root) biomass parameter, but analysis to its sensitivity is not shown. Further, given the difficulty of matching observations of complex N cycle processes with more abstract "big leaf" models such as CLM5, and the ecological complexities governing the abundance their abundance, Meyerholt and Zaehle (2015) also advocate for the dynamic simulation of the distribution of N-fixing plants within ESMs, as opposed to the assumption that maximum N-fixing capacity if static in space and time.
Given that the response of GPP is by definition non-linear with respect to LAI (given the saturating properties of light absorbtion in increasingly dense canopies), we anticipate that the impacts of initial state will be similar in most land surface models. While it is perhaps not a novel point, water limitation should act to modulate CO 2 fertilization, while N fertilization will of course be stronger in areas that are more N limited in the first instance.
We find that the model parameters with substantial control over fertilization responses are dissimilar to those with greatest control over the initial model state. However, the CLM5 utilizes a highly simplified fixed-fraction plant allocation scheme, wherein alternative parameterizations can generate very highly divergent ecosystem properties, and further, the allocation fractions are not easily constrained by real-world estimates, which vary with plant size. Other models with alternative allocation schemes might not be so sensitive to these parameterizations. Nonetheless, the majority of the parameterizations affecting the transient responses had limited impact on the pre-fertilization state. This suggests that model behavior in response to rapid change might indeed be timescale dependent and perhaps requires further analysis.

Future Directions
In this study we investigate the impacts of parametric uncertainty on the carbon and nitrogen responses of the CLM5. In particular, we probe a targeted parameter set comprising those parameters that determine the impacts of CO 2 and N fertilization. The full response of the model to transient climate forcing has considerably more dimensions than these two transient properties. In particular, here we have not considered responses to temperature nor to changes in the supply or demand of water. The suite of parameters that in principle control temperature responses of the model are somewhat distinct from the parameters considered here, given the large number of processes that are directly or implicitly affected by temperature (photosynthesis, respiration, soil decomposition, N mineralization, and cryosphere interaction) and should be the subject of further investigations.
Many kinds of further analyses are of course possible, including ramping transient simulations, global simulations with alternative parameterizations, calibration of critical parameters against manipulation experiments, wider assessment against greater ranges of flux tower data, and tighter constraints of the parameter ranges used. We hope that this initial assessment of the CLM5 parameter space will provide a starting point for such activities.
In principle, integrating more targeted data sets and data products should narrow the range of model parameter sets that are acceptable in comparison with the suite of available and relevant data sets, as well as accelerate the rejection of model structures that are unrealistic. The development of the CLM5 slightly predated the operational usage of the International Land Model Benchmarking Project package (Collier et al., 2018; for assessing model skill across a very broad range of model components. Ideally, this effort should include synthesis of ecosystem manipulation experiments that provide altered boundary conditions in real life, given the large differences observed here between the sets of parameters affecting present-day conditions and those that impact transient responses.
It is common practice in land surface modeling literature to present studies that focus on the modification of a single aspect of model structure, process representation, or calibration and to highlight the significance of individual modifications when scaled up to global land-mediated climate feedback processes. In this analysis, we illustrate the degree to which a multitude of individual parameters can potentially have important and interacting impacts on model results, how these results differ across sites and output variables and with initial vegetation conditions. Future studies of process and parameter modification in CLM5 might use these results as a guide to the relative importance of individual focus parameters and processes. We thus hope to move forward from studies that modify individual processes in isolation to those that take a more holistic view of the complexity inherent in modeling coupled land surface systems.

Conclusions
In this paper we present an analysis of a suite of parameters relevant to vegetation-mediated biogeochemical feedbacks in the extensively modified CLM5. We illustrate how responses of the model to elevated carbon dioxide and nitrogen deposition might be impacted by parameter choice. We find, in particular, that parameters controlling the degree to which ecosystems can fix additional nitrogen have a dominant effect over responses to fertilization, but limited impacts on the control state. Further, we find that the impact of most parameters is complex and affected by the balance of water, nutrient and other limitations to growth specific to individual sites.
We have attempted here to illustrate and understand the connections between the structure and parameterization of the CLM5 and its biogeochemical functionality. This necessarily complex exercise represents one small part of the assessment of CMIP6-era representations of the land surface required to understand how advances in model structures alter our predictions of land surface functionality. The challenges of comprehending, measuring and constraining the parametric and structural uncertainty of increasingly complex process-based and land surface models are daunting. We argue that in order to make progress on this task, the Earth system modeling community should move beyond the paradigm of assessing single points in model parameter space in great depth, ignoring all other plausible instantiations. In the atmospheric sciences, for example, the 'internal variability' caused by initial conditions is now routinely sampled (Kay et al., 2015;Sanderson et al., 2018), in light of its impact on seasonal to decadal timescale predictions. Bonan and Doney (2018), however, illustrate how initial conditions are relatively unimportant for carbon-cycle prediction on land, relative to structural model uncertainty. Given this, and the finding here that responses to forcing are sensitive to parametric variation, sampling the important dimensions of parametric uncertainty that govern biospheric feedbacks in the climate system should be taken seriously when planning future assessments of biogeochemical feedback strength.