A “toy” model of biogeochemical dynamics on climate gradients

Abstract

We used a simple “toy” model to aid in the evaluation of the controls of biogeochemical patterns along a climate gradient. The model includes simplified treatments of water balance (precipitation minus Potential Evapotranspiration), leaching, weathering of cation- and P-bearing minerals, N cycling and loss, biomass production, and biological N fixation. We use δ¹⁵N as a central integrator of biogeochemical processes, because δ¹⁵N integrates multiple pathways of N input, output, and transformation in ecosystems. The model simulated the location and magnitude of a peak in δ¹⁵N on a gradient on Kohala Volcano, Hawai‘i which peaked ~  + 14 ‰ in sites receiving ~ 3.5 cm/month average precipitation (− 1300 mm/year water balance); the model also captured a peak in total P in surface soil at intermediate levels of precipitation and water balance, and other biogeochemical features on the gradient. We then applied the model to understanding the patterns of and mechanisms underlying nutrient limitation to net primary production (NPP) and plant biomass on the gradient, testing for the existence and extent of N and P limitation by simulated additions of N and/or P in the model. Both a simulated symbiotic biological N fixer and a simulated non-fixer were limited by P supply across the gradient; the non-fixer was independently limited by N supply in wetter sites. By running the toy model with and without the influence of temperature, we demonstrated that water is the most important factor shaping biogeochemical patterns on this gradient.

Introduction

A low supply of biologically available nitrogen (N) often constrains primary productivity and other ecosystem processes across a wide range of ecosystems. This constraint—termed N limitation—is demonstrated by the fact that additions of N often increase rates of productivity and other processes, both alone and interactively with other elements. Moreover, N additions enhance productivity in more terrestrial ecosystems and do so to a greater extent than do additions of any other nutrient (Elser et al. 2007; LeBauer and Treseder 2008; Fay et al. 2015; Hogberg et al. 2017). The widespread importance of N limitation begs the question “What controls the availability and cycling of N?”.

Absent supplementation by fertilization, anthropogenic atmospheric deposition, or similar processes, N supply does not exert independent control over terrestrial ecosystems. Rather, the pools and fluxes of N are consequences of multiple other processes operating on multiple time scales and hierarchical levels of control. These controls range from proximate ones (such as biological N fixation and its controls, including the supply of weathering-derived elements like P) that may themselves be influenced by other controlling processes, to ultimate, independent controls that may be remote from immediate influence on N cycling. Following the framework of Jenny (1980), terrestrial ecosystems can be conceptualized as being controlled (ultimately) by a few independent “state factors”, the most important of which Jenny identified as climate, relief or topography, parent material, organisms (defined as the species pool that could be present in a site), time, and human activity. In this framework, state factors ultimately control the development and properties of soils and the composition and functioning of biological communities (the species actually present on a site, and their dynamics). These soil and biotic processes interact and (at times together with the state factors themselves) control the inputs, losses, pools, and dynamics of N within terrestrial ecosystems.

In this paper, we evaluate N dynamics and other soil properties along a well-defined and well characterized climate gradient (climosequence, in Jenny’s terminology) in the Hawaiian Islands, using natural abundance N (δ¹⁵N) as an integrated measure of N dynamics and a “toy” model that simulates N dynamics and other ecosystem properties. A number of earlier studies evaluated N dynamics on climate gradients—often using natural abundance N isotope measurements as an integrated measure of cycling processes (Austin and Sala 2002; Houlton et al. 2006; Wang et al. 2014), and global syntheses have taken a similar approach (Amundson et al. 2003; Craine et al. 2009). These studies generally find an increase in δ¹⁵N (to increasingly enriched values) in progressively lower-rainfall sites. We complement that approach with a “toy” model designed to provide explanations for that increase here. A “toy” model is defined as a simplified set of objects and equations that can be used to understand a mechanism; we find developing such a model particularly useful for testing our understanding of the implications of biogeochemical processes. Developing a simple model forces us to write equations that express our understanding, running such a model shows whether our understanding of mechanisms (when pulled together in the toy model) actually yields the patterns we observe in the field. The number of variables and constants that we simulate in the model means that changes in reasonable parameter values can influence the position and magnitude of the less-central results of the model (which include essentially everything but δ¹⁵N); accordingly, we are more interested in evaluating patterns of variation than in absolute values. Here, we evaluate biogeochemical properties and processes along a well-defined, well-studied climate gradient on Kohala Volcano, Hawai‘i, along which rainfall varies from 285 to 3238 mm/year, precipitation minus potential evapotranspiration (water balance) from − 1581 to + 1839 mm/year, and temperature from 16.1 to 23.3 °C (Vitousek and Chadwick 2013; von Sperber et al. 2017). The main vegetation/land use on this gradient is now pasture, dominated by the closely related C₄ pasture grasses Pennisetum cladestinum and Cenchrus ciliaris; for most of its history, the gradient has been dominated by native forest that varied from wet rainforest through mesic forest to dry forest and shrubland as a function of climate (Chadwick et al. 2007). One result of empirical studies on this gradient is that there are striking non-linearities and discontinuities in biogeochemical properties that emerge along this continuous gradient, where biogeochemical properties change substantially with a small increment in forcing at a particular place on a continuous gradient. We term these discontinuities “pedogenic thresholds” (Chadwick and Chorover 2001). Complementarily, there are substantial portions of the gradient where biogeochemical properties and processes change little despite substantial variation in forcing. We term these regions “soil process domains” (Vitousek and Chadwick 2013).

Initially, we compared model outputs with soil and ecosystem properties that we had measured along the gradient; we then used the model to evaluate the mechanisms underlying constraints to biological N fixation and nitrogen (N) and phosphorus (P) limitation to primary production on the gradient. Finally, we applied this model to another climate gradient on Mauna Kea Volcano, Hawai‘i; there, rainfall varies from 439 to 5202 mm/year, water balance from − 1267 to + 3810 mm/year, and mean annual temperature from 5.7 to 14.2 °C (Bateman et al. 2019) (all climate information from Giambelluca et al. 2013). The gradients differed in that the Mauna Kea (MK) gradient developed in ~ 20,000 year old substrate (Bateman et al. 2019), while the Kohala (KH) gradient was developed in ~ 150,000 year old substrate (Chadwick et al. 2003), and in that the driest sites were warm leeward sites on KH while the driest sites were cold high elevation sites on MK. Temperature, precipitation, and water balance were very highly correlated within each gradient, but we could use the model to separate the effects of temperature vs precipitation (and water balance) on soil properties and N dynamics (especially δ¹⁵N natural abundance) on both gradients.

Methods

Toy model

We wrote the program that runs the toy model in Matlab (version 2018a, Mathworks Inc.); the code is available as on-line supplementary information associated with this article, together with a table in which variables and constants are listed and defined. The basic structure of the model is outlined in Fig. 1; in practice, the program that runs the model consists of an outer loop in average rainfall, which (for KH) varies from 2 cm/month of rain to 27 cm/month in 1 cm/month increments, and an inner loop in time (with time steps assumed to be months). The time loop at each average rainfall starts with all mineral pools fully unweathered and no plant biomass or soil organic matter; in the most widely used implementation, it runs for 60,000 time steps. For the first 12,000 time steps of this implementation, it runs with the same rainfall each month; after 12,000 time steps, we introduce random month-to-month variation (with the appropriate mean rainfall for that position on the gradient) for the remaining 48,000 time steps. Outputs of interest are averaged for the last 10,000 time steps at each point in the rainfall loop. These particular time steps were selected as a compromise between the time it takes the model to equilibrate and the need to integrate the results of a stochastic simulation (see below) across long time intervals, and the time required to run the model. In the program, we described and implemented simplified versions of several aspects of ecosystem dynamics, including:

Fig. 1

The overall structure of the toy model. More detailed summaries for the submodels focused upon water, P, N, and productivity/biomass are displayed in subsequent figures. DOC, DON, and DOP are dissolved organic C, N, and P

Water

The structure of the portion of the toy model concerned with water is summarized in Fig. 2. As described above, the program (for KH) includes a loop from ~ 2 cm/month average precipitation to 27 cm/month. Month-to-month variation in precipitation is calculated as:

Fig. 2

The structure of the toy model with respect to water pools and fluxes

For: v ≤ 20

z = v − 2

$$PPT\left( {v,y} \right) = meanPPT\left( {\text{v}} \right) + \left[ {meanPPT\left( {\text{v}} \right)* \left( {11 - \frac{{\text{z}}}{1.8}} \right){ *}\left( {{\text{rand}} - \left( {0.7 - \frac{{{\text{z}} - 2}}{90}} \right)} \right)} \right]$$

For v > 20 (up to the maximum rainfall of 27 cm/month)

$$PPT\left( {v,y} \right) = meanPPT({\text{v)}} + \left[ {meanPPT\left( {\text{v}} \right)*\left( {rand - 0.5} \right)} \right]$$

where PPT(v, y) is precipitation for that month at that location on the gradient in cm, meanPPT(v) is mean monthly precipitation (which loops from 2 to 27 cm/month), v is the variable that loops (mean monthly precipitation), z is a variable indexed to v, y is the variable representing time steps, and rand is a pseudorandom number evenly distributed in the interval 0 to 1. This formulation includes greater month-to-month variability in precipitation in drier sites, which typically is observed in practice and in this case is designed to match the greater month-to-month variation in precipitation in low-rainfall sites observed across the same climate gradient over a 25 year period on a neighboring ranch (Pono von Holt, personal communication; Marshall et al. 2017). The coefficients in the equation were selected empirically to match the pattern in precipitation variability that was observed. The full program can be run without month-to-month variation in precipitation, so the importance of that variation for simulated ecosystem dynamics can be assessed.

Soil water-holding capacity (WHC) to an implicit depth of 50 cm is assumed to be 30 cm, reflecting the high WHC of volcanic soils. If soil water (stored water plus that month’s rainfall) exceeds 30 cm, the excess leaches from the surface soil. Production, decomposition, and mineral weathering are influenced by soil water content; they are not constrained when soil moisture is at WHC, and decrease linearly with decreases in soil moisture below WHC, with plant production and transpiration stopping at a water content ~ 17% of WHC, as an expression of a permanent wilting point.

Potential evapotranspiration (PET) is estimated by calculating soil surface evaporation (maximum 2 cm/month) and transpiration (maximum 13 cm/month). The total (15 cm/month, or 1800 mm/year) is similar to the relatively constant Priestley-Taylor PET reported across the KH and MK gradients in Giambelluca et al. (2013). Both surface evaporation and transpiration are constrained in the model when there is insufficient soil moisture to support potential evaporation or transpiration; transpiration ceases when the soil water content is < 17% of WHC, while evaporation ceases when at 0 soil water content.

We plot the results of both observations and simulations against both rainfall and an estimate of water balance (ppt-PET). We use water balance because it represents a calculation of the primary driver of nutrient depletion in terrestrial ecosystems (the amount of water that can leach through soil); we also use rainfall because it is measured directly and because PET varies relatively little across these climate gradients while rainfall varies substantially, so that water balance largely reflects variation in ppt.

We also simulate a deeper layer of soil, implicitly from 0.5 to 2 m depth. Water-holding capacity in this layer is assumed to be 50 cm; leaching of water from the deeper layer occurs when the sum of stored water plus water inputs from the surface layer exceeds 50 cm. We assume that plants take up water from this deep soil layer when there is insufficient water in the surface soil to support potential transpiration, and when there is sufficient water in the deep soil layer to do so. For this deep layer, we assume that the permanent wilting point occurs at 20% of water holding capacity.

Temperature

Temperature is highly correlated with precipitation on both KH and MK, so we calculated mean annual temperature from monthly precipitation for both gradients.

$$MAT_{KH} \left( v \right) = 25 - \left[ {\frac{{meanPPT\left( {\text{v}} \right)}}{2.7}} \right]$$

$$MAT_{MK} \left( v \right) = 5 + \left[ {\frac{{meanPPT\left( {\text{v}} \right)}}{4.5}} \right]$$

where MAT_KH(v) and MAT_MK(v) are mean annual temperature for each average rainfall (each site) along the KH and MK gradients, respectively. For KH, this formulation has the effect of making MAT ~ 24 °C at the warm, dry, low elevation end of the gradient and ~ 15 °C that at the cooler, wet, higher elevation end, slightly warmer (< 1 °C) at the dry end of the gradient and slightly cooler (~ 1 °C) at the wet end than the values in Giambelluca et al. (2013). For MK, this has the effect of making MAT 5.9 °C at the cold, dry, high elevation end of the gradient and 15 °C at the cool, wet, lower elevation end, again slightly more extreme than the values reported in Giambelluca et al. (2013). We did not simulate month to month variation in temperature; annual variation in mean monthly temperature typically is < 4 °C in Hawai‘i (though even that small variation could be important, especially at the low-temperature end of the MK gradient).

We calculate the effect of temperature on productivity and on mineral weathering as a linear function of temperature, with a base value (at which temperature does not constrain productivity or mineral weathering) of 25 °C. The equation for this scalar is presented below in the discussion of plant production. The effect of temperature on decomposition is calculated as an exponential (Q₁₀) function, with a Q₁₀ value of 2.5 (so that for every 10 °C decrease in temperature below 25 °C (below the highest temperature on either gradient), there is a 2.5-fold decrease in decomposition. Linear variation in productivity and exponential variation in decomposition with temperature, and these parameter values, are consistent with observations in Hawaiian ecosystems (Raich et al. 1997; Vitousek 2004).

Cations

We calculate weathering of cations (abbreviated Ca, and treated like calcium) from primary minerals and leaching from soil in a manner similar to that used in the toy model in Vitousek et al. (2016). Briefly, we assume that there are two generic classes of primary minerals (rather than the three classes in Vitousek et al. 2016) that differ in their lability (susceptibility to weathering).

Weathering of cation-containing minerals belonging to class x (of the two classes of minerals) is calculated as:

$$Weather\,Ca_{x} = {\text{weathercoeff}}_{x} *climate scalar*\left( {MaxCaconc_{x} - Caconc} \right)$$

where WeatherCa_x is the quantity of Ca weathered from minerals of class x (units of increase in concentration), weathercoeff_x is a kinetic coefficient expressing the rate at which Ca can be weathered from that class of minerals (units of quantity per time), the climate scalar expresses the combined (multiplicative) influence of temperature and soil moisture content on weathering, $MaxCa{conc}_{x}$ is the Ca concentration in solution above which no further weathering of that class of mineral takes place (units of concentration), and Caconc is the concentration of Ca in the soil at the beginning of the weathering process (for that class of minerals). We assume that weathercoeff and MaxCaconc for a given class of minerals are fully correlated; a class of minerals that is recalcitrant to weathering has both a slower kinetic coefficient and a lower equilibrium concentration of Ca than a more labile mineral. We further assume that the class of minerals that is more labile to weathering is depleted first, and that as long as it remains in the soil it suppresses the weathering of the more recalcitrant class of minerals by raising Caconc above $MaxCa{conc}_{x}$ for the more recalcitrant class of minerals. Ca is lost from the soil via leaching at a rate that depends on the Ca concentration in the soil and on the amount of water leached from the soil. We assume that an amount of leaching equal to the soil water stored at WHC removes all of the labile Ca in soil, with a proportional reduction in removal for lower quantities of water leaching. Weathering is simulated in the deep soil layer as well, but we do not simulate cation uptake by plants.

Phosphorus

The structure of the portion of the toy model that is focused on P is summarized in Fig. 3. Helfenstein et al. (2018) recently analyzed and summarized variation in P dynamics across the KH gradient. P in the toy model weathers in a fashion similar to Ca, again with two classes of P-bearing minerals that differ in lability. The P minerals are assumed to have lower kinetic coefficients for weathering and lower equilibrium concentrations than do the Ca minerals. P differs from Ca in multiple additional ways. We resolve an organic pool of soil P, and calculate mineralization of P from that pool as:

$${\text{P}}_{min} = Soilorg{\text{P}} - \frac{{Soilorg{\text{C}}}}{{{\text{CP}}_{crit} }}$$

Fig. 3

The structure of the toy model with respect to P pools and fluxes. Class 1 primary minerals are more labile to weathering than are class 2 minerals, and class 1 consists of both soil minerals and minerals derived from atmospheric deposition. Classes of P containing primary minerals and their weathering in the deep soil are not shown, although they were simulated. PO₄ can reach the deep soil via leaching from the surface soil, or be generated by weathering in the deep soil. DOP is the leaching loss of dissolved organic P

where P_min is P mineralization at that time step, SoilorgP is the pool of organic P in the soil, SoilorgC is the pool of organic C in the soil, and CP_crit is the ratio of C/P below which P is mineralized and above which it remains in SoilorgP (here CP_crit is 200). If this calculation of P_min yields a negative value, we set P_min to 0 for that time step. Also, we calculate losses of dissolved organic P (DOP) from the system at a rate dependent on a loss coefficient, on the size of the dissolved organic P pool, and on the amount of water leaching (with leaching of all soluble DOP (as defined by the loss coefficient) at or above an amount of leaching equal to WHC (so that the entire soil solution is replaced in that time step), and a linear decrease in DOP loss with leaching < WHC.

Additionally, we calculate plant uptake and cycling of P. We assume the loss coefficient for DOP is lower than the loss coefficient for C or N, as observed by Neff et al. (2000) elsewhere in Hawai‘i. We calculate atmospheric deposition of a P-containing mineral, which we assume to be derived from long-distance transport of continental dust from Asia (Chadwick et al. 1999; Vitousek 2004), and (because of the relative immobility of P in soil) we assume that a maximum of 10% of the P pool that has been weathered and is present in biologically available form as phosphate can be leached from the soil, and that 50% of this available P pool is transferred into a secondary mineral P pool. Finally, we add the atmospherically-deposited P to the more labile class of P-containing primary minerals, while secondary mineral P is relatively recalcitrant to weathering and breaks down (releasing P) at the same rate as does the more recalcitrant class of P minerals. We use the model to estimate P fractions of unweathered primary mineral P, labile P (the inorganic P in biologically-available form after weathering, plant uptake, leaching, and secondary mineral formation), organic P, and secondary mineral P in the surface soil.

We calculate weathering and leaching of P in the deep soil layer similarly to the surface layer, though we do not simulate an organic or a secondary mineral P pool in the deep layer. Uptake of P by plants in the deep soil layer is calculated from the amount of deep soil water taken up, the proportion of total deep soil water that this uptake represents, and the quantity of dissolved inorganic P in the deep soil layer.

Nitrogen

The structure of the portion of the toy model that is focused on N is summarized in Fig. 4. N enters ecosystems through atmospheric deposition (which we treat as a linear function of precipitation, as for Ca and P) and through biological N fixation (discussed below). N also can be added as fertilizer, to test for N limitation to net primary production (NPP) (or plant biomass). N is lost through ammonia volatilization, nitrate loss (which could be due to nitrate leaching or denitrification), transformation-dependent trace gas losses (defined below), and dissolved organic N (DON) leaching. For each input or loss pathway, we calculate inputs/outputs of both all N and of ¹⁵N separately, and we assume a level of isotopic abundance (for losses of N, isotopic fractionation) associated with each pathway of input/loss. For inputs, we assume no isotopic discrimination associated with fertilizer or biological N fixation (so ¹⁵N enters at the average value of 0.003663 of total N, for a δ¹⁵N of 0‰). For atmospheric deposition, we assume deposited N is 5 per mil depleted (− 5‰), based on observations of N isotopes in non-N-fixing lichens without access to soil N in Hawai‘i (Vitousek et al. 1989). For the within-ecosystem cycling of N, we calculate N mineralization as described above for P (with a critical C:N ratio (CNcrit) of 20), also setting N mineralization for that time step to 0 when the calculation yields negative N mineralization, and we calculate nitrification as a variable fraction (varying with climate, as described below) of the ammonium-N pool remaining in the soil after N mineralization and after plant N uptake by the non-N-fixer.

Fig. 4

The structure of the toy model with respect to N pools and fluxes. Nitrate (NO₃⁻) is leached from the surface layer of soil to a deeper layer in the soil; that nitrate can then be lost from the system entirely, or taken up by plants. N-fixing plants are assumed not to take up soil N, as discussed in the text. Denitrification is not included here, as discussed in the text; transformation-dependent trace gas fluxes arise from nitrification. DON is leaching loss of dissolved organic N

Pathways of N loss will be discussed individually. Ammonia volatilization is pH-dependent; pH is not simulated directly in the model, so we use the mineral Ca pool in the soil as an index of soil pH. Our reasoning is that weathering is the main buffer of acidity in the soil here (Chadwick and Chorover 2001), and high concentrations of mineral Ca in the soil reflect ongoing weathering and so lower acidity than in older or more fully leached areas. Ammonia volatilization is calculated as 10% of the ammonium pool in the surface soil layer (after plant uptake and nitrification) times the index of soil pH (so that when Ca and so pH is low, ammonia volatilization approaches 0); we assume that the heavy isotope of N (¹⁵N) is discriminated against by 1.5% (so a fractionation of -15 ‰ for ammonia volatilization), making the remaining N in the system ¹⁵N-enriched.

Nitrate loss could occur by either leaching of nitrate through the soil or through denitrification. We believe that on the KH gradient, most of the nitrate loss is via leaching. Our reasoning is that (1) nitrate pools in the soil are very low in the high-rainfall sites on the climate gradient (von Sperber et al. 2017), an observation consistent with either leaching or denitrification; and (2) δ¹⁵N of the N in the small amount of nitrate in high-rainfall sites is strongly ¹⁵N-depleted (von Sperber et al. 2017) a result consistent with nitrification of only a fraction of the soil ammonium pool and subsequent leaching losses, but inconsistent with denitrification of a substantial fraction of the nitrate, because denitrification would leave the residual nitrate–N ¹⁵N-enriched (Houlton et al. 2006). Accordingly, we calculate rates of nitrification that decline from low-rainfall into high-rainfall sites, and we estimate nitrate losses as the fraction of soil water at WHC that leaches each month (up to leaching = WHC, above that point we assume no further increase in leaching) times the nitrate pool in the soil. We assume that the heavy isotope of N (¹⁵N) is discriminated against by 1.0% (so a fractionation of − 10 ‰ for nitrification and subsequent nitrate losses), thereby making the remaining N in the system ¹⁵N-enriched. This assumption is consistent with von Sperber et al.’s (2017) measurements of ¹⁵N in nitrate; this assumed fractionation is less than reports of kinetic isotope fractionation during nitrification (Mariotti et al. 1981), but the co-occurrence of nitrification (by multiple pathways) and denitrification in real ecosystems (Granger and Wankel 2016), together with the extent to which reactions go to completion, makes applying kinetic isotope fractionations to real ecosystems problematic. Episodic denitrification could be more important than it was during von Sperber et al.’s sampling, and so greater than we assume here; if so, it is possible that the fractionation associated with nitrate loss could be greater. While Incomplete nitrification leads to an ¹⁵N depleted nitrate pool (as von Sperber et al. (2017) observed), and denitrification would then fractionate that pool a little less than did nitrification (Mariotti et al. 1981), with a lower overall fractionation for the whole process from ammonium to N₂ than for nitrate leaching. However, more rapid and complete nitrification (at the limit, of the whole ammonium pool) would not lead to fractionation of nitrate relative to ammonium, and denitrification could then cause a greater overall isotopic fractionation associated with nitrate loss.

Transformation-dependent losses of N-containing trace gases represent our effort to represent losses of N that occur during nitrification, we estimate this flux at 1% of nitrification (Firestone and Davidson 1989). We consider those losses to differ in kind from ammonia volatilization and nitrate leaching, both of which occur from accumulated pools of biologically-available N that remain in the soil after organisms have acquired the N that they can. In this sense, we regard both ammonia volatilization and nitrate leaching to be losses that are controllable by organisms, in that those organisms (if their activity would otherwise be constrained by a shortage of N) could prevent those losses; we regard transformation-dependent losses to be uncontrollable, in that they will occur if there is nitrification even if N is in short supply in the ecosystem. In practice, simulated transformation-dependent N trace gas losses represent a small fraction of N losses from unfertilized sites across the KH gradient.

Losses of DON also are assumed to represent uncontrollable losses of N (Hedin et al. 1995), in that the organic N compounds that leach through soils are large-molecular-weight recalcitrant compounds that are not readily taken up by organisms (Perakis and Hedin 2001, 2002). We calculate losses of DON from all the sites on the KH gradient as representing a maximum of 0.1% of soil organic N per month times the fraction of soil water at WHC that is leached each month (again with no further increase in DON loss if leaching is greater than WHC during that month).

Uptake and loss of N from a deep soil layer is calculated similarly to P uptake and loss, as described above. We calculate an ¹⁵N balance only for the surface layer; consequently, we treat the uptake of N from the deep soil layer as an input of ¹⁵N to the surface layer (having calculated the leaching of ¹⁵N from surface to deep as a fractionating loss of ¹⁵N from that layer).

Plant biomass, productivity, and biological N fixation

The structure of the portion of the toy model that is concerned with plant growth and C cycling is summarized in Fig. 5. We simulate plant productivity following the procedure in Vitousek and Field (1999), adjusted to a monthly rather than an annual time step by dividing their maximum potential productivity and the N and P supply necessary to support that productivity by 12. We simulate the growth both of a non-N-fixing plant and of a symbiotic N-fixing plant; the two differ in that the N-fixer has an N- and P-rich stoichiometry relative to the non-fixer (McKey 1994; Menge 2011) (C:N:P for the non-fixer is 500:10:1, while that for the N-fixer is 333.3:10:1), and in that the N-fixer is an obligate fixer (it derives all of its N from fixation) (Menge et al. 2009). Both the fixer and non-fixer are assumed to be perennial plants whose tissues (and N and P contents) turn over at a rate of 1%/month; accordingly, plant biomass is a multiple-month integration of net primary production (NPP). For the non-fixer, potential productivity per month is assumed to represent an (implicitly light-limited) maximum, reduced in practice by the supply of N or P (whichever is in shorter supply relative to plant demand) and by a climate scalar calculated based on multiplicative constraint by soil water content and by temperature. That constraint is calculated as:

$$prodscalar1\left( y \right) = \left( {\frac{Soilwater1\left( y \right)}{3}} \right)*\left( {\frac{Temp\left( v \right)}{{25}}} \right)$$

Fig. 5

The structure of the toy model with respect to plant productivity, plant biomass, decomposition, and N fixation. The plant pools and fluxes are dependent on the water, P, and N submodels summarized in Figs. 2, 3 and 4. Available pools of these resources in soil are represented by boxes bounded by dashed lines, and the uptake of these resources by plants is indicated by dashed arrows, as is biological N fixation. DOC is the leaching loss of dissolved organic C

where prodscalar1(y) is the climatic constraint to potential productivity, y is the month of simulation, soilwater1(y) is the water content of the soil (in cm of water/10—additionally, prodscalar1(y) is set to 0 when soilwater1(y) < 0.5, representing a permanent wilting point), 3 represents the WHC of the soil/10 (30 cm of water), and Temp(v) is the mean temperature (in °C) at that point on the climate gradient.

We calculate a cost of N acquisition for the non-fixer as in Vitousek and Field (1999); that cost is assumed to be low when N is abundant (in acquiring the first half of N available in the soil), but the cost increases as the N pool is depleted.

For the N-fixer, we assume a fixed cost of N fixation (30% of plant NPP goes to support N fixation). The non-fixer makes use of all of the potential productivity up to the point where its cost for acquiring N exceeds that of the N-fixer. Beyond that point, the N fixer (alone) grows and derives all of its N from fixation. The growth of the N fixer is independent of soil N supply; however, the N-fixer can be constrained by the remaining productive capacity of the site (implicit light limitation) after use of a portion of that capacity by the non-fixer, by climate, and/or by P supply (Vitousek and Field 1999). In addition, the growth of both plants is constrained when their biomass is low (they cannot jump all the way to their potential productivity in one time-step), and the growth of the N-fixer is further constrained by a function representing shading by the non-fixer (which reduces the energy available to support N fixation). Both of these constraints are calculated as in Vitousek and Field (1999), adjusted to the monthly time step of the present toy model.

For both plants, we calculate several possible rates of production based on potential productivity (implicitly light-limited) and the productivity that can be supported by biologically-available P and N—for non-fixers (the productivity of N-fixer can be constrained by P but not N). All of these estimates of productivity can be further constrained by soil water, temperature, a low initial biomass, and for the N fixer by shading by the non-fixer. The minimum of these estimates of productivity is then selected and applied to calculate uptake of N and P based on the fixed C:N:P stoichiometry of each plant. Inputs of N through biological N fixation are calculated as the NPP of the N-fixer times its (fixed) N:C ratio; biologically fixed N becomes available in the system as a whole when litter from the N-fixer enters the soil organic C and N pools, and is mineralized as described above. Also, atmospheric inputs of N are assumed to be biologically available.

The uptake of water from the deep soil layer is determined by calculating the deficit (if any) in actual transpiration relative to potential transpiration. We assume water equivalent to that deficit is taken up from the deep soil layer—if it is available there. If the water required to support potential transpiration can be taken up from the surface layer, then there is no uptake of water from the deep layer. We also calculate the amount of N and P taken up from the deep layer in proportion to their concentration (in available forms) there, and the fraction of the pool of deep soil water that is taken up. In effect, we are assuming that deep roots in the soil have the primary function of obtaining water, when it is in short supply in the surface soil, and that their acquisition of nutrients is incidental to that primary function. We then take the minimum of the potential increment in productivity supported by the uptake of water from deep soil, and of N or P from deep soil. If there is additional P taken up from deep soil in excess of the non-fixers requirement (its ratio of C:N:P), then that P is acquired by the N fixer, and its increment in NPP and in N fixation is calculated from the stoichiometry (C:N:P ratio) of the N fixer.

Results and discussion

We simulated ecosystem dynamics along the KH gradient. Results (presented below) generally fit reasonably well with our observations of patterns in biogeochemical dynamics along that gradient. The main areas in which we used observations from the KH gradient to determine the parameter values used in the model were in the identification of leaching losses of nitrate rather than denitrification as the dominant pathway of nitrate loss in wetter sites (discussed above) and in the distribution of biological N fixers in the driest sites. Our observations showed that symbiotic N-fixers were sparse and relatively inactive (in fixing N) in these driest sites (Burnett 2019; Burnett et al. submitted for publication), but initial simulations yielded a greater-than-observed abundance of N-fixers there. The difference may reflect a lack of introduced pasture legumes that are adapted to very dry sites, and the elimination of native N-fixers by over a century of grazing; in other words, the result may be specific to this gradient rather than a general feature of climate gradients. In any case, to simulate KH we reduced the NPP of N-fixers and so the rate of N inputs via biological N fixation in dry sites in accordance with these observations. Below we will present results both with and without this reduction in N fixation in dry sites.

Isotopes of N

The modeled and observed (the latter from von Sperber et al. 2017) distribution of δ¹⁵N along the KH gradient are shown in Fig. 6. Observations here (and throughout) are of the 0–30 cm depth in the soil, while the model implicitly calculates the 0–50 cm depth layer in the soil. Because δ¹⁵N integrates multiple processes which can fractionate it in different directions, it is not easy to model correctly—but the agreement in pattern between results of the toy model and field observations is striking. Observations and simulation results differ most in wetter sites on the gradient, where observations yield soil δ¹⁵N that is a few ‰ enriched in ¹⁵N, while simulations yield little enrichment (or slight depletion) in δ¹⁵N. The difference could be due to an underestimate of denitrification (and its associated fractionation of the lost N), or to partial preservation of ¹⁵N-enriched organic matter on short-range ordered clay minerals in wetter sites (Kramer et al. 2012). Both observations and simulations then show an increase in enrichment in drier sites to a peak near + 14‰, then a decline in enrichment below ~ 3.5 cm/month precipitation (− 1400 mm/year water balance). This decline in enrichment from moderately dry sites into the driest sites also is observed for ¹⁵N in nitrate (von Sperber et al. 2017) and for ¹⁵N in the tissues of non-N-fixing plants (Burnett 2019; Burnett et al. submitted) on this gradient.

Fig. 6

Observed (points) and modeled (line) δ¹⁵N on the Kohala (KH) climate gradient. Observations from von Sperber et al. (2017). In this and subsequent figures, we label the lower X-axis as Mapped and Modeled Mean Monthly Precipitation (in cm/month) where both observations and model results based on stochastic month-to-month variation in precipitation are reported, as Modeled Mean Monthly Precipitation when only model results with stochastic variation in precipitation are reported, as Mapped Monthly Precipitation where only empirical results are presented, and as Modeled Monthly Precipitation when only results from a simulation that did not include stochastic month-to-month are presented. The upper X-axis shows annual water balance (precipitation (ppt) minus Potential Evapo Transpiration (PET)) in mm/year; it follows the same conventions for classes of results as does the lower x-axis. Mapped precipitation and water balance results from Giambelluca et al. (2013)

The model includes stochastic variation in precipitation; accordingly, we summarize the reproducibility of δ¹⁵N for multiple model runs in Table 1. Also, the model includes multiple constants for which we assumed reasonable values but for which there are alternative reasonable values; accordingly, we summarize the sensitivity of simulated δ¹⁵N to these assumptions in Fig. 7. We focus these explorations on δ¹⁵N because as discussed in the Introduction, δ¹⁵N is an integrated measure of processes involved in N dynamics, and we are particularly interested in N cycling and its controls here.

Table 1 The variation in δ¹⁵N observed among 10 independent runs of the model

Fig. 7

A sensitivity analysis for δ¹⁵N in the toy model. Results for δ¹⁵N across the KH climate gradient under standard conditions are shown with a solid line (this is the mean of 10 runs of the model, reported in Table 1); the consequences of lower values of some key variables that could affect δ¹⁵N are displayed with various permutations of thin dashed lines, while the effects of higher values of those variables are shown with thick dark gray dashed lines. X axis as in Fig. 6. (Color figure online)

As Table 1 demonstrates, repeated model runs yield similar results, with some greater (but still small) variation in lower-rainfall sites where rainfall variability (the source of stochasticity in the model) is greater (Table 1). This result is unsurprising given that the results presented (Fig. 6, Table 1) represent mean values for the last 10,000 time steps (months) of simulation. The sensitivity of δ¹⁵N to differing assumptions concerning parameter values (Fig. 7) is also relatively small; δ¹⁵N follows the same pattern across the gradient for a wide range of assumptions about the values of driving variables—with the exceptions that assumptions about δ¹⁵N in inputs and losses of N can have substantial effects on the results (Fig. 7), as can the empirical observation that there are few and relatively inactive N fixers in the driest sites (Fig. 13a).

Other model/observation comparisons

Modeled and observed variations in Ca and P along the KH gradient are summarized in Fig. 8. While the model evaluates a generic cation that is more or less based on Calcium (Ca), the observations here are of the Ca (and P) remaining in the soil relative to the amount that was present initially in parent material [elements remaining were calculated using Niobium (Nb) as an immobile index element (Brimhall et al. 1992; Kurtz et al. 2000)]. Observations and simulations are not precisely parallel, in that the observations account for the gain or loss in mass of the parent material due to additions to soil or depletion of labile minerals and so for the expansion of soil volume with the addition of organic matter and hydration, or collapse in soil volume as material is lost, while the simulations deal with a fixed depth interval. Nevertheless, for P the correspondence between model and observations is good, perhaps because there is little net expansion or collapse of soils in this range of substrate age. In both simulation and observations, P concentrations are consistent and intermediate (neither enriched nor depleted) in the driest sites. There is then a substantial enrichment in P in intermediate-rainfall sites (from 7 to 15 cm/month precipitation, or − 1000 to 0 mm/year water balance, for the simulation; from 11 to 15 cm/month precipitation or − 500 to 0 mm/year water balance for observations) This peak in P reflects biological uplift (Jobbagy and Jackson 2001; Vitousek and Chadwick 2013) (uptake of water and P from the deep soil layer, deposition of P in litter on the soil surface); it defines the location of an intensive indigenous agricultural system that occupied this portion of the climate gradient, in the centuries before European contact (Vitousek et al. 2004; Marshall et al. 2017), which drew upon this uplifted P. P then declines monotonically to low concentrations in wetter sites, in both observations and simulations.

Fig. 8

Observed (points) and modeled (lines) for P (solid line, filled symbols) and Ca (dashed line, hollow symbols) remaining within surface soils on the KH climate gradient. Observations from Vitousek and Chadwick (2013); observed % of these elements remaining is integrated to a soil depth of 30 cm and is calculated by comparison with parent material that contributes to this surface soil, determined using Nb as an immobile index element. Modeled % remaining is to a depth of 50 cm and is calculated by comparing the total P and Ca (organic and mineral forms for P) remaining in soil with the initial primary mineral pool in that depth interval. Values > 100% for P remaining reflect biological uplift of P from deeper soils (Vitousek and Chadwick 2013). X axes as in Fig. 6

The correspondence between simulation and observations is not as strong for Ca, except that both model and observations show Ca to be lost more readily than P, and both model and observations both show Ca to be strongly depleted in the highest-rainfall sites (more so for Ca than for P). In dry sites (< 7 cm/month precipitation or <  − 800 mm/year water balance) observed Ca is less than simulated Ca; possible causes of this difference are the shorter time of simulations relative to soil development (but see below on time courses), or aeolian transport of already-weathered soils from upslope. In intermediate rainfall sites, there is a sharp threshold where Ca concentrations decline rapidly at ~ 7 cm/month precipitation in the model or ~ 14 cm/month precipitation (~ − 200 mm/year water balance) in observations (Fig. 8). Here we suspect the difference between observations and simulation is that we did not simulate biological uplift of Ca, but it occurs in practice, as was demonstrated by Bullen and Chadwick (2016). Observations and simulations are similar in showing 3 soil process domains (bounded by two pedogenic thresholds) on the gradient (Vitousek and Chadwick 2013); one threshold occurs in moderately dry sites where rock-derived nutrients remain at rainfall levels below the threshold and are depleted at higher rainfall (Chadwick and Chorover 2001); the other occurs in wetter sites at the upper-rainfall bound of sites where the biological uplift of P is important (Vitousek and Chadwick 2013).

Model results and observations for soil organic N are shown in Fig. 9. Model and observations are not directly comparable for soil N; not only do simulation and observations reflect different soil depths, as described above, but also the observations are of the concentration of organic N (in % dry mass) while the simulation is of the N pool in mass/area (here kg/ha). We think the difference in depth is relatively unimportant, but the comparison of N concentrations and N inventories could be more important. Nevertheless, the correspondence in pattern along the KH gradient for simulations and observations is relatively strong (Fig. 9). Both show low N in the driest sites, an increase to a peak near 20 cm/month precipitation (~ + 600 mm/year water balance), and then a decline into the wettest sites. We note that Fig. 9 excludes two points with very high N concentrations in relatively wet sites; these high concentrations are outliers relative to the overall pattern (Von Sperber et al. 2017).

Fig. 9

Observed (points) and modeled (line) soil N on the KH climate gradient. The left y-axis is the modeled organic N pool on an areal basis (kg/ha); the right y-axis is the soil N concentration observed in each site on the gradient (from von Sperber et al. 2017). Most of the observed N concentration is organic N, although any ammonium-N present in the soil would be included in this pool. X axes as in Fig. 6

Simulations of unmeasured fluxes

In addition to properties for which we can compare observations to simulations (even with caveats), we can simulate a number of unmeasured properties. For example, we model N losses by multiple pathways (Fig. 10a) with interesting results: low-rainfall sites (from 3 to 7 cm/month ppt, or − 1200 to − 1400 mm/year water balance) lose most of their N via ammonia volatilization; intermediate-rainfall sites from 7 to 12 cm/month precipitation (or ~ − 1200 to ~ − 300 mm/year water balance) loses most N as nitrate; and wetter sites lose most N as DON. (Transformation-dependent trace gas fluxes peak in the domain where nitrate losses are most important, but their magnitude is at most < 5% of that of nitrate losses.) These differing pathways of N loss shape the pattern in δ¹⁵N (Fig. 6), they are also important in that while ammonia volatilization and nitrate loss both remove an accumulated pool of biologically available N that could be retained by organisms, DON losses represent uncontrollable losses of N and thus meet one of the two requirements for substantial and/or persistent N limitation to biological activity, as discussed below (Vitousek and Field 1999).

Fig. 10

Modeled features for which there are no field observations along the KH climate gradient. a Pathways of N loss. The thinner black solid line on the left represents ammonia volatilization, the thick dark gray line in the center represents nitrate loss, and the thinner black dashed line to the right represents Dissolved Organic Nitrogen (DON) loss. As described in the text, we treat nitrate loss as representing nitrate leaching and not denitrification. b Simulated uptake of water from a deep soil layer, calculated as described in the text. c Simulated P fractions in the surface soil (to 50 cm depth); primary mineral P represents primary minerals that have not yet broken down, labile P is soluble, biologically available P that corresponds to the PO₄ pool in Fig. 3, and secondary mineral P is mineral P that formed in the soil. X axes as in Fig. 6. (Color figure online)

Other simulation results of interest include the uptake of water from a deep soil layer (Fig. 10b). We found that in dry sites, there was insufficient deep soil water to support transpiration from the 0.5 to 2 m depth in soil; in wet sites, all of potential evapotranspiration could be supported from the surface horizon, so again, there was no uptake of water from deep soil. However, at intermediate rainfall (in the simulation peaking at ~ 14 cm/month ppt, or ~  − 200 mm/year water balance), water uptake from deep soil was an important fraction, up to almost 20%, of total transpiration. Most likely, most of this deep-water uptake occurred during the many millennia when deep-rooted mesic native forests occupied the middle portion of the gradient (Chadwick et al. 2007)—before an indigenous cropping system was established here, even longer before pastures dominated the landscape.

We could also estimate soil P fractions from the toy model (Fig. 10c). In this simulation, after 60,000 time steps (months) of simulation, by far the most abundant form of soil P in surface soils in all but the wettest sites was secondary mineral P. This result was interesting in that we assumed (and utilized coefficients reflecting) that P-containing primary minerals would weather more slowly than would Ca-containing primary minerals. However, by having half of the available P in the soil go into secondary mineral forms, we reduced concentrations of labile mineral P substantially and so enhanced weathering of P-containing primary minerals by reducing end-product inhibition of their weathering. The end result was that P minerals broke down more rapidly than did Ca minerals, even though the kinetic and equilibrium weathering coefficients of P minerals were lower. We believe this is a meaningful result.

As a comparison of Fig. 10b and c shows, the simulated peak in total and secondary mineral P occurred at a lower rainfall than did the simulated maximum in deep soil water uptake. We believe that while more P is taken up and deposited as litter on the soil surface where deep water uptake is at its maximum, more of this P leaches back to deep soils where rainfall is greater.

Precipitation variability

We also used the toy model to evaluate the implications of month-to-month variation in precipitation by comparing simulated soil and ecosystem properties with month-to-month variation in precipitation (as described in “Methods” section) versus with a constant precipitation each month (~ equivalent to the mean of the variable precipitation at each point on the gradient (Fig. 11). Month-to-month variation in precipitation means that some months have much greater than average precipitation, and so leaching of nutrients occurs at levels of (average) rainfall below those where it would occur with consistent rainfall with the same mean. Because such high-rainfall months occur infrequently (and because leaching losses may require consecutive months of above-average rainfall), rock-derived elements continue to weather and be lost over long time scales (versus equilibrating quickly to constant rainfall) (Vitousek et al. 2016). Also, because of the way N and P mineralization are modeled here (above), variable rainfall can lead to asynchrony between element supply via mineralization and plant demand for elements for growth (Vitousek and Field 2001). Results of these simulations showed that a realistic degree of monthly variation in precipitation has a substantial effect on simulated ecosystem properties, including δ¹⁵N (Fig. 11a), patterns of variation in soil P (especially the peak in total P in surface soils at intermediate rainfall) (Fig. 11b), and plant biomass (Fig. 11c). For δ¹⁵N and soil P, the simulations with month-to-month variation in precipitation yielded results that were much closer to observations than were simulations without such variation (Figs. 6, 8 and 11a, b); there were no observations for plant biomass. Without month-to-month variation in precipitation, there was no leaching loss of P or Ca at rainfall below the point where water balance = 0. That lack of losses persisted indefinitely, while with month-to-month variation there is a continuing depletion in rock-derived elements over time (Vitousek et al. 2016). Without variation in precipitation, we simulated some increase in P and Ca in dry sites due to atmospheric deposition and a lack of leaching, but that increase was much less than the simulated biological uplift of P (Figs. 8, 11b).

Fig. 11

The simulated consequences of removing month-to-month variation in precipitation along the KH gradient. Dashed lines show results of simulations with stochastic month-to-month variation in precipitation included, solid lines indicate results of simulations with constant precipitation each month. The long-term average of monthly precipitation is approximately the same as the constant precipitation, at each point on the climate gradient. a δ¹⁵N, as in Fig. 6. b P and Ca remaining from parent material, as in Fig. 7. c Plant biomass. X axes as in Fig. 6. (Color figure online)

Nutrient limitation

With the toy model, we can evaluate simulated nutrient limitation to plant growth, by adding N and/or P to our simulations and evaluating plant response. For this test, we run the simulation for 45,000 time steps (with month-to-month random variation in precipitation after 12,000 time steps), then add N and/or P at each time step for the last 15,000 time steps (averaging and presenting results for the last 10,000 time steps as above). Since plant biomass (in the model) is a multiple-month integration of net primary production (NPP), biomass and longer-term average NPP can be regarded as synonymous. This simulation yields long-term rather than immediate responses to nutrient additions—for example, if P limits the growth of and N fixation by N fixers and N limits the growth of the non-fixer, the time scale of our simulation means that modeled biomass of the non-fixer in our simulations will respond to both N and P additions (since P stimulates N fixation and that in turn will increase N supply at this long time scale). These dynamics are illustrated in Fig. 12, in which simulated additions of P cause a rapid increase in the growth of the N-fixer; the consequent increase in N fixation drives a slower but sustained increase in growth of the non-fixer.

Fig. 12

The time course of responses in non-N-fixer and N-fixer plant biomass C to P additions within the toy model at the highest rainfall and most positive level of water balance. P additions were simulated each month beginning in month 45,000; the N fixer responded quickly to added P with increased growth (as described in the text, biomass of both plants is a many-month integration of primary productivity) and increased N fixation; this added N drove a slower increase in growth of the non-fixer. Growth of the non-fixer in turn partially suppressed the growth of the N-fixer, because of the non-fixer’s modeled priority in resource use. (Color figure online)

We are particularly interested in understanding what the toy model can tell us of the drivers of simulated N limitation here, because we regard the widespread importance of N limitation to be one of the greatest unsolved mysteries in biogeochemistry. Organisms with the capacity to fix N biologically are ubiquitous, and symbiotic N fixers in particular have the capacity to quickly overcome any deficit in N (as well as possessing what would appear to be a large competitive advantage in N-limited ecosystems). We see the co-occurrence of two conditions as being essential to N limitation that is more than marginal and/or ephemeral. First, there must be constraints to N fixation that keep it from responding fully to N deficiency. Second, there must be losses of N from the system that cannot be prevented by organisms within the system, even when they are limited by N. The first condition is widely appreciated (Vitousek and Howarth 1991; Menge et al. 2008); the second condition is important because even low rates of N inputs (such as atmospheric deposition of fixed N in unpolluted regions) would otherwise eventually accumulate to the point that even the slow turnover of a large pool of N would suffice to meet the needs of organisms and reverse N limitation.

Biomass of the non-N-fixing plant, the N-fixing plant, and the sum of both without nutrient additions, with N additions, with P additions, and with N + P additions in the toy model are summarized in Fig. 13. Month-to-month variation in precipitation causes asynchrony between the release of N and P by mineralization and their uptake by plants; consequently, losses of both N and P via leaching can occur across a wide range of precipitation. These losses reduce biomass of the non-fixer, and cause the simulated non-fixer to have its growth ultimately limited by P at all levels of precipitation (Fig. 13a) (although this long-term P limitation to the non-fixer often is a consequence of P limitation to N fixation—Fig. 12). As a consequence of P limitation to the N-fixer, the non-fixer is limited by N in all sites with simulated average precipitation of < 4 cm/month (~ − 1200 mm/year water balance) and > 12 cm/month (water balance > − 300 mm/year) (Fig. 12a), as is demonstrated by the non-fixer’s positive response to added N in those sites. The N-fixer is limited by P supply everywhere on the gradient (Fig. 12b); in the model, its growth and biomass are suppressed by that of the non-fixer when and where climatic conditions are favorable to the non-fixer and where added N (from N fixation or fertilization) enhances growth of the non-fixer (Figs. 12, 13a, b). We interpret these results to mean that (1) growth and activity of the N fixer is constrained by P supply across the gradient; (2) in the highest-rainfall sites, losses of DON represent uncontrollable losses of N; together with N fixation that is constrained by P supply, this means that the necessary conditions for N limitation are met); and (3) in intermediate-rainfall sites (12–22 cm/month simulated precipitation, or − 300 to + 1400 mm/year simulated water balance) the asynchrony between N mineralization and demand for N by the non-N-fixing plant drives losses of N. While these losses are of nitrate, a form of N that organisms can use readily, they occur at times when the abundance of biologically-available N is greater than is demand for that N, so they too represent uncontrollable losses of N that, together with constraints to N fixation (here again constrained by the supply of P across the gradient) (Fig. 13b), meet the necessary conditions for N limitation. This is true because losses of nitrate that occur when its supply is greater than demand reduce N pools, and cause (at other times) the transient situation in which demand is greater than supply and as a consequence N limits production of the non-fixer. Similar dynamics cause N limitation in the driest sites (in this case involving losses of N via ammonia volatilization (Fig. 10a).

Fig. 13

Simulated long-term N and P limitation (with N and P limitation determined independently and interactively, by adding N and/or P at each time step from time step 45,000 through the end of the simulation). a Non-N-fixing plant. b Symbiotic biological N-fixing plant. c Sum of non-fixing and N-fixing plant. X axes as in Fig. 6. (Color figure online)

Model caveats

In Fig. 14 we address two caveats to generalization of results of the toy model. First, the most consequential (to the integrative measure δ¹⁵N) empirical adjustment to the toy model is the incorporation of the observation that symbiotic biological N fixing plants are sparse to absent in the driest sites. The observation that shrubby or herbaceous symbiotic N-fixers are absent in the driest sites describes the current state of the KH gradient accurately (Burnett 2019, Burnett et al. submitted), and while the phreatophic leguminous tree Prosopis pallida is present in the driest sites, we avoided the increasingly scattered individuals of that tree in increasingly drier sites in our sampling. The lack of shrubby or herbaceous N-fixers in dry sites at KH may reflect the paucity of native N fixers in Hawai‘i, the vulnerability of those that were present to grazing, and the lack of introduced pasture legumes in the driest sites. To test the implications of the possibility that the observed lack of N fixers in the driest portion of the KH gradient is not a general feature of climate gradients, we ran the toy model with and without any systematic constraint to biological N fixation in the driest sites (Figs. 6, 14a); we observed thereby that without the empirically-based constraint of low N fixation in the driest sites, δ¹⁵N did not decrease sharply in the driest sites, unlike field observations. The solid line in Fig. 14a may thus reflect a more general pattern than the KH-specific result summarized by the dashed line in Fig. 14a, though it should be noted that a climate gradient across North China (Wang et al. 2014) yielded a peak in δ¹⁵N at the same point (though of a lower magnitude) than the one that was observed on the KH gradient. We know of no reason to think that symbiotic N fixation is constrained by a low diversity of legumes in China (unlike in Hawai‘i where legumes may well be so constrained).

Fig. 14

Adjustments to the standard conditions for model runs for the KH gradient. a Testing the implications of an empirically-based reduction in N fixation in the driest sites. The dashed line shows simulated δ¹⁵N with this reduction in N fixation in place, the solid line shows simulated δ¹⁵N with that reduction removed. b Testing the consequences of a greater number of time steps in a simulation (600,000 time steps, versus the standard of 60,000 time steps). Most soil and ecosystem properties were little affected (directly) by the longer model run; however, patterns in rock derived elements (shown here for P and Ca) differed, in that month-to-month variation in precipitation led to continued cumulative losses of rock-derived elements. Solid lines show results for 600,000 time steps, dashed lines for 60,000 time steps. X axes as in Fig. 6. (Color figure online)

Second, we evaluated the consequences of extending the simulations, comparing the standard run of 60,000 time steps with a run of 600,000 time steps (again introducing month-to-month random variation in precipitation after 12,000 time steps, and averaging and presenting results from the last 10,000 time steps). We observed little direct effect of the longer run on most soil and ecosystem properties; the major exception was for total P and Ca pools in the soil (Fig. 14b). There was a shift in the depletion of rock-derived elements with the longer simulation towards drier sites for both P and Ca, and an increase in magnitude and a shift to a drier position on the gradient for the peak in P (Fig. 14b). The shift in the depletion of weathering-derived elements towards drier sites reflects the accumulation of rare leaching events in dry sites within the simulation and consequently the fact that the average rainfall where the pedogenic threshold that represents the depletion of nutrients shifts slowly towards drier sites over many time steps (Vitousek et al. 2016), as has been observed empirically on older substrates (Vitousek et al. 2014). The shift in position of the peak in P represents a similar process (cumulative leaching of P from surface soils in wetter sites); the increase in magnitude of the peak reflects the longer time over which P from deep in the soil has been lifted up to the surface. The fact that with 10× the number of time steps the magnitude of the peak has only doubled shows that the increase does not proceed linearly or indefinitely; the fact that the enrichment in dry sites (resulting from atmospheric deposition where there is no leaching) is much less than the peak in P (Fig. 14b) shows that (in the simulation) the peak is due to biological uplift of P from deep in the soil more than it is due to cumulative atmospheric deposition of P.

Temperature effects on biogeochemistry

We further used the toy model to explore the importance of temperature control versus water (precipitation or water balance calculated as ppt-PET) across the KH climate gradient (Fig. 15), where temperature and precipitation are so highly correlated (r > 0.99) that it is impossible to disentangle controls empirically. We tested the importance of temperature on the KH gradient by running the toy model with and without an influence of temperature incorporated—with both temperature and soil moisture influencing the scalars that are used to calculate weathering, NPP, and decomposition/mineralization in the model, and with just soil moisture influencing those scalars. The results of this test (Fig. 15) are clear; soil water controls most of the patterns in biogeochemical dynamics along the KH gradient. Temperature does have an influence on plant biomass/NPP in the cooler, wetter, higher-elevation end of the gradient (Fig. 15d), and consequently an influence on other biogeochemical processes there. However, the presence and position of high and low points in (and so the pattern of) biogeochemical features across the gradient are the same with and without the influence of temperature (Fig. 15).

Fig. 15

Consequences of simulating biogeochemical patterns on the KH gradient with (dashed lines) and without (solid lines) the incorporating the influence of temperature, while retaining an influence of soil water. a Simulated δ¹⁵N in ‰. b Simulated P and Ca pools, in % of the initial mineral pool remaining in the top 50 cm of soil. c Simulated soil organic N pool, in kg/ha. d Simulated biomass of non-fixing and N-fixing plants, in kg of C/ha. While some pools (notably plant biomass) and fluxes vary in the presence or absence of an influence of temperature, patterns of variation in biogeochemistry along the KH gradient are not influenced substantially by temperature in these simulations. X axes as in Fig. 6. (Color figure online)

We also attempted to extend the toy model to the Mauna Kea (MK) climate gradient that is described in Bateman et al. (2019), motivated in part by the very different pattern of variation in δ¹⁵N as a function of precipitation (or water balance) that we observed on that gradient compared to KH (Fig. 16). While δ¹⁵N on the KH peaked at ~ 3.5 cm/month precipitation (~ − 1400 mm/year water balance), δ¹⁵N on the MK gradient peaked at a much lower enrichment of ~  + 8‰ (versus + 14‰ for KH) at ~ 19 cm/month precipitation (~ + 1000 mm/year water balance) (Fig. 16). Because the MK gradient has cooler temperatures throughout than does the KH gradient, and because the driest sites on MK have the coolest temperature of all (mean annual temperatures of < 6 °C) (Giambelluca et al. 2013), versus the driest sites being the warmest on KH, we expected that the MK gradient would provide an opportunity to evaluate the effects of temperature on nutrient dynamics along a climate gradient. It is worth noting that while PET generally decreases at low temperatures, the situation on the MK gradient is one in which the cold dry sites are located above a trade wind temperature inversion, in a region where humidity is very low (Giambelluca and Nullet 1991). The low humidity offsets the low temperature there, so that PET varies little across the whole of the MK gradient.

Fig. 16

Observations of δ¹⁵N along the KH gradient (solid symbols) and the MK gradient (hollow symbols). δ¹⁵N integrates multiple sources of N input/output and multiple transformations in the N cycle. Results for KH from von Sperber et al. (2017); δ¹⁵N for MK is reported for the first time here. X axes as in Fig. 6

Previous comparisons of MK to KH have emphasized their difference in substrate age (~ 20 kyr for MK, versus ~ 150 kyr for KH) (Bateman et al. 2019). However, we could use the toy model to evaluate the potential importance of climate vs substrate age in causing differences between these gradients. To test the role of time directly, we compared the consequences of running the model for the standard 60,000 time steps versus running it for 600,000 time steps (Fig. 14). As described above, the additional time did not affect δ¹⁵N (or other biological processes involving N) substantially. However, the additional time did drive depletion of rock-derived cations and P to lower levels precipitation (and water balance) (Fig. 14b). We conclude that while differences between KH and MK in the weathering and loss of rock-derived elements could reflect differences in substrate age, differences in N probably reflect differences in climate.

In practice, the toy model did not simulate biogeochemical dynamics on the MK gradient very well. We applied the model formulation that worked reasonably well for KH to the MK gradient, adjusting climate to that of the MK gradient, and compared results of that simulation with one in which the influence of temperature was not included (Fig. 17). Neither including nor excluding the effects of temperature resulted in a close match to observations on the MK gradient (Figs. 16, 17), except for P with incorporation of a temperature effect (Fig. 17b). The strong conclusion we could reach is that low temperatures on the MK gradient are capable of changing simulated biogeochemical patterns, not just the peak rates of some biological processes (Fig. 17); this effect was particularly strong for plant biomass. (Fig. 17d).

Fig. 17

Consequences of applying the model used to simulate biogeochemical dynamics on KH to the MK gradient, changing only the range of precipitation/water balance values existing across the gradient and the range of temperatures to those of the MK gradient. Simulations that include the influence of temperature are shown with dashed lines, those that exclude the influence of temperature are shown with solid lines. a Simulated δ¹⁵N in ‰ (observed δ¹⁵N is reported in Fig. 16). b Simulated and observed P and Ca remaining on the MK gradient, in % of the initial mineral pool remaining in the top 50 cm of soil, and observed Ca and P remaining, in % of the pools in parent material in the top 30 cm of soil (both calculated as in Fig. 8). c Simulated soil organic N pool in the upper 50 cm of soil, in kg/ha, and N concentration in the upper 30 cm of soil in % dry mass. Observations in b and c from Bateman et al. (2019). d Simulated biomass of non-fixing and N-fixing plants, in kg of C/ha. For most biogeochemical properties, simulations for MK do not do as well as those for KH in reproducing observations, where observations are available. The presence/absence of the influence of temperature in the range observed on the MK gradient can change the simulated pattern, not just the simulated magnitude, of biogeochemical variation on that gradient. X axes as in Fig. 6. (Color figure online)

Conclusions

A simple model of variation in biogeochemical processes on climate gradients provided insight into the control of biogeochemical processes on the KH climate gradient. Because the model fit aspects of N cycling on the KH gradient, including δ¹⁵N, particularly well (Figs. 6, 9), it was useful for evaluating the conditions under which N limitation to primary production could take place (Fig. 13). N availability limited the simulated NPP and biomass of a non-N-fixer where (in the model) there was a constraint to the ability of the simulated N-fixer to respond fully to N deficiency, and where there were losses of N from the system that could not be prevented by biota in the system. In the simulation, N fixation was constrained by P supply across the gradient, in part because the inclusion of a secondary mineral pool of P in the simulation depleted primary mineral P rapidly. The pathway of the uncontrollable N loss shifted from losses of ammonia or nitrate due to asynchrony between plant demand for N and the supply of N in dry and intermediate sites on the gradient, to leaching losses of DON in the wettest sites (Figs. 10a, 13). In the model, the asynchrony between N supply and demand in dry and some intermediate-rainfall sites on the gradient was caused by month-to-month variation in precipitation, which meant that N was in excess of organisms’ requirements at some times and its losses at those times when supply exceeded demand constrained plant growth at times when demand exceeded supply.

We also used the toy model to evaluate whether temperature or rainfall (or water balance calculated as ppt-PET) had a greater effect on biogeochemical patterns. Temperature and rainfall were highly correlated on both gradients, but the relationship ran in opposite directions, with the driest sites being warm low-elevation leeward sites on KH and the driest sites being cold high elevation alpine sites on MK. Results demonstrated that the KH gradient is primarily a gradient in rainfall (and water balance); removing any influence of temperature on biogeochemical processes there did not affect biogeochemical patterns (Fig. 15), though simulated temperature did affect the magnitude of some pools and fluxes there. In contrast, the inadequacy of the toy model’s application to MK precludes strong conclusions about temperature vs rainfall effects there, but it is clear that temperatures in the range of those observed on the MK gradient are sufficient to affect the pattern, not just the magnitude, of biogeochemical variation there (Fig. 17).