Multiple soil nutrient competition between plants, microbes, and mineral surfaces: model development, parameterization, and example applications in several tropical forests

Abstract

Soil is a complex system where biotic (e.g., plant roots, micro-organisms) and abiotic (e.g., mineral surfaces) consumers compete for resources necessary for life (e.g., nitrogen, phosphorus). This competition is ecologically significant, since it regulates the dynamics of soil nutrients and controls aboveground plant productivity. Here we develop, calibrate and test a nutrient competition model that accounts for multiple soil nutrients interacting with multiple biotic and abiotic consumers. As applied here for tropical forests, the Nutrient COMpetition model (N-COM) includes three primary soil nutrients (NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$ and PO${}_x$; representing the sum of PO${}_{\mathrm{4}}^{\mathrm{3}-}$, HPO${}_{\mathrm{4}}^{\mathrm{2}-}$ and H${}_{\mathrm{2}}$PO${}_{\mathrm{4}}^{-})$ and five potential competitors (plant roots, decomposing microbes, nitrifiers, denitrifiers and mineral surfaces). The competition is formulated with a quasi-steady-state chemical equilibrium approximation to account for substrate (multiple substrates share one consumer) and consumer (multiple consumers compete for one substrate) effects. N-COM successfully reproduced observed soil heterotrophic respiration, N${}_{\mathrm{2}}$O emissions, free phosphorus, sorbed phosphorus and NH${}_{\mathrm{4}}^{+}$ pools at a tropical forest site (Tapajos). The overall model uncertainty was moderately well constrained. Our sensitivity analysis revealed that soil nutrient competition was primarily regulated by consumer–substrate affinity rather than environmental factors such as soil temperature or soil moisture. Our results also imply that under strong nutrient limitation, relative competitiveness depends strongly on the competitor functional traits (affinity and nutrient carrier enzyme abundance). We then applied the N-COM model to analyze field nitrogen and phosphorus perturbation experiments in two tropical forest sites (in Hawaii and Puerto Rico) not used in model development or calibration. Under soil inorganic nitrogen and phosphorus elevated conditions, the model accurately replicated the experimentally observed competition among nutrient consumers. Although we used as many observations as we could obtain, more nutrient addition experiments in tropical systems would greatly benefit model testing and calibration. In summary, the N-COM model provides an ecologically consistent representation of nutrient competition appropriate for land BGC models integrated in Earth System Models.

Introduction

Atmospheric CO${}_{\mathrm{2}}$ concentrations have risen sharply since the pre-industrial era, primarily due to anthropogenic fossil fuel combustion and land use and land cover change (Houghton, 2003; Le Quéré et al., 2013; Marland et al., 2008). Terrestrial ecosystems mitigate the increasing atmospheric CO${}_{\mathrm{2}}$ trend by absorbing roughly a quarter of anthropogenic CO${}_{\mathrm{2}}$ emissions (Le Quéré et al., 2009). However, it is still an open question whether the terrestrial CO${}_{\mathrm{2}}$ sink can be sustained (Sokolov et al., 2008; Zaehle et al., 2010), given that plant productivity is generally limited by soil nutrients (Elser et al., 2007; LeBauer and Treseder, 2008; Vitousek and Howarth, 1991) and soil nutrients could be quickly depleted through biogeochemical (Chauhan et al., 1981; Nordin et al., 2001; Shen et al., 2011) and hydrological (Dise and Wright, 1995; Perakis and Hedin, 2002) processes. Therefore, a holistic representation of soil nutrient dynamics is critically important to model the responses of terrestrial ecosystem CO${}_{\mathrm{2}}$ uptake to climate change.

Until recently, land models integrated in Earth system models (ESMs) have largely ignored the close coupling between soil nutrient dynamics and the carbon cycle, although the impacts of soil nutrients (primarily nitrogen and phosphorus) regulating carbon–climate feedback are clearly required in ecosystem biogeochemistry and land models (Zaehle and Dalmonech, 2011; Zhang et al., 2011). For example, none of the land models in C${}^{\mathrm{4}}$MIP (Coupled Climate Carbon Cycle Model Intercomparison Project phase 4) had coupled carbon and nitrogen dynamics (Friedlingstein et al., 2006). The current generation of CMIP5 (Anav et al., 2013) models used for the recent IPCC (Intergovernmental Panel on Climate Change) assessment had only two members (CLM4CN: Thornton et al., 2007; and BNU-ESM: Ji et al., 2014) that considered nitrogen regulation of terrestrial carbon dynamics. However, as discussed below, several recent studies have shown that these models had large biases in most of the individual processes important for simulating nutrient dynamics. We therefore believe that, at the global scale, no credible representation of nutrient constraints on terrestrial carbon cycling yet exists in ESMs.

Further, none of the CMIP5 ESMs included a phosphorus cycle, which is likely important for tropical forest carbon budgets (Vitousek and Sanford, 1986). The recent IPCC report highlights the importance of nitrogen and phosphorus availability on land carbon storage, even though the phosphorus limitation effect is uncertain (Stocker et al., 2013). Since the next generation of ESMs participating in the CMIP6 synthesis will continue to focus on the impacts of a changing climate on terrestrial CO${}_{\mathrm{2}}$ and abiotic exchanges with the atmosphere (Provides, 2014), developing ecologically realistic and observationally constrained representations of soil nutrient dynamics and carbon–nutrient interactions in ESMs is critical.

The importance of nutrient limitations in terrestrial ecosystems has been widely demonstrated by nitrogen and phosphorus fertilization experiments (Elser et al., 2007). For instance, plant net primary production (NPP) is enhanced in plots with nutrient addition (LeBauer and Treseder, 2008). Similarly, plant growth can be stimulated due to atmospheric nitrogen deposition (Matson et al., 2002). Boreal forests are strongly limited by nitrogen availability (Vitousek and Howarth, 1991), because low temperatures reduce nitrogen mineralization (Bonan and Cleve, 1992) and N${}_{\mathrm{2}}$ fixation (DeLuca et al., 2002, 2008). In contrast, tropical forests are often phosphorus limited (Vitousek et al., 2010), since tropical soils are old and phosphorus derived from parent material weathering has been depleted through long-term pedogenesis processes (Vitousek and Farrington, 1997; Walker and Syers, 1976). In natural ecosystems without external nutrients inputs (e.g., N deposition), soil nitrogen or phosphorus (or both) are likely insufficient to satisfy both plant and microorganism demands (Vitousek and Farrington, 1997). Plants have to compete with microorganisms and mineral surfaces (Kaye and Hart, 1997; Schimel et al., 1989) to obtain sufficient nutrients to sustain their biological processes (e.g., photosynthesis, respiration). Therefore, it is critical to improve the representation of nutrient competition to accurately model how terrestrial ecosystems will respond to perturbations in soil nutrient dynamics (e.g., from elevated nitrogen deposition or CO${}_{\mathrm{2}}$ fertilization-induced nutrient requirements).

Intense competition between plants and microorganisms is a well-observed phenomenon in nutrient-limited systems (Hodge et al., 2000a; Johnson, 1992; Kaye and Hart, 1997). Previously, plants were thought to be initial losers in nutrient competition, due to the fact that microbes are more intimately associated with substrates (Woodmansee et al., 1981). However, increasing observational evidence indicates that plants compete effectively with soil microorganisms (Schimel and Bennett, 2004) under certain circumstances, sometimes even outcompeting them and suppressing microbial growth (Hu et al., 2001; Wang and Lars, 1997). ${}^{15}$N isotope studies have also demonstrated that plants can capture a large fraction of added nitrogen (Hodge et al., 2000b; Marion et al., 1982). In the short term (days to months), plants maintain their competitiveness mainly through (1) establishing mycorrhizal fungi associations (Drake et al., 2011; Rillig et al., 1998), which help plants acquire organic and inorganic forms of nitrogen (Hobbie and Hobbie, 2006; Hodge and Fitter, 2010) and (2) root exudation of extracellular enzymes that decompose rhizosphere soil organic matter (Phillips et al., 2011). In the relatively longer term (months to years), morphological adjustment occurs; for example, plants allocate more carbon to fine roots to explore laterally and deeper (Iversen et al., 2011; Jackson et al., 2009). Finally, over the course of years to decades, plant succession can occur (Medvigy et al., 2009; Moorcroft et al., 2001) and the new plant demography will need to be considered to represent nutrient controls on this timescale.

Given these patterns from the observational literature, nutrient competition is either absent or over-simplified in existing ESMs. One common representation of plant–microbe competition is that plants compete poorly against microbes in resource acquisition. For example, the O-CN land model (Zaehle and Friend, 2010) assumes that soil-decomposing microbes have the priority to immobilize soil mineral nitrogen. After microbes meet their demands, the remaining nitrogen is then available for plant uptake.

Model structure. Boxes represent pools, solid arrows represent aqueous fluxes and dashed arrows represent gaseous pathways out or into the system. Three essential chemical elements (carbon, C; nitrogen, N; and phosphorus, P) are simulated in N-COM (blue, red and green represent C, N and P pools and processes, respectively).

[Figure omitted. See PDF]

Another treatment in ESM land models is that microbial and plant nutrient acquisition competitiveness is based on their relative demands. For example, CLM4CN (Thornton et al., 2007) assumes that the plant and microbial nitrogen demands are satisfied simultaneously. Under nitrogen infertile conditions, all nitrogen demands in the system are down-regulated proportional to the individual demands and subject to available soil mineral nitrogen. This approach led to unrealistic diurnal cycles of gross primary production (GPP), with midday depressions in GPP occurring because of predicted diurnal depletion of the soil mineral nitrogen pool. Emergent impacts of this conceptualization of nutrient constraints on GPP resulted in poor predictions compared to observations, with smaller than observed plant C growth responses to N deposition (Thomas et al., 2013a) and larger than observed responses to N fertilization (Thomas et al., 2013b). Further, most biogeochemistry models not integrated in ESMs also adopt one of these approaches. For instance, Biome-BGC (Running and Coughlan, 1988), CENTURY (Parton et al., 1988), CASA (Carnegie–Ames–Stanford Approach; Potter et al., 1993) and the Terrestrial Ecosystem Model – TEM (McGuire et al., 1992) assume that available nutrients preferentially satisfy the soil microbial immobilization demand.

We believe the two conceptualizations of competition used in ESMs substantially over-simplify competitive interactions between plants and microbes and lead to biases in carbon cycle predictions. To begin to address the problems with these simplified approaches, Tang and Riley (2013) showed that complex consumer–substrate networks can be represented with an approach – called equilibrium chemical approximation (ECA) kinetics – that simultaneously resolves multiple demands for multiple substrates, and demonstrated that the approach was consistent with observed litter decomposition observations. ECA kinetics has also recently been applied to analyze the emergent temperature response of soil organic matter (SOM) decomposition, considering equilibrium, non-equilibrium and enzyme temperature sensitivities and abiotic interactions with mineral surfaces (Tang and Riley, 2014). We extend on that work here by presenting an implementation of ECA kinetics to represent competition for multiple soil nutrients in a multiple consumer environment. We note that this paper demonstrates a method to handle instantaneous competition in the complex soil–plant network, but a robust competition representation for climate-scale models will require representation of dynamic changes in plant allocation and plant composition.

The aim of this study is to provide a reliable nutrient competition approach applicable for land models integrated in ESMs. However, before integration into an ESM, the competition model needs to be carefully calibrated and independently tested against observational data. This paper will therefore focus on model development and evaluation at several tropical forest sites where observations are available. Our objectives are to (1) develop a soil biogeochemistry model with multiple nutrients (i.e., NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$ and PO${}_x$; represented as the sum of PO${}_{\mathrm{4}}^{\mathrm{3}-}$, HPO${}_{\mathrm{4}}^{\mathrm{2}-}$ and H${}_{\mathrm{2}}$PO${}_{\mathrm{4}}^{-})$ and multiple nutrient consumers (i.e., decomposing microbes, plants, nitrifiers, denitrifiers, and mineral surfaces) competition using ECA kinetics (Tang and Riley, 2013; Zhu and Riley, 2015); (2) constrain the model with in situ observational data sets of soil carbon, nitrogen, and phosphorus dynamics using a Markov chain Monte Carlo (MCMC) approach; and (3) test model performance against nitrogen and phosphorus fertilization studies.

Method

Model development

The Nutrient COMpetition model (N-COM) is designed as a soil biogeochemistry model (Fig. 1) to simulate soil carbon decomposition, nitrogen and phosphorus transformations, abiotic interactions and plant demands. Although our ultimate goal is to incorporate N-COM into a decomposition model that represents active microbial activity as the primary driver of decomposition, we start here by presenting the N-COM approach using a Century-like (Koven et al., 2013; Parton et al., 1988) structure, with additions to account for phosphorus dynamics. In our approach, we calculate potential immobilization using literature-derived parameters (e.g., VMAX, $K_M)$ in a Michaelis–Menten (MM) kinetics framework. The potential immobilization is subsequently modified using the ECA competition method.

Five pools of soil organic carbon (C), nitrogen (N) and phosphorus (P) are considered: coarse wood debris (CWD), litter, fast soil organic matter (SOM) pool, medium SOM pool, and slow SOM pool. Litter is further divided into three sub-groups: metabolic, cellulose and lignin. The soil organic C, N and P decomposition ($F_{\mathrm{C},j}^{\mathrm{dec}}$, $F_{\mathrm{N},j}^{\mathrm{dec}}$, $F_{\mathrm{P},j}^{\mathrm{dec}})$ follow first-order decay: $$\begin{array}{}& & F_{\mathrm{C},j}^{\mathrm{dec}}=k_j{\mathrm{C}}_j{r}_{\mathit{\theta }}{r}_T& & F_{\mathrm{N},j}^{\mathrm{dec}}=k_j{\mathrm{N}}_j{r}_{\mathit{\theta }}{r}_T& & F_{\mathrm{P},j}^{\mathrm{dec}}=k_j{\mathrm{P}}_j{r}_{\mathit{\theta }}{r}_T,\end{array}$$ where $k_j$ is the rate constant of soil organic matter decay (s${}^{-\mathrm{1}})$; C${}_j$, N${}_j$ and P${}_j$ are pool sizes (g m${}^{-\mathrm{2}})$ of carbon, nitrogen and phosphorus, respectively ($j$ from 1 to 7 represents the soil organic matter pools: CWD, metabolic litter, cellulose litter, lignin litter, fast soil organic carbon (SOC), median SOC, slow SOC); $r_T$ and $r_{\mathit{\theta }}$ (dimensionless) are soil temperature and moisture environmental regulators.

Decomposed carbon ($F_{C,i}^{\mathrm{dec}}$; upstream $i$th pool) either (1) enters a downstream pool ($j$th) or (2) is lost as CO${}_{\mathrm{2}}$. Soil organic carbon (downstream $j$th pool) temporal change is calculated as $${\displaystyle \frac{{\mathrm{dC}}_j}{\mathrm{d}t}}=-F_{\mathrm{C},j}^{\mathrm{dec}}+\sum _{i=\mathrm{1}}^{\mathrm{N}}{F}_{\mathrm{C},ij}^{\text{move}},$$ where $\sum _{i=\mathrm{1}}^{\mathrm{N}}{F}_{\mathrm{C},ij}^{\text{move}}$ is the summation of carbon fluxes that move from the upstream pool ($i)$ to the downstream pool ($j)$ due to the decomposition of upstream SOC. For each upstream carbon pool ($i=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{7}$), the fractions integrated into downstream pools ($j=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{7}$) is summarized in a 7 $\times$ 7 matrix $f_{ij}$ (Table 2). The percentage of decomposed carbon that is respired as CO${}_{\mathrm{2}}$ is represented by $g_i$ (Table 2). Simultaneously, soil organic N and P changes follow C decomposition: $$\begin{array}{}& & {\displaystyle \frac{{\mathrm{dN}}_j}{\mathrm{d}t}}=-F_{\mathrm{N},j}^{\mathrm{dec}}+\sum _{i=\mathrm{1}}^{\mathrm{N}}{F}_{\mathrm{N},ij}^{\text{move}}+\sum _{i=\mathrm{1}}^{\mathrm{N}}{F}_{{\mathrm{NH}}_{\mathrm{4}},ij}^{\text{immob}}+\sum _{i=\mathrm{1}}^{\mathrm{N}}{F}_{{\mathrm{NO}}_{\mathrm{3}},ij}^{\text{immob}}& & {\displaystyle \frac{{\mathrm{dP}}_j}{\mathrm{d}t}}=-F_{\mathrm{P},j}^{\mathrm{dec}}+\sum _{i=\mathrm{1}}^{\mathrm{N}}{F}_{\mathrm{P},ij}^{\text{move}}+\sum _{i=\mathrm{1}}^{\mathrm{N}}{F}_{\mathrm{P},ij}^{\text{immob}},\end{array}$$ where $F_{\mathrm{N},ij}^{\text{move}}$ and $F_{\mathrm{P},ij}^{\text{move}}$ are fluxes of nitrogen and phosphorus moving from the upstream ($i)$ to downstream ($j)$ pools. $F_{{\mathrm{NH}}_{\mathrm{4}},ij}^{\text{immob}}$, $F_{{\mathrm{NO}}_{\mathrm{3}},ij}^{\text{immob}}$ and $F_{\mathrm{P},ij}^{\text{immob}}$ are immobilization fluxes of soil mineral nitrogen and phosphorus. $F_{\mathrm{N},j}^{\mathrm{dec}}$ and $F_{\mathrm{P},j}^{\mathrm{dec}}$ represent soil organic matter decomposition losses.

Equations (5) and (6) state that changes in the $j$th organic N or P pool are the summation of three terms: (1) organic N and P lost during soil organic matter mineralization ($-F_{\mathrm{N},j}^{\mathrm{dec}}$ and $-F_{\mathrm{P},j}^{\mathrm{dec}})$; (2) a fraction of the $i$th organic N or P pool (upstream) enters into the $j$th pool (downstream; $F_{\mathrm{N},ij}^{\text{move}}$ and $F_{\mathrm{P},ij}^{\text{move}})$; and (3) soil microbial immobilization ($F_{{\mathrm{NH}}_{\mathrm{4}},ij}^{\text{immob}}$, $F_{{\mathrm{NO}}_{\mathrm{3}},ij}^{\text{immob}}$ and $F_{\mathrm{P},ij}^{\text{immob}})$. Immobilization occurs only when the newly entering organic N is insufficient to sustain the soil C : N (or C : P) ratio (more details described in Appendix A).

The inorganic nitrogen pools (NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$; Eqs. 7–8) are altered by production (organic N mobilized by microbes), consumption (uptake by plants and microbes, gaseous or aqueous losses) and transformation (nitrification and denitrification). Inorganic P (PO${}_x)$ is assumed to be either taken up by plants and decomposing microbes or adsorbed to mineral surfaces (Eq. 9). Plants utilize all forms of phosphate (e.g., PO${}_{\mathrm{4}}^{\mathrm{3}-}$, HPO${}_{\mathrm{4}}^{\mathrm{2}-}$, and H${}_{\mathrm{2}}$PO${}_{\mathrm{4}}^{-})$, but for simplicity we use the symbol PO${}_x$ to represent the sum of all possible phosphate forms throughout the paper: $$\begin{array}{cc}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{\mathrm{d}t}}=\sum _{j=\mathrm{1}}^{\mathrm{N}}\sum _{i=\mathrm{1}}^{\mathrm{N}}{F}_{{\mathrm{NH}}_{\mathrm{4}},ij}^{\text{mob}}-F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{nit}}-F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}\\ & & -F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{immob}}+F^{\text{BNF}}+F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{dep}}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{\mathrm{d}t}}=-F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{den}}+(\mathrm{1}-f^{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}})F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{nit}}-F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{plant}}-F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{immob}}\\ & & -F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{leach}}+F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{dep}}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{PO}}_x\right]}{\mathrm{d}t}}=\sum _{j=\mathrm{1}}^{\mathrm{N}}\sum _{i=\mathrm{1}}^{\mathrm{N}}{F}_{\mathrm{P},ij}^{\text{mob}}-F_{\mathrm{P}}^{\text{plant}}-F_{\mathrm{P}}^{\text{immob}}\\ & & -F_{\mathrm{P}}^{\text{surf}}-F_{\mathrm{P}}^{\text{leach}}+F^{\text{weather}},\end{array}$$ where $F_{{\mathrm{NH}}_{\mathrm{4}},ij}^{\text{mob}}$ and $F_{\mathrm{P},ij}^{\text{mob}}$ are gross mineralization rates for nitrogen and phosphorus. $F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{nit}}$ is the nitrification flux, part of which is lost through a gaseous pathway ($f^{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}})$ and the rest is incorporated into the NO${}_{\mathrm{3}}^{-}$ pool. $F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{den}}$ is the denitrification flux, which transforms nitrate to N${}_{\mathrm{2}}$O and N${}_{\mathrm{2}}$ which then leave the soil system. Plant uptake of soil NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$ and PO${}_x$ are represented as $F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}$, $F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{plant}}$ and $F_{\mathrm{P}}^{\text{plant}}$, respectively. Soil-decomposing microbial immobilization of soil NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$ and PO${}_x$ are represented as $F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{immob}}$, $F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{immob}}$ and $F_{\mathrm{P}}^{\text{immob}}$. The $F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{leach}}$ and $F_{\mathrm{P}}^{\text{leach}}$ are leaching losses of soil NO${}_{\mathrm{3}}^{-}$ and PO${}_x$. External inputs into soil inorganic N pools include atmospheric ammonia deposition ($F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{dep}})$, atmospheric nitrate deposition ($F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{dep}})$ and biological nitrogen fixation ($F^{\text{BNF}})$. External sources of phosphate come from parent material weathering ($F^{\text{weather}})$.

Finally, the dynamics of sorbed P (P${}_{\text{S}})$, occluded P (P${}_{\text{O}})$ and parent material P (P${}_{\text{P}})$ are modeled as $$\begin{array}{}& & {\displaystyle \frac{\mathrm{d}\left[P_{\text{S}}\right]}{\mathrm{d}t}}=F_{\mathrm{P}}^{\text{surf}}-F_{\mathrm{P}}^{\text{occl}}& & {\displaystyle \frac{\mathrm{d}\left[P_{\text{O}}\right]}{\mathrm{d}t}}=F_{\mathrm{P}}^{\text{occl}}& & {\displaystyle \frac{\mathrm{d}\left[P_{\text{P}}\right]}{\mathrm{d}t}}=-F^{\text{weather}}+F_{\mathrm{P}}^{\text{dep}},\end{array}$$ where the pool of sorbed P is balanced by the adsorption flux ($F_{\mathrm{P}}^{\text{surf}})$ and occlusion flux ($F_{\mathrm{P}}^{\text{occl}})$. Parent material is lost by weathering ($F^{\text{weather}})$ and is slowly replenished by external atmospheric phosphorus inputs ($F_{\mathrm{P}}^{\text{dep}}$, such as dust). More detailed information on the modeled C, N and P fluxes is documented in Appendix A.

A summary of the modeled consumer–resource competition network.

Resources	Consumers
NH${}_{\mathrm{4}}^{+}$	Plant	Decomposing microbe	Nitrifier
NO${}_{\mathrm{3}}^{-}$	Plant	Decomposing microbe	Denitrifier
PO${}_x$	Plant	Decomposing microbe	Mineral surface

Multiple-consumer multiple-resource competition network

The soil biogeochemistry model presented in Sect. 2.1 has multiple potential nutrient consumers (plants, SOM decomposing microbes, nitrifiers, denitrifiers, mineral surfaces) and multiple soil nutrients (NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$, PO${}_x)$. The consumer–resource network is summarized in Table 1. As in many land BGC models (CLM, Century, etc.), we have not explicitly included the mineral surface adsorption of NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$, since we assume ammonia is quickly protected by mineral surfaces from leaching (no leaching term in Eq. 7) but then released for plant and microbial uptake when the biotic demand arises. An improved treatment of these dynamics would necessitate a prognostic model for pH, which is beyond the scope of this analysis. Unlike sorbed P (which can be occluded), there is no further abiotic loss of sorbed ammonia. Therefore, the free ammonia pool is interpreted in the current model structure as a potential free ammonia pool (free $+$ sorbed).

Model parameters and baseline values.

C associated
$g_i$	Percentage of carbon remains in the	–	(1.0; 0.45; 0.5; 0.5;	Koven et al. (2013)
	soil after decomposition of $i$th SOM		0.83; 0.45; 0.45)
$f_{ij}$	fraction of SOM leave from $i$th	–	(0, 0, 0.76, 0.24, 0, 0, 0;	Koven et al. (2013)
	pool and enter into $j$th pool		0, 0, 0, 0, 1, 0, 0;
			0, 0, 0, 0, 1, 0, 0;
			0, 0, 0, 0, 0, 1, 0;
			0, 0, 0, 0, 0, 0.995, 0.005;
			0, 0, 0, 0, 0.93, 0, 0.07;
			0, 0, 0, 0, 1, 0, 0)
CN	Soil organic matter CN ratio	–	(13,16,7.9)	Parton et al. (1988)
CP	Soil organic matter CP ratio	–	(110,320,114)	Parton et al. (1988)
TURN${}_{\text{SOM}}$	Soil organic matter turn over	year	(4.1, 0.066, 0.25, 0.25,	Koven et al. (2013)
	(CWD, metabolic lit, cellulose lit,		0.17, 5, 270)
	lignin lit, fast SOM, medium SOM, slow SOM)
N associated
$k_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}$	Reaction rate of plant NH${}_{\mathrm{4}}^{+}$ carrier enzyme	day${}^{-\mathrm{1}}$	120${}^{\mathrm{a}}$	Jackson et al. (1997);
				Min et al. (2000)
$K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}$	Half-saturation constant for plant NH${}_{\mathrm{4}}^{+}$ uptake	g m${}^{-\mathrm{2}}$	0.09	Kuzyakov and Xu (2013)
$K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}$	Half-saturation constant for	g m${}^{-\mathrm{2}}$	0.02	Kuzyakov and Xu (2013)
	decomposing microbe NH${}_{\mathrm{4}}^{+}$ immobilization
$k_{\text{nit}}$	Maximum fraction of NH${}_{\mathrm{4}}^{+}$	day${}^{-\mathrm{1}}$	10 %	Parton et al. (2001)
	pool that could be utilized by nitrifiers
$K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}}$	Half-saturation constant for nitrifier NH${}_{\mathrm{4}}^{+}$ consumption	g m${}^{-\mathrm{2}}$	0.076	Drtil et al. (1993)
$k_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{plant}}$	Reaction rate of plant NO${}_{\mathrm{3}}^{-}$ carrier enzyme	day${}^{-\mathrm{1}}$	2${}^{\mathrm{a}}$	Jackson et al. (1997);
				Min et al. (2000)
$K_M^{\text{plant},{\mathrm{NO}}_{\mathrm{3}}}$	Half-saturation constant for plant NO${}_{\mathrm{3}}^{-}$ uptake	g m${}^{-\mathrm{2}}$	0.07	Kuzyakov and Xu (2013)
$K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}$	Half-saturation constant for decomposing microbe NO${}_{\mathrm{3}}^{-}$	g m${}^{-\mathrm{2}}$	0.04	Kuzyakov and Xu (2013)
	microbe NO${}_{\mathrm{3}}^{-}$ immobilization
$K_M^{\text{den},{\mathrm{NO}}_{\mathrm{3}}}$	Half-saturation constant for denitrifier NO${}_{\mathrm{3}}^{-}$ consumption	g m${}^{-\mathrm{2}}$	0.011	Murray et al. (1989)
$\left[E_{\text{N}}^{\text{plant}}\right]$	Plant nitrogen carrier enzyme abundance for nitrogen uptake	g m${}^{-\mathrm{2}}$	$C_{\text{froot}}\cdot {0.0000125}^{\mathrm{a}}$	Tang and Riley (2013);
				Trumbore et al. (2006)
$\left[E_{\text{N}}^{\text{mic}}\right]$	Decomposing microbes nitrogen carrier	g m${}^{-\mathrm{2}}$	${\frac{F_N^{\text{immob,pot}}}{1000}}^{\mathrm{b}}$	Tang and Riley (2013)
	enzyme abundance for nitrogen immobilization
$\left[E_{\text{N}}^{\text{nit}}\right]$	Nitrifier nitrogen carrier enzyme	g m${}^{-\mathrm{2}}$	1.2E${}^{-\mathrm{3}}$	Raynaud et al. (2006)
	abundance for NH${}_{\mathrm{4}}^{+}$ assimilation
$\left[E_{\text{N}}^{\text{den}}\right]$	Denitrifier nitrogen carrier enzyme	g m${}^{-\mathrm{2}}$	1.2E${}^{-\mathrm{3}}$	Raynaud et al. (2006)
	abundance for NO${}_{\mathrm{3}}^{-}$ assimilation
$f^{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$	Fraction of nitrification flux lost as N${}_{\mathrm{2}}$O	–	6E${}^{-\mathrm{4}}$	Li et al. (2000)
P associated
$k_{\text{weather}}$	Parent material P weathering rate	g P m${}^{-\mathrm{2}}$ year${}^{-\mathrm{1}}$	0.004	Wang et al. (2010)
$k_{\text{occl}}$	P occlude rate	month${}^{-\mathrm{1}}$	1.0E${}^{-\mathrm{6}}$	Yang et al. (2014)
$k_P^{\text{plant}}$	Reaction rate of plant PO${}_x$ carrier enzyme	day${}^{-\mathrm{1}}$	12${}^{\mathrm{a}}$	Colpaert et al. (1999)
$K_M^{\text{plant,P}}$	Half-saturation constant for plant PO${}_x$ uptake	g m${}^{-\mathrm{2}}$	0.067	Cogliatti and Clarkson (1983)
$K_M^{\text{mic,P}}$	Half-saturation constant for	g m${}^{-\mathrm{2}}$	0.02	Chen (1974)
	decomposing microbe PO${}_x$ immobilization
VMAX${}_{\text{P}}^{\text{surf}}$	Maximum mineral surface PO${}_x$ adsorption	g m${}^{-\mathrm{2}}$	133	Wang et al. (2010)
$K_M^{\text{surf,P}}$	Half-saturation constant for mineral surface PO${}_x$ adsorption	g m${}^{-\mathrm{2}}$	64	Wang et al. (2010)
$\left[E_{\text{P}}^{\text{plant}}\right]$	Plant phosphorus carrier enzyme abundance for PO${}_x$ uptake	g m${}^{-\mathrm{2}}$	$C_{\text{froot}}\cdot {0.0000125}^{\mathrm{a}}$	Tang and Riley (2013);
				Trumbore et al. (2006)
$\left[E_{\text{P}}^{\text{mic}}\right]$	Decomposing microbes phosphorus	g m${}^{-\mathrm{2}}$	${\frac{F_{\mathrm{P}}^{\text{immob,pot}}}{400}}^{\mathrm{b}}$	Tang and Riley (2013)
	carrier enzyme abundance for PO${}_x$ immobilization
$\left[E_{\text{P}}^{\text{surf}}\right]$	Mineral surface “effective enzyme” abundance for PO${}_x$ adsorption	g m${}^{-\mathrm{2}}$	VMAX${}_{\text{P}}^{\text{surf}}-\left[PS\right]$	Tang and Riley (2013)

${}^{\mathrm{a}}$ The scaling factor for plant nutrient enzyme abundance is 0.0000125. This number is inferred by assuming that growing season plant nutrient carrier enzymes are roughly the same order of magnitude compared with decomposing microbes'. Typical values for soil decomposing microbe biomass and tropical forest fine root biomass are 0.1 (Tang and Riley, 2013) and 400 gC m${}^{-\mathrm{2}}$ (Trumbore et al., 2006). A typical value of scaling factor that scales microbial biomass to enzyme abundance is 0.05 (Tang and Riley, 2013). Therefore,$C_{\text{froot}}\cdot x=C_{\text{mic}}\cdot 0.05$ or $400\cdot x=0.1\cdot 0.05$. We have $x=0.0000125$. Further, we have $k_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}\cdot \left[E_{\text{N}}^{\text{plant}}\right]={\text{VMAX}}_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}$. We know that typical values for VMAX${}_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}$ and $\left(E_{\text{N}}^{\text{plant}}\right)$ are 0.6 g m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}}$ (Min et al., 2000) and 0.005 g m${}^{-\mathrm{2}}$. Then we have $k_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}=120$ day${}^{-\mathrm{1}}$. Similarly, we have $k_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{plant}}\cdot \left(E_{\text{N}}^{\text{plant}}\right)={\text{VMAX}}_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{plant}}$, $k_P^{\text{plant}}\cdot \left[E_{\text{P}}^{\text{plant}}\right]={\text{VMAX}}_P^{\text{plant}}$. Knowing that typical values for VMAX${}_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{plant}}$ and VMAX${}_{\text{P}}^{\text{plant}}$ are 0.01 (Min et al., 2000) and 0.06 (Colpaert et al., 1999) g m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}}$, we have $k_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{plant}}=\mathrm{2}$ and $k_P^{\text{plant}}=12$ day${}^{-\mathrm{1}}$. ${}^{\mathrm{b}}$ For decomposing microbes, we have VMAX${}_N^{\text{mic}}=k_N^{\text{mic}}\cdot \left[E_N^{\text{mic}}\right]$. Typical values for VMAX${}_N^{\text{immob}}$ and $\left[E_{\text{N}}^{\text{mic}}\right]$ are 5 g m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}}$ (Kuzyakov and Xu, 2013) and 0.005 g m${}^{-\mathrm{2}}$ (Tang and Riley, 2013). Therefore, we have $k_N^{\text{mic}}=1000$. Since our model calculates potential N immobilization rates and approximates them as VMAX${}_N^{\text{mic}}$. The changes of potential N immobilization rates at each time step imply the changes of enzyme abundance through $\left[E_{\text{N}}^{\text{mic}}\right]=\frac{F_N^{\text{immob,pot}}}{k_N^{\text{mic}}}=\frac{F_N^{\text{immob,pot}}}{1000}$. Similarly, we have that VMAX${}_{\text{P}}^{\text{immob}}$ and $\left[E_{\text{P}}^{\text{mic}}\right]$ are 2 g m${}^{-\mathrm{2}}$ day${}^{-\mathrm{1}}$ (Chen, 1974) and 0.005 g m${}^{-\mathrm{2}}$. Therefore, $k_P^{\text{mic}}=400$ and $E_{\text{P}}^{\text{mic}}=\frac{F_{\mathrm{P}}^{\text{immob,pot}}}{400}$.

Competition between different consumers in acquiring different resources is summarized in Table 1. Each consumer–substrate competition reaction is represented by $$\mathrm{S}+E\underset{k_{\mathrm{1}}^{-}}{\overset{k_{\mathrm{1}}^{+}}{\rightleftharpoons }}\mathrm{C}\stackrel{k_{\mathrm{2}}^{+}}{⟶}\mathrm{P}+E.$$ The enzyme ($E$: e.g., nutrient carrier enzyme produced by plants and microbes) and substrate (S: e.g., NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-})$ reaction (reversible reaction) forms a substrate-enzyme complex (C). The following irreversible reaction leads to product (P: meaning the nutrients has been taken up) and releases enzyme ($E$) back into soil media. For the whole complex reaction network, nutrient uptakes are formulated as $$\begin{array}{}& & F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}=k_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}\cdot {\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]\cdot \left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}(\mathrm{1}+{\frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}}}^{\left(\mathrm{1}\right)}+{\frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{K_M^{\text{plant},{\mathrm{NO}}_{\mathrm{3}}}}}^{\left(\mathrm{2}\right)}+{\frac{\left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}}}^{\left(\mathrm{3}\right)}+{\frac{\left[E_{\text{N}}^{\text{mic}}\right]}{K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}}}^{\left(\mathrm{4}\right)}+{\frac{\left[E_{\text{N}}^{\text{nit}}\right]}{K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}}}}^{\left(\mathrm{5}\right)})}}& & F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{immob}}=k_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{immob}}\cdot {\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]\cdot \left[E_{\text{N}}^{\text{mic}}\right]}{K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}(\mathrm{1}+\frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[E_{\text{N}}^{\text{mic}}\right]}{K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[E_{\text{N}}^{\text{nit}}\right]}{K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}}})}}& & F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{nit}}=k_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{nit}}\cdot {\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]\cdot \left[E_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{nit}}\right]}{K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}}(\mathrm{1}+\frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left(\right[E_{\text{N}}^{\text{mic}}]}{K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[E_{\text{N}}^{\text{nit}}\right]}{K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}}})}}& & F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{plant}}=k_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{plant}}\cdot {\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]\cdot \left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NO}}_{\mathrm{3}}}(\mathrm{1}+\frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{K_M^{\text{plant},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{mic}}\right]}{K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{den}}\right]}{K_M^{\text{den},{\mathrm{NO}}_{\mathrm{3}}}})}}& & F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{immob}}=k_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{immob}}\cdot {\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]\cdot \left[E_{\text{N}}^{\text{mic}}\right]}{K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}(\mathrm{1}+\frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{mic}}\right]}{K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{den}}\right]}{K_M^{\text{den},{\mathrm{NO}}_{\mathrm{3}}}})}}& & F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{den}}=k_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{den}}\cdot {\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]\cdot \left[E_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{den}}\right]}{K_M^{\text{den},{\mathrm{NO}}_{\mathrm{3}}}(\mathrm{1}+\frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{K_M^{\text{den},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{mic}}\right]}{K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{den}}\right]}{K_M^{\text{den},{\mathrm{NO}}_{\mathrm{3}}}})}}& & F_{\mathrm{P}}^{\text{plant}}=k_P^{\text{plant}}\cdot {\displaystyle \frac{\left[{\mathrm{PO}}_x\right]\cdot \left[E_{\text{P}}^{\text{plant}}\right]}{K_M^{\text{plant,P}}(\mathrm{1}+\frac{\left[{\mathrm{PO}}_x\right]}{K_M^{\text{plant,P}}}+\frac{\left[E_{\text{P}}^{\text{plant}}\right]}{K_M^{\text{plant,P}}}+\frac{\left[E_{\text{P}}^{\text{mic}}\right]}{K_M^{\text{mic,P}}}+\frac{\left[E_{\text{P}}^{\text{surf}}\right]}{K_M^{\text{surf,P}}})}}& & F_{\mathrm{P}}^{\text{mic}}=k_P^{\text{mic}}\cdot {\displaystyle \frac{\left[{\mathrm{PO}}_x\right]\cdot \left[E_{\text{P}}^{\text{mic}}\right]}{K_M^{\text{mic,P}}(\mathrm{1}+\frac{\left[{\mathrm{PO}}_x\right]}{K_M^{\text{mic,P}}}+\frac{\left[E_{\text{P}}^{\text{plant}}\right]}{K_M^{\text{plant,P}}}+\frac{\left[E_{\text{P}}^{\text{mic}}\right]}{K_M^{\text{mic,P}}}+\frac{\left[E_{\text{P}}^{\text{surf}}\right]}{K_M^{\text{surf,P}}})}}& & F_{\mathrm{P}}^{\text{surf}}=k_P^{\text{surf}}\cdot {\displaystyle \frac{\left[{\mathrm{PO}}_x\right]\cdot \left[E_{\text{P}}^{\text{mic}}\right]}{K_M^{\text{surf,P}}(\mathrm{1}+\frac{\left[{\mathrm{PO}}_x\right]}{K_M^{\text{surf,P}}}+\frac{\left[E_{\text{P}}^{\text{plant}}\right]}{K_M^{\text{plant,P}}}+\frac{\left[E_{\text{P}}^{\text{mic}}\right]}{K_M^{\text{mic,P}}}+\frac{\left[E_{\text{P}}^{\text{surf}}\right]}{K_M^{\text{surf,P}}})}},\end{array}$$ where the $F$ represent the nutrient uptake fluxes and $k$ is the base reaction rate that enzyme–substrate complex forms product ($k_{\mathrm{2}}^{+}$ in Eq. 13). [$E$] and $K_M$ denote enzyme abundance and half-saturation constants (substrate–enzyme affinity). Superscripts and subscripts refer to consumers and substrates, respectively. These equations account for the effect of (1) multiple substrates (e.g., NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-})$ sharing one consumer, which inhibits the effective binding between any specific substrate and the consumer (terms ${}^{\left(\mathrm{1}\right)}$ and ${}^{\left(\mathrm{2}\right)}$ in Eq. 14) and (2) multiple consumers (e.g., plants, decomposing microbes and nitrifiers) sharing one substrate (e.g., NH${}_{\mathrm{4}}^{+})$, which lowers the probability of effective binding between any consumer and NH${}_{\mathrm{4}}^{+}$ (terms ${}^{\left(\mathrm{3}\right)}$, ${}^{\left(\mathrm{4}\right)}$, and ${}^{\left(\mathrm{5}\right)}$ in Eq. 14).

For our reaction network (Eqs. 13–22), we make the following four assumptions: (1) Plant roots and decomposing microbes possess two types of nutrient carrier enzymes (nutrient transporters). One is for nitrogen (NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$; $E_{\text{N}}^{\text{plant}},E_{\text{N}}^{\text{mic}})$, and the other is for phosphorus, including different forms of phosphate ($E_{\text{P}}^{\text{plant}},E_{\text{P}}^{\text{mic}})$. (2) Nutrient carrier enzyme abundance is scaled with biomass (fine root or microbial biomass). Scaling factors are 0.0000125 (for plants) and 0.05 (for decomposing microbes; Table 2). (3) Mineral surface “effective enzyme” abundance [$E_{\text{P}}^{\text{surf}}]$ is approximated by the available sorption surface area (VMAX${}_{\text{P}}^{\text{surf}}-$[PS]). (4) Nitrifiers and denitrifiers are not explicitly simulated, therefore we assume that their biomass and associated nutrient transporter abundance are fixed ($E_{\text{N}}^{\text{nit}},E_{\text{N}}^{\text{denit}})$.

For simplicity, we group the “decomposing microbes/nitrifier/denitrifier/mineral surface nutrient carrier enzyme ($E$)” and their “base reaction rate $k$” into one single variable “VMAX” (see Appendix B for full derivation). Furthermore, we defined “potential rates (potential immobilization, nitrification, denitrification, adsorption rates)” and used them as proxies of “VMAX”. Therefore, Eqs. (15), (16), (18), (19), (21), and (22) become $$\begin{array}{}& & F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{immob}}=F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{immob,pot}}\cdot {\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}(\mathrm{1}+\frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[E_{\text{N}}^{\text{mic}}\right]}{K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[E_{\text{N}}^{\text{nit}}\right]}{K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}}})}}& & F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{nit}}=F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{nit,pot}}\cdot {\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}}(\mathrm{1}+\frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[E_{\text{N}}^{\text{mic}}\right]}{K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[E_{\text{N}}^{\text{nit}}\right]}{K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}}})}}& & F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{immob}}=F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{immob,pot}}\cdot {\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}(\mathrm{1}+\frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{mic}}\right]}{K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{den}}\right]}{K_M^{\text{den},{\mathrm{NO}}_{\mathrm{3}}}})}}& & F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{den}}=F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{den,pot}}\cdot {\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{K_M^{\text{den},{\mathrm{NO}}_{\mathrm{3}}}(\mathrm{1}+\frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{K_M^{\text{den},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{mic}}\right]}{K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{den}}\right]}{K_M^{\text{den},{\mathrm{NO}}_{\mathrm{3}}}})}}& & F_{\mathrm{P}}^{\text{mic}}=F_{\mathrm{P}}^{\text{immob,pot}}\cdot {\displaystyle \frac{\left[{\mathrm{PO}}_x\right]}{K_M^{\text{mic,P}}(\mathrm{1}+\frac{\left[{\mathrm{PO}}_x\right]}{K_M^{\text{mic,P}}}+\frac{\left[E_{\text{P}}^{\text{plant}}\right]}{K_M^{\text{plant,P}}}+\frac{\left[E_{\text{P}}^{\text{mic}}\right]}{K_M^{\text{mic,P}}}+\frac{\left[E_{\text{P}}^{\text{surf}}\right]}{K_M^{\text{surf,P}}})}}& & F_{\mathrm{P}}^{\text{surf}}=F_{\mathrm{P}}^{\text{surf,pot}}\cdot {\displaystyle \frac{\left[{\mathrm{PO}}_x\right]}{K_M^{\text{surf,P}}(\mathrm{1}+\frac{\left[{\mathrm{PO}}_x\right]}{K_M^{\text{surf,P}}}+\frac{\left[E_{\text{P}}^{\text{plant}}\right]}{K_M^{\text{plant,P}}}+\frac{\left[E_{\text{P}}^{\text{mic}}\right]}{K_M^{\text{mic,P}}}+\frac{\left[E_{\text{P}}^{\text{surf}}\right]}{K_M^{\text{surf,P}}})}}.\end{array}$$ In this case, the potential rates are treated as maximum reaction rates (VMAX), because they are calculated without nutrient constraints or biotic and abiotic interactions. For example, potential P immobilization rate ($F_{\mathrm{P}}^{\text{immob,pot}})$ is based on the total phosphorus demand that can perfectly maintain the soil CP stoichiometry during soil organic matter decomposition (Eq. A9). This potential immobilization rate represents the maximum phosphorus influx that the soil could take up at that moment. The maximum adsorption rate ($F_{\mathrm{P}}^{\text{surf,pot}})$ is the time derivative of the Langmuir equation (Eq. A12), which is a theoretically maximal adsorption rate excluding all other biotic and abiotic interactions. The potential rates (VMAX) are updated by the model rather than calibrated, except for VMAX${}_{\text{P}}^{\text{surf}}$. VMAX${}_{\text{P}}^{\text{surf}}$ denotes the maximum adsorption capacity (not maximum adsorption rate), which affects the potential adsorption rate ($F_{\mathrm{P}}^{\text{surf,pot}})$.

Observational data sets used for calibration. Number of observations for each data stream is included in brackets.

Processes	Data sets		Location	References
C associated	Soil heterotrophic		Tapajos National Forest,	(Silver et al., 2012)
	respiration (20)		Para, Brazil
N associated	Soil NH${}_{\mathrm{4}}^{+}$ (5)	N${}_{\mathrm{2}}$O efflux (20)	Tapajos National Forest,	(Silver et al., 2012)
			Para, Brazil
P associated	Soil free phosphate (3)	Sorb phosphate (3)	Tapajos National	(McGroddy et al., 2008)
			Forest, Para, Brazil

The model is run on an hourly time step, initialized with state variables and critical parameters (Table 2). Since the model is designed to be a component of the Community and ACME Land Models (CLM, ALM; which are essentially currently equivalent), we used CLM4.5 site-level simulations to acquire temporally-resolved: (1) soil temperature factors on decomposition ($r_T)$; (2) soil moisture factors on decomposition ($r_{\mathit{\theta }})$; (3) the anoxic fraction of soil pores ($f^{\text{anox}}$ in Appendix Eqs. A10–11); (4) annual NPP (NPP${}_{\text{annual}}$ in Appendix Eq. A13); (5) NH${}_{\mathrm{4}}^{+}$ deposition ($F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{dep}})$; (6) NO${}_{\mathrm{3}}^{-}$ deposition ($F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{dep}})$; and (7) hydrologic discharge ($Q_{\text{dis}}$ in Appendix Eq. A14). External inputs of mineral phosphorus are derived from Mahowald et al. (2005, 2008).

Model parameterization and sensitivity analysis

We constrained model parameters and performed sensitivity analyses using a suite of observations distinct from the observations we used subsequently to test the model against the N and P manipulation experiments. Because tropical systems can be either nitrogen or phosphorous limited (or both; Elser et al., 2007; Vitousek et al., 2010), we chose observations from a tropical forest site to constrain the N and P competition in our model (Tapajos National Forest, Para, Brazil; Table 3).

Calibrated parameters are reported in terms of (1) mean/standard deviation by fitting to a Gaussian distribution; (2) 25 and 75 % quantile. Both variance-based and quantile-based parameters uncertainty reduction are provided; (3) Gelman–Rubin convergence criterion.

Parameters	${\mathit{\mu }}_{\text{prior}}$	${\mathit{\sigma }}_{\text{prior}}$	${\mathit{\mu }}_{\text{posterior}}$	${\mathit{\sigma }}_{\text{posterior}}$	UR	$Q_{\text{prior}}^{25}$	$Q_{\text{prior}}^{75}$	$Q_{\text{posterior}}^{25}$	$Q_{\text{posterior}}^{75}$	UR	Gelman–Rubin
											criterion
TURN${}_{\text{SOM}}$	(3.7, 0.06,	(3.9, 0.06,	(5.2, 0.07,	(0.33, 0.01,	(92, 83,	(5.33, 0.086,	(19.32, 0.31,	(5.05, 0.63,	(5.39, 0.076,	(97, 94,	(1.69, 1.03,
(CWD, metabolic, cellulose,	0.23, 0.23,	0.24, 0.24,	0.17, 0.17,	0.01, 0.005,	96, 98,	0.33, 0.33,	1.18, 1.18,	0.16, 0.17,	0.18, 0.18,	97, 99,	1.75, 1.01,
lignin lit, fast, medium SOM)	0.16, 4.6)	0.18, 4.8)	0.14, 3.6)	0.008, 0.37)	96, 92)	0.22, 6.5)	0.8, 23.5)	0.13, 3.2)	0.14, 3.9)	98, 96	1.06, 1.55)
$k_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}$	109	114	58	14	88	156.1	565.4	52.8	60.0	98	1.87
$K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}$	0.082	0.086	0.173	0.018	79	0.12	0.42	0.16	0.18	93	2.86
$K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}$	0.018	0.019	0.071	0.0067	65	0.026	0.094	0.065	0.076	85	1.94
$k_{\text{nit}}$	0.091	0.095	0.37	0.038	60	0.13	0.47	0.36	0.39	91	2.02
$K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}}$	0.069	0.072	0.082	0.012	83	0.10	0.36	0.07	0.09	94	1.95
$k_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{plant}}$	1.8	1.9	7.6	1.7	13	2.60	9.42	6.11	9.14	56	1.01
$K_M^{\text{plant},{\mathrm{NO}}_{\mathrm{3}}}$	0.064	0.067	0.085	0.0064	90	0.09	0.33	0.08	0.09	97	3.17
$K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}$	0.036	0.038	0.096	0.014	63	0.05	0.19	0.09	0.10	92	2.64
$K_M^{\text{den},{\mathrm{NO}}_{\mathrm{3}}}$	0.0101	0.0105	0.022	0.0034	68	0.014	0.052	0.019	0.024	87	1.03
$k_P^{\text{plant}}$	11	11.5	59	0.75	93	15.61	56.54	58.86	59.81	98	1.06
$K_M^{\text{plant,P}}$	0.061	0.064	0.11	0.015	77	0.09	0.32	0.10	0.12	94	1.52
$K_M^{\text{mic,P}}$	0.018	0.019	0.037	0.0047	75	0.026	0.094	0.034	0.039	93	2.86
VMAX${}_{\text{P}}^{\text{surf}}$	121	127	182	30	76	173.0	626.6	156.5	206.3	89	2.25
$K_M^{\text{surf,P}}$	64	58	200	50	18	83.2	301.5	162.6	233.0	68	1.05

In the parameter estimation procedure, several data streams are assimilated into the N-COM model, including measurements of soil NH${}_{\mathrm{4}}^{+}$ concentrations, soil free phosphate concentrations, sorbed phosphate concentrations and N${}_{\mathrm{2}}$O and CO${}_{\mathrm{2}}$ flux measurements. The data sets are summarized in Table 3 and cover a wide range of N and P biogeochemistry dynamics. A set of model parameters is selected for calibration (Table 4), which comprise nutrient competition kinetics parameters ($k$ and $K_M)$ as well as the fast soil carbon turnover time (TURN${}_{\text{SOM}})$. Because we had only a short-term CO${}_{\mathrm{2}}$ respiration flux record, we were unable to calibrate the longer turnover time parameters. However, since we test the calibrated model against short-term fertilization responses, this omission will not affect our evaluation. Longer records from eddy covariance flux towers and ${}^{14}$C soil measurements are required to constrain the longer turnover time pool values.

We employed the MCMC approach (Ricciuto et al., 2008) to assimilate the observations into N-COM. MCMC directly draws samples from a pre-defined parameter space and tries to minimize a pre-defined cost function: $$J=\left(\mathit{M}\mathbf{\right(}\mathit{\theta }\mathbf{)}-\mathit{D}{)}^T{R}^{-\mathrm{1}}(M\left(\mathit{\theta }\right)-D),$$ where $\mathit{M}\mathbf{\left(}\mathit{\theta }\mathbf{\right)}$ and $\mathit{D}$ are vectors of model outputs and observations including time series of different simulated variables (e.g., soil CO${}_{\mathrm{2}}$ and N${}_{\mathrm{2}}$O effluxes and soil concentrations of NH${}_{\mathrm{4}}^{+}$, free PO${}_x$ and sorbed PO${}_x)$; $\mathit{\theta }$ is a vector of model parameters (${\mathit{\theta }}_i)$; and $i$ from 1 to 20 represents the parameters that are calibrated (Table 4). $R^{-\mathrm{1}}$ is the inverse of data error covariance matrix. We assumed that diagonal elements are 40 % of observed values and off-diagonal elements are zeros. We further assumed that the prior parameter follows a lognormal distribution. $\mathit{\mu }$ and $\mathit{\sigma }$ were 0.91 and 0.95 of their initial values, respectively (Table 4). We then ran MCMC to sample 50 000 parameter pairs (Fig. A1). The second half of the samples was fit to a Gaussian distribution. We also employed the Gelman–Rubin criterion to quantitatively show whether or not the MCMC chain converged. The calibrated model parameters are reported in terms of means and standard deviations. Uncertainty reduction (UR) is calculated based on (1) variance (Eq. 30a) and (30b) 25 and 75 % quantile (Eq. 30b): $$\begin{array}{}& & {\text{UR}}_{\mathit{\sigma }}=(\mathrm{1}-{\displaystyle \frac{{\mathit{\sigma }}_{\text{posterior}}}{{\mathit{\sigma }}_{\text{prior}}}})\cdot 100\mathit{\%}& & {\text{UR}}_Q=(\mathrm{1}-{\displaystyle \frac{Q_{\text{posterior}}^{75}-Q_{\text{posterior}}^{25}}{Q_{\text{prior}}^{75}-Q_{\text{prior}}^{25}}})\cdot 100\mathit{\%},\end{array}$$ where ${\mathit{\sigma }}_{\text{prior}}$ is prior parameter uncertainty, which is 95 % of the parameter initial value. The ${\mathit{\sigma }}_{\text{posterior}}$ is calibrated parameter uncertainty, which is calculated by fitting the calibrated model parameters to a Gaussian distribution. $Q^{75}$ and $Q^{25}$ are 75 and 25 % percentage quantile of each parameter. UR is a useful metric (Zhu and Zhuang, 2014), because it quantitatively reveals the reduction in the range of a particular parameter after calibration with MCMC. It does not, however, indicate that the parameter itself is more consistent with observed values of the parameter. A large value of UR implies a more robust model.

In addition, we conducted a sensitivity study to identify the dominant controlling factors regulating nutrient competition in N-COM. Three scenarios were considered: (1) baseline climate and soil conditions; (2) elevated soil temperature (by 5 ${}^{\circ }$C); and (3) elevated soil moisture (by 50 %). SOBOL sampling (Pappas et al., 2013), a global sensitivity technique, is employed to calculate the sensitivities of output variables with respect to various inputs:

$$S_i={\displaystyle \frac{{\text{VAR}}_{p_i}\left(E_{p_{\sim i}}\right(Y\left|p_i\right))}{\text{VAR}\left(Y\right)}},$$ where $S_i$ is the first-order sensitivity index of the $i$th parameter and ranges from 0 to 1. By comparing the values of $S_i$, we were able to evaluate which processes affect the pattern of nutrient competition. $Y$ represents the model outputs of plant NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$ or PO${}_x$ uptake; $p_i$ is the target parameter; $p_{\sim i}$ denotes all parameters that are associated with nutrient competition except the target parameter; and VAR(.) and $E(.)$ represent variance and mean, respectively.

Short-term (24 or 48 h) fertilization experiments off NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$ or PO${}_{\mathrm{4}}^{\mathrm{3}-}$ additions used to evaluate the performance of the N-COM competition scheme.

Data sets	Added	Competitors			Duration	References
	nutrient				(hour)
PO${}_{\mathrm{4}}^{\mathrm{3}-}$ fertilization	10 $\mathrm{\mu }$g g${}^{-\mathrm{1}}$	I. Mineral	II. Decomposing		48	Olander and Vitousek (2005)
		surface	microbe
NH${}_{\mathrm{4}}^{+}$ fertilization	4.6 $\mathrm{\mu }$g g${}^{-\mathrm{1}}$	I. Plant	II. Decomposing	III. Nitrifier	24	Templer et al. (2008)
			microbe
NO${}_{\mathrm{3}}^{-}$ fertilization	0.92 $\mathrm{\mu }$g g${}^{-\mathrm{1}}$	I. Plant	II. Decomposing		24	Templer et al. (2008)
			microbe

Model application

After calibration, we applied the N-COM model to several tropical forest nutrient fertilization studies not included in the calibration data set, where isotopically labeled nitrogen or phosphorous fertilizer was injected into the soil. The fertilization experiments measured the fate of added nutrients; for example, identifying the fraction of added N or P that goes into the plant, is immobilized by microbes, or is stabilized by mineral surfaces. These measurements offer an effective baseline to test whether the N-COM model captures short-term nutrient competition.

Because we have focused in this paper on applications in tropical forests, we choose three tropical forest fertilization experiments with (1) PO${}_{\mathrm{4}}^{\mathrm{3}-}$, (2) NH${}_{\mathrm{4}}^{+}$ and (3) NO${}_{\mathrm{3}}^{-}$ additions (Table 5). The PO${}_{\mathrm{4}}^{\mathrm{3}-}$ fertilization experiment (Olander and Vitousek, 2005) was conducted in three Hawaiian tropical forests along a soil chronosequence (300, 20 000 and 4 100 000 year old soils) that were fertilized with 10 $\mathrm{\mu }$g g${}^{-\mathrm{1}}$ ${}^{32}$PO${}_{\mathrm{4}}^{\mathrm{3}-}$, respectively, and microbial demand vs. soil sorption was measured. We did not evaluate the role of plants in phosphorus competition for the Hawaii sites, since plant phosphorus uptake was not measured in those field studies. Our model discriminates the Hawaii sites along the chronosequence by setting distinct initial pool sizes (derived from Olander and Vitousek, 2004, 2005) of soil organic carbon, nitrogen and phosphorus, and soil parent material phosphorus.

Distribution of prior and calibrated model parameters.

[Figure omitted. See PDF]

We also used measurements from NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ fertilization studies located at the Luquillo tropical forest in Puerto Rico (Templer et al., 2008). In that study, 4.6 $\mathrm{\mu }$g g${}^{-\mathrm{1}}$ ${}^{15}$NH${}_{\mathrm{4}}^{+}$ was added into the highly weathered tropical forest soil and the consumption of ${}^{15}$NH${}_{\mathrm{4}}^{+}$ by plant roots, decomposing microbes and nitrifiers were measured. In the same study, 0.92 $\mathrm{\mu }$g g${}^{-\mathrm{1}}$ ${}^{15}$NO${}_{\mathrm{3}}^{-}$ was added to the soil and the plant uptake and microbial immobilization was measured. The measurements were made 24 or 48 h after the fertilizers were added.

For the model scenarios, we (1) spun up the N-COM model for 100 years; (2) perturbed the soil nutrient pool by the same amount as the fertilization; (3) ran the model for 24 or 48 h and calculated how much of the added nutrients were absorbed by plants, microbes or mineral surfaces; and (4) compared our model simulations with the observed data to assess model predictability. The 100-year spin-up simulation aimed at eliminating the effects of imposed initial inorganic pool sizes on fertilization experiments, rather than accumulating soil organic matter in the system, since we initialized the soil organic carbon pools from CLM4.5 steady-state predictions.

Results and discussion

Calibrated model parameters

Our best estimates (second half of the MCMC chain) of the selected model parameters based on the observations at the Tapajos National Forest, Para, Brazil are shown in Fig. 2. We found that calibrated parameter samples were not heavily tailed and they generally follow Gaussian distributions (Fig. A3). In order to quantitatively compare the calibrated parameter distributions with prior distributions, we fit parameter samples to a Gaussian distribution and estimated its means and standard deviations (Table 4).

Even though the parameter mean was improved, the uncertainty may still be relatively large. In other words, a prognostic prediction based on these calibrated parameters could be relatively uncertain (Scholze et al., 2007), due to large uncertainty associated with the calibrated parameters. Therefore, we calculated the variance-based UR (UR${}_{\mathit{\sigma }}$; Eq. 30a) to evaluate model improvement in terms of parameter uncertainty. We found that parameters' uncertainties were reduced by 13–98 %. This calculation might either overestimate or underestimate the UR${}_{\mathit{\sigma }}$, due to the fact that the calibrated parameters did not strictly follow Gaussian distributions. But the actual UR${}_{\mathit{\sigma }}$ should not be far from our estimates, because these samples were not widely spread across the potential parameter space (Fig. 2). The least constrained parameter was $k_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{plant}}$ (reaction rate of plant nitrogen carrier enzyme with NO${}_{\mathrm{3}}^{-}$ substrate). Two other NO${}_{\mathrm{3}}^{-}$ dynamics related parameters were also not well constrained: UR${}_{\mathit{\sigma }}$ of $K_M^{\text{mic},{\mathrm{NO}}_{\mathrm{3}}}$ (half-saturation constant for decomposing microbe NO${}_{\mathrm{3}}^{-}$ immobilization) and $K_M^{\text{den},{\mathrm{NO}}_{\mathrm{3}}}$ (half-saturation constant for denitrifier NO${}_{\mathrm{3}}^{-}$ consumption) were only 63 and 68 %, respectively. Compared with NH${}_{\mathrm{4}}^{+}$ or PO${}_x$ competition-related parameters, we concluded that parameters associated with NO${}_{\mathrm{3}}^{-}$ competition were the least constrained in the model. This result was primarily due to the lack of NO${}_{\mathrm{3}}^{-}$ pool size data, and secondarily due to the fact that NO${}_{\mathrm{3}}^{-}$ was not the major nitrogen source for plant or decomposing microbes. We also provide quantile-based UR for reference (Table 4). The above-mentioned conclusions still hold with quantile-based UR${}_Q$, although the quantile-based UR${}_Q$ is generally higher than variance-based UR${}_{\mathit{\sigma }}$. One parameter was calibrated to be at the upper boundary of its prior ranges ($k_P^{\text{plant}})$, implying that this tropical plant is highly efficient in phosphorus uptake. Although we do not have direct kinetic parameter observations for the specific tropical species involved in our study, an inferred high phosphorus uptake efficiency is reasonable for tropical species that have adapted to these phosphorus-deficient environments (Begum and Islam, 2005; Föhse et al., 1988).

Convergence of model parameters is reported with the Gelman–Rubin criterion (univariate potential scale reduction factor; Table 4 and Fig. A2). Using this criterion, seven (out of twenty) parameters are found to converge (Gelman–Rubin $\le 1.1$). The lack of convergence (in addition, 20-dimensional multivariate potential scale reduction factor is 12.04) of the remaining parameters is partly due to data paucity. In particular, starting from different initial values, MCMC calibrations may result in different models that give rise to similar model–data misfit (i.e., “equifinality”; Tang and Zhuang, 2008). In this regard, high-frequency measurements may improve model calibration (see more discussion in Sect. 3.3). The non-convergence of model parameters implies an imperfect model. Therefore, for large-scale model application, more work on data collection, parameter tuning and uncertainty analysis is needed. However, even with these caveats, the model predictability is reasonably good when applied to the tropical forest fertilization experiments described in Sect. 3.4.

We re-organize the right-hand sides of Eqs. (14)–(22) to be the product of potential nutrient uptake rate and an ECA limitation term; for example for plant NH${}_{\mathrm{4}}^{+}$ uptake, $$\begin{array}{cc}& & F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}=k_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}\cdot {\text{ECA}}_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}& {\text{ECA}}_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}}\\ & & ={\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]\cdot \left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}(\mathrm{1}+\frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{K_M^{\text{plant},{\mathrm{NO}}_{\mathrm{3}}}}+\frac{\left[E_{\text{N}}^{\text{plant}}\right]}{K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[E_{\text{N}}^{\text{mic}}\right]}{K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}}+\frac{\left[E_{\text{N}}^{\text{nit}}\right]}{K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}}})}}.\end{array}$$ Other “consumer-substrate reactions” have similar forms. Under a nutrient-abundant situation (e.g., fertilized agriculture ecosystem), the relative competitiveness of each consumer (ECA) is dominated by its specific enzyme abundance [$E$]. Under such conditions, substrate affinity is no longer a controlling factor. In contrast, under nutrient-limited conditions (e.g., many natural ecosystems), ECA is dominated by the specific enzyme abundance as well as the substrate affinity ([$E$]$/K_M)$. Therefore, consumers could either enable an alternative high-affinity nutrient transporter system (low $K_M)$ or exude more enzyme to enhance competitiveness. For example, at the whole-soil scale it has been shown that root spatial occupation ($C_{\text{froot}})$ determines a plant's competitiveness when low soil nutrient diffusivity is limiting nutrient supply (Raynaud and Leadley, 2004). Consistently, our results highlighted the dominant role of nutrient carrier enzyme abundance ($E$ proportional to $C_{\text{froot}})$ in controlling competition. If we further assumed that plants, decomposing microbes and nitrifiers enzyme abundances were approximately equal, we will have that their relative competitiveness in acquiring NH${}_{\mathrm{4}}^{+}$ was about $\mathrm{4}:10:\mathrm{9}$ ($\mathrm{1}/K_M^{\text{plant},{\mathrm{NH}}_{\mathrm{4}}}:\mathrm{1}/K_M^{\text{mic},{\mathrm{NH}}_{\mathrm{4}}}:\mathrm{1}/K_M^{\text{nit},{\mathrm{NH}}_{\mathrm{4}}})$. However, such results could not be easily generalized to other ecosystems, because they heavily relied on the traits (affinity) of specific competitors. For a different ecosystem, those traits would be drastically different due to the change of, e.g., plant species composition and microbial community structure. Even for the same ecosystem, those traits could be highly heterogeneous. For example, the community structure of decomposing microbes could be different in rhizosphere and bulk soil (with different $K_M)$. However, in this work we assumed a well-mixed environment (one soil column), in order to be consistent with large-scale ecosystem models. Although beyond the scope of the current study, the consequences of ignoring the rhizosphere vs. bulk soil heterogeneity warrants further investigation. Large-scale models aim to quantify ecosystem level dynamics, although they are usually driven by parameters inferred from in situ field observations. In the absence of a model that explicitly represents this spatial heterogeneity, it is difficult to quantify the impacts of using inferred rhizosphere decomposer affinities on model predictions of the whole soil (Schimel et al., 1989). Furthermore, the assumption of well-mixed environment in large-scale model is an inevitable flaw, because of large computational demands and a lack of scale-aware parameters and model structures for large-scale models to run fine-scale simulations.

Model sensitivity analysis with SOBOL sampling. For each metric, three scenarios are shown: baseline (Control), elevated soil temperature by 5 ${}^{\circ }$C ($+T_s)$, and elevated soil moisture by 50 % ($+\mathit{\theta })$, respectively. The length of bar (plot in polar coordinate) is the sensitivity (unit-less) of model output with respect to model input variables. Our results showed that the plant nutrient uptake was mostly regulated by internal consumer–substrate uptake kinetics rather than the external environmental conditions (e.g., $T_{\text{s}}$, $\mathit{\theta })$.

[Figure omitted. See PDF]

Although in this study ECA was applied to a large-scale model, the competition framework is readily applicable to fine-scale models that consider soil heterogeneity. In fine-scale models, bulk soil nutrient competition can occur only among different microbes because they are ubiquitous in the soil (e.g., nitrifier vs. microbial decomposer), while rhizosphere nutrient competition occurs among plants and microbes (e.g., nitrifier vs. microbial decomposer vs. roots). This distinction implies that the competitiveness parameters we infer here for N-COM, which does not currently explicitly represent bulk versus rhizosphere processes, subsume the range of fine-scale processes controlling nutrient uptake. More research is required to link these different model spatial scales, theory and parameterizations.

Our modeling framework highlights the important concept that “competitiveness” is a dynamic property of the competition network, and more importantly that it is linked to competitor functional traits (affinity and nutrient carrier enzyme abundance). This concept is in contrast to the prevailing assumption underlying all major large-scale ecosystem models, which either assume “relative demand competitiveness for different nutrient consumers” (Thornton et al., 2007) or “soil microbes outcompete plants” (McGuire et al., 1992; Parton et al., 1988). Imposing such pre-defined orders of competitiveness neglects the diversity of nutrient competitors (plants and microbes) and their differences in nutrient uptake capacity expressed by relevant functional traits. Our model framework offers a theoretically consistent approach to account for the diversity of nutrient competition in different competitor networks.

Model performance at Tapajos National Forest, Para, Brazil. Overall, the calibrated model (blue line) improved predictions over the prior model (gray line) when compared to observations. Green areas indicate the calibrated model uncertainties.

[Figure omitted. See PDF]

Model sensitivity analysis

Through sensitivity analysis, we separately investigated the factors controlling plant NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$ and PO${}_x$ competition (Fig. 3). Each sensitivity analysis consisted of three scenarios: (1) normal conditions (control); (2) elevated soil temperature ($+T_s)$; and (3) elevated soil moisture ($+\mathit{\theta })$. The sensitivity analysis indicates that the model is highly sensitive to kinetics parameters (e.g., $K_M)$. Furthermore, the model is consistently sensitive to the same parameters across all temperature and moisture conditions. The environment affects the nutrient competition primarily through altering the nutrient abundance. Enhanced soil temperature and soil moisture accelerated soil organic carbon turnover, thereby releasing more inorganic nutrient into the soil (gross mineralization). However, the impacts on plant nutrient uptake are limited (Fig. 3) because the enhanced soil organic matter decay also requires higher immobilization fluxes to sustain the soil organic matter CNP stoichiometry. The enhancement of net mineralization would be limited, and therefore would not change soil nutrient status dramatically.

Model performance

The prior and calibrated models were compared against observational data sets of pool sizes of soil free phosphate, sorbed phosphate, and NH${}_{\mathrm{4}}^{+}$, CO${}_{\mathrm{2}}$ efflux and N${}_{\mathrm{2}}$O efflux (Fig. 4). We note that although we attempted to acquire as many data sets that contained these five observations as possible, more observations in tropical ecosystems would clearly improve the parameter estimates. For example, in the experiment we analyzed, only three measurements of soil free phosphate were made during 1999. Many detailed dynamics are therefore missing and could impact our parameter estimates. The prior model predicted an increasing trend of soil free PO${}_x$, which resulted from underestimates of plant P uptake (by underestimating of $k_{\text{P}}^{\text{plant}})$ and soil microbial P immobilization (by overestimating $K_M^{\text{mic,P}})$. The calibrated model captured the seasonal dynamics of soil free PO${}_x$ reasonably well: increases during the wet season and gradual decreasing during the dry season (August to November). The prior model also largely underestimated the seasonal variability of nitrogen dynamics and underestimated the NH${}_{\mathrm{4}}^{+}$ pool size due to overestimation of plant NH${}_{\mathrm{4}}^{+}$ uptake ($k_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{plant}})$. In addition, it also underestimated the denitrification N${}_{\mathrm{2}}$O emissions, because of an underestimation of NH${}_{\mathrm{4}}^{+}$ to NO${}_{\mathrm{3}}^{-}$ transformation rate ($k_{\text{nit}})$. Consequently, there was not enough NO${}_{\mathrm{3}}^{-}$ substrate to react with denitrifiers and release N${}_{\mathrm{2}}$O. The calibrated model, however, accurately reproduced the seasonal dynamics of both NH${}_{\mathrm{4}}^{+}$ pool sizes and soil N${}_{\mathrm{2}}$O emissions. There were small differences between the prior and calibrated model predictions of soil CO${}_{\mathrm{2}}$ emissions. The CO${}_{\mathrm{2}}$ and N${}_{\mathrm{2}}$O effluxes were more frequently observed at Tapajos National Forest during 1999 to 2001, compared with phosphorus data. Most of the measurements were collected during the wet season. Therefore the modeled CO${}_{\mathrm{2}}$ and N${}_{\mathrm{2}}$O emissions were largely improved by assimilating these data sets.

The model performance implies that after assimilating multiple data sets, our model predictions were improved over the prior model. However, it is clear that more observations of the metrics applied in our MCMC approach would benefit the model calibration. Unfortunately, because of our focus on tropical sites, we were unable to acquire more data sets that had the full suite of measurements required. Data sets of soil nutrient pool sizes (e.g., NO${}_{\mathrm{3}}^{-})$ and higher frequency sampling of those sparse measurements (e.g., PO${}_x$) would significantly benefit the model uncertainty reduction.

Model perturbation experiments compared with nitrogen and phosphorus fertilization field experimental data. The blue dots show the difference between control and perturbed simulations, which mean how much newly added nutrient each consumer takes up. The red circles are recovered isotopically labeled nutrient within each consumer. Since plants phosphorus uptake was not measured at Hawaii sites, we did not include the plants in the perturbation study.

[Figure omitted. See PDF]

Model testing against nitrogen and phosphorus fertilization studies

To test the calibrated N-COM model, we conducted short-term numerical competition experiments (24 or 48 h simulations) by manually imposing an input flux into nutrient pools equivalent to the N and P fertilization experiments described above and in Table 5. The simulated results were compared with observations from the field manipulations.

In the P addition experiments across the Hawaiian chronosequence, the partitioning of phosphate between microbes and mineral surfaces was well represented by the N-COM model in the intermediate (20 K yr) and old (4.1 M yr) sites (Fig. 5b and c), with no significant differences between model predictions and observations. In the youngest Hawaiian site (300 years; Fig. 5a), the relative partitioning was correctly simulated, but the predicted PO${}_{\mathrm{4}}^{\mathrm{3}-}$ magnitudes were lower than observations. Our simulations indicated that at the young soil site the added P exceeded microbial demand, resulting in lower predicted microbial P uptake than observed. This discrepancy reflected a possible deficiency of first-order SOC decay models (as we used here), which implicitly treat microbes as a part of soil organic matter. Since microbial nutrient immobilization is strictly regulated by the SOC turnover rate in this type of model, external nutrient inputs will no longer affect microbial nutrient uptake if the inputs exceed potential microbial demand. We therefore believe that explicit microbe–enzyme models might be able to better explain the strong microbe PO${}_{\mathrm{4}}^{\mathrm{3}-}$ uptake signal observed at the young Hawaii fertilization experiment site. Microbial models explicitly simulate the dynamics of microbial biomass, which might be able to capture the expected rapid growth of microbial communities under conditions of improved substrate quality (Kaspari et al., 2008; Wieder et al., 2009).

In the Puerto Rican Luquillo forest nitrogen addition experiments, partitioning of added ammonium between plants and heterotrophic bacteria was well captured by the N-COM model, with no significant differences between model predictions and observations (Fig. 5d). However, the model underestimated nitrifier NH${}_{\mathrm{4}}^{+}$ uptake. NO${}_{\mathrm{3}}^{-}$ competition in this site was also relatively accurately predicted (Fig. 5e), although the measurements did not include denitrification. Model estimates of plant NO${}_{\mathrm{3}}^{-}$ uptake and microbial NO${}_{\mathrm{3}}^{-}$ immobilization were consistent with the observed ranges, but we highlight the large observational uncertainties, particularly for microbial NO${}_{\mathrm{3}}^{-}$ uptake.

In the pseudo-first-order decomposition model we applied here to demonstrate the ECA competition methodology, the soil organic matter C : N : P ratio also limited microbial N and P uptake. For this type of decomposition model, stoichiometric differences between soil organic matter and microbes are not dynamically simulated. Such a simplification of soil and microbial stoichiometry favors large spatial scale model structures over long temporal periods, but hampers prediction of microbial short-term responses to N and P fertilization. For example, the observed difference between microbial and soil C : P ratios can be as large as 6-fold (Mooshammer et al., 2014; Xu et al., 2013). Were that the case in the observations we applied, the potential soil P demand calculated based on a fixed soil organic matter C : P ratio could be only 17 % of that based on microbial C : P ratio.

Implications of ECA competition treatment

Terrestrial ecosystem growth and function are continuously altered by climate (e.g., warming, drought; Chaves et al., 2003; Springate and Kover, 2014), external nutrient inputs (e.g., N deposition; Matson et al., 1999, 2002), and atmospheric composition (e.g., CO${}_{\mathrm{2}}$ concentration; Norby et al., 2010; Oren et al., 2001; Reich et al., 2006). Improved understanding of the underlying mechanisms regulating ecosystem responses to environmental changes has been obtained through in situ level to large-scale and long-term manipulation experiments. For example, decade-long Free-Air Carbon Dioxide Enrichment (FACE) experiments have revealed that nitrogen limitation diminished the CO${}_{\mathrm{2}}$ fertilization effect of forest (Norby et al., 2010) and grassland (Reich and Hobbie, 2013) ecosystems. However, fewer efforts have been made towards incorporating the observed process-level knowledge into ESMs. Therefore, a major uncertainty that has limited the predictability of ESMs has been the incomplete representation of soil nutrient dynamics (Zaehle et al., 2014). Even though new soil nutrient cycle paradigms were proposed during recent decades (Korsaeth et al., 2001; Schimel and Bennett, 2004), they were restricted to either conceptual models or only applied to explain laboratory experiments.

Many large-scale terrestrial biogeochemistry models (e.g., O-CN, CASA, TEM) have adopted the classical paradigm that microbes decompose soil organic matter and release NH${}_{\mathrm{4}}^{+}$ as a “waste” product (Waksman, 1931). The rate of this process is defined as “net N mineralization”, and is adopted as a “measure” of plant available inorganic N (Schimel and Bennett, 2004). This classical paradigm overlooked the fact that “net N mineralization” actually comprised two individual processes – gross N mineralization and microbial N immobilization. Implicitly, the classical paradigm assumes that the microbes have priority to assimilate as much of the available nutrient pool as possible. Soil nutrients were only available for plant uptake if there were not enough free energy materials (e.g., dissolved soil organic carbon) to support microbial metabolism. As a result, soil microbes were considered “victors” in the short-term nutrient competition. Some other large-scale terrestrial biogeochemistry models (e.g., CLM4CN), simplify the concept of nutrient competition differently. They calculate the plant N uptake and soil N immobilization separately; and then down-regulate the two fluxes according to the soil mineral N availability. As a result, plant and soil microbe competitiveness for nutrients is determined by their relative demand.

Climate-scale land models have over-simplified or ignored competition between plants, microbes, and abiotic mechanisms. In reality, under high nutrient stress conditions, plants can exude nutrient carrier enzymes or facilitate mycorrhizal fungi associations to enhance competitiveness for nutrient acquisition (Drake et al., 2011; Hobbie and Hobbie, 2006; Treseder and Vitousek, 2001). In addition, plants can adjust C allocation to construct more fine roots, which scavenge nutrients over larger soil volumes (Iversen et al., 2011; Jackson et al., 2009; Norby et al., 2004). Soil spatial heterogeneity might also contribute to the success of plant nutrient competition (Korsaeth et al., 2001). Therefore, most ecosystem biogeochemistry models with traditional treatments of nutrient competition likely underestimate plant nutrient uptake.

Nutrient competition should be treated as a complex consumer–substrate reaction network: multiple “consumers”, including plant roots, soil heterotrophic microbes, nitrifiers, denitrifiers and mineral surfaces, each competing for substrates of organic and inorganic nitrogen and phosphorus as nutrient supply. In such a model structure, the success of any consumer in substrate acquisition is affected by its consumer–substrate affinity (Nedwell, 1999). Such competitive interactions have been successfully applied to microbe–microbe and plant–microbe substrate competition modeling (Bonachela et al., 2011; Lambers et al., 2009; Maggi et al., 2008; Maggi and Riley, 2009; Moorhead and Sinsabaugh, 2006; Reynolds and Pacala, 1993) for many years.

Here, we applied the consumer–substrate network in a broader context of plant, microorganism and abiotic mineral interactions. We analyzed the consumer-substrate network using a first-order accurate equilibrium chemistry approximation (ECA; Tang and Riley, 2013; Zhu and Riley, 2015). Our sensitivity analysis confirmed that the consumer–substrate affinity and nutrient carrier enzyme abundance were the most important factors regulating relatively short-term competitive interactions. The ECA competition treatment represents ecosystem responses to environmental changes and has the potential to be linked to a microbe-explicit land biogeochemistry model. The approach allows competition between plants, microbes and mineral surfaces to be prognostically determined based on nutrient status and capabilities of each consumer.

Conclusions

In this study, we developed a soil biogeochemistry model (N-COM) that resolves the dynamics of soil nitrogen and phosphorus, plant uptake of nutrients, microbial uptake and abiotic interactions. We focused on the implementation, parameterization and testing of the nutrient competition scheme that we plan to incorporate into the ESM land models CLM and ALM. We described the multiple-consumer and multiple-nutrient competition network with the Equilibrium Chemical Approximation (ECA; Tang and Riley, 2013) considering two inhibitive effects: (1) multiple substrates (e.g., NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-})$ sharing one consumer inhibits the effective binding between any specific substrate and the consumer and (2) multiple consumers (e.g., plants, decomposing microbes, nitrifiers) sharing one substrate (e.g., NH${}_{\mathrm{4}}^{+})$ lowers the probability of effective binding between any consumer and that substrate. We calibrated the model at a tropical forest site with highly weathered soil (Tapajos National Forest, Para, Brazil), using multiple observational data sets with the MCMC approach. The calibrated model compared to multiple categories of observational data was substantially improved over the prior model (Fig. 4). The seasonal dynamics of soil carbon, nitrogen and phosphorus were moderately well captured. However, our results would likely be more robust if more temporally resolved observations of carbon, nitrogen and phosphorous were available. Although the calibrated model is the best one we can derive based on limited data, several model parameters were not well converged. We therefore conclude that more work on data collection, parameter tuning and uncertainty analysis is needed.

To test the resulting model using the calibrated parameters, we applied N-COM to two other tropical forests (Hawaii tropical forest and Luquillo tropical forest) not used in the calibration process and conducted nutrient perturbation studies consistent with fertilization experiments at these sites. The results showed that N-COM simulated the nitrogen and phosphorus competition well for the majority of the observational metrics. However, the model underestimated NH${}_{\mathrm{4}}^{+}$ uptake by nitrifiers, probably due to the loosely constrained nitrification parameters that were the result of NO${}_{\mathrm{3}}^{-}$ pool size data paucity during calibration at the Brazil site (Table 4). Data sets of soil nutrient pool sizes and CO${}_{\mathrm{2}}$ and N${}_{\mathrm{2}}$O effluxes with high-frequency sampling would significantly benefit the model uncertainty reduction.

To date, many terrestrial ecosystem biogeochemistry models assume that microbes outcompete plants and immobilize nutrients first (Wang et al., 2007; Zaehle and Friend, 2010; Zhu and Zhuang, 2013), although CLM currently assumes constant and relative demand competitiveness of plants and microbes. Few models, to our knowledge, consider the role of abiotic interactions in the competitive interactions. In the case of microbes outcompeting plants, the plant is only able to utilize the nutrients that exceed microbial demands during that time step. The leftover nutrients are defined as net mineralization, which is a widely adopted concept in soil biogeochemistry modeling (Schimel and Bennett, 2004). These models oversimplify plant–microbe interactions by imposing dubious assumptions (e.g., microbes always win against plants). We showed that (in Sect. 3.1) “competitiveness” is a dynamic rather than fixed property of the competition network, and more importantly, it should be linked to competitor functional traits (affinity and nutrient carrier enzyme abundance).

This study is an important step towards implementing more realistic nutrient competition schemes in complex climate-scale land models. Traditional ESMs generally lack realistic soil nutrient competition, which likely biases the estimates of terrestrial ecosystem carbon productivity and biosphere–climate feedbacks. This study showed the effectiveness of ECA kinetics in representing soil multiple-consumer and multiple-nutrient competition networks. Offline calibration and independent site-level testing is critically important to ensuring the newly incorporated model will perform reasonably when integrated in a complex ESM. To this end, we provide a universal calibration approach using MCMC, which could in the future be used to further constrain N-COM across plant functional types, climate and soil types.

CNP fluxes

The fluxes of carbon, nitrogen and phosphorus coming from the upstream pool ($i)$ to the downstream pool ($j)$ due to SOM decomposition are calculated as $$\begin{array}{}& & F_{\mathrm{C},ij}^{\text{move}}=f_{ij}{F}_{\mathrm{C},i}^{\mathrm{dec}}{g}_i& & F_{\mathrm{N},ij}^{\text{move}}=f_{ij}{F}_{\mathrm{C},i}^{\mathrm{dec}}min\left({\displaystyle \frac{\mathrm{1}}{{\mathrm{CN}}_i}},,,{\displaystyle \frac{g_i}{{\mathrm{CN}}_j}}\right)& & F_{\mathrm{P},ij}^{\text{move}}=f_{ij}{F}_{\mathrm{C},i}^{\mathrm{dec}}min\left({\displaystyle \frac{\mathrm{1}}{{\mathrm{CP}}_i}},,,{\displaystyle \frac{g_i}{C{P}_j}}\right),\end{array}$$ where $g_i$ is the percentage of carbon remaining in the soil after decomposition of the $i$th SOM pool (i.e., CUE, with the rest being released as CO${}_{\mathrm{2}})$; $f_{ij}$ is the fraction of SOM leaving the $i$th pool and entering the $j$th pool; and $F_{\mathrm{C},i}^{\mathrm{dec}}$ is the first-order decay of the $i$th SOM pool. CN and CP are soil C : N and C : P ratios, respectively.

If the upstream-decomposed soil organic nitrogen (phosphorus) is more than enough to sustain the downstream C : N (C : P) ratio, then the excess nitrogen (phosphorus) enters the soil NH${}_{\mathrm{4}}^{+}$ (PO${}_x)$ pool. PO${}_x$ represents the sum of PO${}_{\mathrm{4}}^{\mathrm{3}-}$, HPO${}_{\mathrm{4}}^{\mathrm{2}-}$ and H${}_{\mathrm{2}}$PO${}_{\mathrm{4}}^{-}$ that could be utilized by plants and microorganisms, and adsorbed by mineral surfaces: $$\begin{array}{}& & F_{\mathrm{N},ij}^{\text{mob}}=f_{ij}{F}_{\mathrm{C},i}^{\mathrm{dec}}max\left({\displaystyle \frac{\mathrm{1}}{{\mathrm{CN}}_i}},-,{\displaystyle \frac{g_i}{{\mathrm{CN}}_j}},,,\mathrm{0}\right)& & F_{\mathrm{P},ij}^{\text{mob}}=f_{ij}{F}_{\mathrm{C},i}^{\mathrm{dec}}max\left({\displaystyle \frac{\mathrm{1}}{{\mathrm{CP}}_i}},-,{\displaystyle \frac{g_i}{{\mathrm{CP}}_j}},,,\mathrm{0}\right),\end{array}$$ where $F_{\mathrm{N},ij}^{\text{mob}}$ and $F_{\mathrm{P},ij}^{\text{mob}}$ are the nitrogen and phosphorus gross mineralization rates. Equations (A4) and (A5) ensure that gross mineralization is not less than zero. In contrast, if nitrogen (phosphorus) is insufficient, soil microbes immobilize free NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ (PO${}_x)$: $$\begin{array}{}& & F_{\mathrm{N},ij}^{\text{immob,pot}}=f_{ij}{F}_{\mathrm{C},i}^{\mathrm{dec}}max\left({\displaystyle \frac{g_i}{{\mathrm{CN}}_j}},-,{\displaystyle \frac{\mathrm{1}}{{\mathrm{CN}}_i}},,,\mathrm{0}\right)& & F_{{\mathrm{NH}}_{\mathrm{4}},ij}^{\text{immob,pot}}=F_{\mathrm{N},ij}^{\text{immob,pot}}\cdot {\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}\right]}{\left[{\mathrm{NH}}_{\mathrm{4}}\right]+\left[{\mathrm{NO}}_{\mathrm{3}}\right]}}& & F_{{\mathrm{NO}}_{\mathrm{3}},ij}^{\text{immob,pot}}=F_{\mathrm{N},ij}^{\text{immob,pot}}\cdot {\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{\left[{\mathrm{NH}}_{\mathrm{4}}\right]+\left[{\mathrm{NO}}_{\mathrm{3}}\right]}}& & F_{\mathrm{P},ij}^{\text{immob,pot}}=f_{ij}{F}_{\mathrm{C},i}^{\mathrm{dec}}max\left({\displaystyle \frac{g_i}{{\mathrm{CP}}_j}},-,{\displaystyle \frac{\mathrm{1}}{{\mathrm{CP}}_i}},,,\mathrm{0}\right),\end{array}$$ where $F_N^{\text{immob,pot}}$, $F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{immob,pot}}$, $F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{immob,pot}}$ and $F_{\mathrm{P}}^{\text{immob,pot}}$ are microbial N, NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$ and PO${}_x$ immobilization rates. ($\mathrm{NH}$${}_{\mathrm{4}}$) and ($\mathrm{NO}$${}_{\mathrm{3}}$) are the free NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ pools, respectively. We assume that microbes have no preference for NH${}_{\mathrm{4}}^{+}$ or NO${}_{\mathrm{3}}^{-}$ (Eqs. A7, A8). If soil nutrients are limited, a limitation factor will be applied to those potential soil decomposition CNP fluxes (Eqs. A1–A9) to maintain the soil organic matter CNP stoichiometry.

MCMC chain. Blue line represents the MCMC samples that are used to infer our model posterior parameters. Two other replicated MCMC calibrations (with different random number seeds) were conducted (yellow and red lines), in order to check the convergence of MCMC calibration.

[Figure omitted. See PDF]

Gelman–Rubin convergence criterion (solid lines) calculated from three chains in Fig. A1. Baseline value is set to 1.1 (dashed lines). When the Gelman–Rubin criterion is smaller than or equal to 1.1, the chains are thought to converge.

[Figure omitted. See PDF]

Besides decomposing microbe nutrient immobilization, other potential nutrient uptakes are $$\begin{array}{}& & F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{nit,pot}}=\left(\right[{\mathrm{NH}}_{\mathrm{4}}]\cdot k_{\text{nit}}\cdot r_{\mathit{\theta }}\cdot r_T\cdot (\mathrm{1}-f^{\text{anox}})& & F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{den,pot}}=min\left(f\right(\text{decomp}),f(\left[{\mathrm{NO}}_{\mathrm{3}}\right]\left)\right)\cdot f^{\text{anox}}& & F_{\mathrm{P}}^{\text{surf,pot}}={\displaystyle \frac{{\text{VMAX}}_{\mathrm{P}}^{\text{surf}}\cdot K_M^{\text{surf,P}}}{(K_M^{\text{surf,P}}+[{\mathrm{PO}}_x]{)}^{\mathrm{2}}}}\cdot {\displaystyle \frac{\mathrm{d}\left[{\mathrm{PO}}_x\right]}{\mathrm{d}t}},\end{array}$$ where $F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{nit,pot}}$, $F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{den,pot}}$ and $F_{\mathrm{P}}^{\text{surf,pot}}$ are potential rates for NH${}_{\mathrm{4}}^{+}$ nitrification, NO${}_{\mathrm{3}}^{-}$ denitrification and mineral surface PO${}_x$ adsorption. $k_{\text{nit}}$ is the maximum fraction of free NH${}_{\mathrm{4}}^{+}$ pool that could be utilized by nitrifiers. The potential nitrification rate is controlled by soil temperature ($r_T)$, soil moisture ($r_{\mathit{\theta }})$ and soil oxygen status ($\mathrm{1}-f^{\text{anox}})$. The potential denitrification rate ($F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{den,pot}})$ is either constrained by substrate availability ($f$(decomp)) or NO${}_{\mathrm{3}}^{-}$ availability ($f\left(\right[{\mathrm{NO}}_{\mathrm{3}}\left]\right)$; Del Grosso et al., 2000), taking into account the soil anaerobic condition ($f^{\text{anox}})$. $F_{\mathrm{P}}^{\text{surf,pot}}$ is derived from the Langmuir adsorption model (Barrow, 1978), where adsorbed P is equal to VMAX${}_{\text{P}}^{\text{surf}}\cdot \frac{\left[{\mathrm{PO}}_x\right]}{K_M^{\text{surf,P}}+\left[{\mathrm{PO}}_x\right]}$. Taking the time derivative leads to the adsorption rate (Wang et al., 2010).

Posterior model parameters (blue bars) fitted to Gaussian distribution (red line).

[Figure omitted. See PDF]

Soil NH${}_{\mathrm{4}}^{+}$ content is altered by inputs from deposition ($F_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{dep}})$ and biological N${}_{\mathrm{2}}$ fixation ($F^{\text{BNF}}$; Cleveland et al., 1999): $$F^{\text{BNF}}=1.8\cdot {\displaystyle \frac{\mathrm{1}-e^{-0.003\cdot {\text{NPP}}_{\text{annual}}}}{365\cdot 86400}},$$ where NPP${}_{\text{annual}}$ is annual net primary production. Controls on biological N${}_{\mathrm{2}}$ fixation are complex and several models have been developed for large-scale land BGC models (Cleveland et al., 1999; Fisher et al., 2010; Hartwig, 1998; Parton et al., 1993; Running et al., 1989; Vitousek and Field, 1999). However, the emergent responses predicted across these model structures are inconsistent (Galloway et al., 2004). Recognizing this important structural uncertainty, we used a simple model where biological N${}_{\mathrm{2}}$ fixation ($F^{\text{BNF}})$ is modeled as a function of annual NPP (Cleveland et al., 1999).

Soil NO${}_{\mathrm{3}}^{-}$ content is modified by external deposition inputs ($F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{dep}})$ and leaching losses ($F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{leach}})$: $$F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{leach}}={\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}\right]}{W}}\cdot Q_{\text{dis}},$$ where soil nitrate concentration ([$\mathrm{NO}$${}_{\mathrm{3}}$]: gN m${}^{-\mathrm{2}})$ divided by soil water content ($W$: gH${}_{\mathrm{2}}$O m${}^{-\mathrm{2}})$ results in the concentration of dissolved nitrate (DIN). The hydrologic discharge ($Q_{\text{dis}}$: gH${}_{\mathrm{2}}$O m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}})$ applied to DIN (gN gH${}_{\mathrm{2}}$O${}^{-\mathrm{1}})$ leads to the leaching loss (gN m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}})$.

Soil PO${}_x$ content is affected by external inputs from parent material weathering ($F^{\text{weather}})$ and leaching losses ($F_{\mathrm{P}}^{\text{leach}})$. Sorbed P (P${}_{\text{S}})$ could be further strongly occluded and become unavailable for plant and microbial uptake. Parent material stock can be increased by atmospheric dust deposition ($F_{\mathrm{P}}^{\text{dep}})$ (Mahowald et al., 2008): $$\begin{array}{}& & F^{\text{weather}}=\left[P_P\right]\cdot k_{\text{weather}}& & F_{\mathrm{P}}^{\text{leach}}={\displaystyle \frac{\left[{\mathrm{PO}}_x\right]}{W}}\cdot Q_{\text{dis}}& & F_{\mathrm{P}}^{\text{occl}}=\left[P_{\text{S}}\right]\cdot k_{\text{occl}},\end{array}$$ where parent material weathering ($F^{\text{weather}})$ is calculated using a weather rate ($k_{\text{weather}})$ and parent material P content ([P${}_{\text{P}}$]). PO${}_x$ leaching loss is modeled with a similar approach to nitrate leaching (Eq. A16). The phosphorus occlusion rate is modeled as the product of a constant rate ($k_{\text{occl}})$ and the sorbed P content ([P${}_{\text{S}}$]).

Derivation of VMAX

The enzyme substrate reaction is $\mathrm{S}+E\underset{k_{\mathrm{1}}^{-}}{\overset{k_{\mathrm{1}}^{+}}{⟷}}\mathrm{C}\stackrel{k}{⟶}\mathrm{P}+E$, where the enzyme ($E$) and substrate (S) reaction is reversible and forms complex (C). The irreversible reaction releases product (P) and liberates enzyme ($E$). At steady state, the formation rate of the enzyme substrate complex is equal to the consumption rate: $$k_{\mathrm{1}}^{+}\left[\mathrm{S}\right]\left[\mathrm{E}\right]=k_{\mathrm{1}}^{-}\left[\mathrm{C}\right]+k\left[\mathrm{C}\right].$$ To simply the equation, we define an affinity parameter: $$K_M={\displaystyle \frac{k_{\mathrm{1}}^{-}+k}{k_{\mathrm{1}}^{+}}}={\displaystyle \frac{\left[\mathrm{S}\right]\cdot \left[E\right]}{\left[\mathrm{C}\right]}}.$$ By definition, the total enzymes [E${}_{\text{tot}}]$ in the system is the sum of free enzymes [E] and enzymes that are bound with the substrate [C]: $$\left[E_{\text{tot}}\right]=\left[E\right]+\left[\mathrm{C}\right].$$ Substituting Eq. (B3) into (B2), we have $$K_M={\displaystyle \frac{\left[\mathrm{S}\right]\cdot \left(\right[E_{\text{tot}}]-[\mathrm{C}\left]\right)}{\left[\mathrm{C}\right]}}.$$ Collecting terms containing [C], we have $$\left[\mathrm{C}\right]\cdot (K_M+[\mathrm{S}\left]\right)=\left[E_{\text{tot}}\right]\cdot \left[\mathrm{S}\right].$$ The production rate is $${\displaystyle \frac{d\left[\mathrm{P}\right]}{\mathrm{d}t}}=k\cdot \left[\mathrm{C}\right].$$ Substituting Eq. (B5) into (B6), we have $${\displaystyle \frac{d\left[\mathrm{P}\right]}{\mathrm{d}t}}=k\cdot \left[E_{\text{tot}}\right]\cdot {\displaystyle \frac{\left[\mathrm{S}\right]}{K_M+\left[\mathrm{S}\right]}}.$$ Comparing Eq. (B7) with the classic Michaelis–Menten equation, it is clear that the definition of maximum production rate is the product of the reaction rate and enzyme abundance in the system: $$\text{VMAX}=k\cdot \left[E_{\text{tot}}\right].$$

Acknowledgements

This research was supported by the Director, Office of Science, Office of Biological and Environmental Research of the US Department of Energy under Contract no. DE-AC02-05CH11231 as part of the Regional and Global Climate Modeling (RGCM) and ACME programs. Edited by: A. Rammig