KEYLINK: towards a more integrative soil representation for inclusion in ecosystem scale models—II: model description, implementation and testing

Structural effects

Pore size distribution determines accessibility for trophic interactions of soil fauna and soil microorganisms (Fig. 2), both by size and by aeration and H₂O; soil fauna changes pore size distribution and produces cracks and fissures in the soil. In the model, pore size distribution is divided into the following five categories:

• Inaccessible pores (<0.1 µm in diameter): pores around inaccessible C (within the micro-aggregate, organo-clay interaction). Water is held here but is not available to plants (measured from wilting point). The volume of inaccessible pores is related to the clay content and type.

• Bacterial pores (0.1–2 µm): the pores within macro-aggregates and pores in loam, accessible only to bacteria. Engineer saprotrophs (e.g., earthworms) can also use SOM in these pores (and in the following pore categories, all except inaccessible pores) because they eat directly all soil.

• Micropores (2–30 µm): pores not accessible to macrofauna, mesofauna and most predators, but accessible to microfauna bacterivores and fungivores, fungi, mycorrhiza and bacteria. Water is held at field capacity but available to plants. In sandy soil and within macro-aggregates (>250 µm), pores fall in this category.

• Mesopores (30 µm–1.5 mm): pores where most soil fauna can penetrate (not macrofauna) between large macro-aggregates (>1 mm) or formed by fine roots. Mesopore volume can be determined in the field from drained water capacity (but this includes macropores). These pores are well aerated also at field capacity, but can dry out below field capacity.

• Macropores (>1.5 mm): cracks or biopores formed by ecosystem engineers. They are of vital importance for soil hydrology as preferential flow through these pores has a major impact on infiltration rate. These are the first pores to have O₂ when water level is above field capacity, but dry out quickly below field capacity.

The initial values of soil porosity in the model simulations can be calculated from measured soil water retention curves, or even using models such as Saxton et al. (1986) that yield field capacity, porosity and wilting point from the C, clay and sand contents, or using measured bulk density (D_b). Following Malamoud et al. (2009), the percentage of total porosity (P_%) can be computed from D_b and soil particle density (D_S) as shown in Eq. (2). D_S can be measured or is calculated from D_m = soil mineral particle density (2.65 g cm⁻³) and D_SOM = organic particle density (1.35 g cm⁻³) as: (1) $$D_{\mathrm{s}}={\displaystyle \frac{100}{{\displaystyle \frac{\mathrm{\%}\mathrm{S}\mathrm{O}\mathrm{M}}{D_{\mathrm{S}\mathrm{O}\mathrm{M}}}+{\displaystyle \frac{100-\mathrm{\%}\mathrm{S}\mathrm{O}\mathrm{M}}{D_{\mathrm{m}}}}}}}$$

(2) $$P_{\mathrm{\%}}={\displaystyle \frac{D_{\mathrm{s}}-D_{\mathrm{b}}}{D_{\mathrm{s}}}100}$$

Water flow

We advise using KEYLINK model in combination with a detailed water model including preferential flow through macropores as well as good representation for matrix flow (s.a. Richards’ equation). However, we show in this paper how it can be used with a simpler representation of water flow but still allowing the important dynamic interactions between pore sizes and hydrology that are fundamental to the model. A spilling bucket approach is used at a daily time-step, where water drains from a layer into the underlying layer when its water content is above field capacity in the soil matrix. However, in contrast to conventional spilling bucket models, we allow water to flow faster through macropores (before the soil matrix is saturated). Net precipitation (P_net) is calculated as: (3) $$P_{\mathrm{n}\mathrm{e}\mathrm{t}}=P-E$$where P is precipitation (mm day⁻¹) and E is evapotranspiration (mm day⁻¹) from measured or modelled data (vegetation model). Infiltration (I) is assumed to be equal to the part of precipitation entering the soil. Infiltration and runoff (P_runoff, mm day⁻¹) must equal P_net.

(4) $$I+P_{\mathrm{r}\mathrm{u}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{f}}=P_{\mathrm{n}\mathrm{e}\mathrm{t}}$$Infiltration is composed of water entering the soil matrix, water filling the macropores and water draining from macropores. Water that enters macropores remains in the macropore domain or enters the layers below. The fraction of infiltration entering macropores depends on the surface area of the macropores (SA_macro), assumed cylindrical. Assume measured or derived maximal infiltration rate (I_maxMat, mm day⁻¹) of the soil matrix. Maximal infiltration rate through macropores (I_maxPor, mm day⁻¹) is calculated from the volume of the pores (PV_macro), assumed not limiting at daily scale, plus infiltration capacity of the layer (n + 1) in which the macropores end.

(5) $$I_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{P}\mathrm{o}\mathrm{r}}=\mathrm{P}{\mathrm{V}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}+I_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{M}\mathrm{a}\mathrm{t}\left(\mathrm{n}+1\right)}$$If P_net > (I_maxPor + I_maxMat) runoff is calculated as: (6) $$P_{\mathrm{r}\mathrm{u}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{f}}=P_{\mathrm{n}\mathrm{e}\mathrm{t}}-\left(I_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{P}\mathrm{o}\mathrm{r}}+I_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{M}\mathrm{a}\mathrm{t}}\right)$$after which calculations continue using P_net − P_runoff as net precipitation.

If I_maxMat < P_net < (I_maxPor + I_maxMat) the soil matrix is filled at a rate equal to the maximum infiltration rate, all other water is lost either through the macropores to the next layer or by filling macropores. If I_maxMat > P_net the soil matrix is filled with water, traditional spilling bucket, but an equivalent volume is lost through macropores to the bottom layer depending on the surface area of the macropores. The total soil water volume of soil layer n, SW_n, is then limited by the total pore volume of the layer and the water already filling the pores, and is calculated as: (7) $${\mathrm{S}\mathrm{W}}_{\mathrm{n}}=\mathrm{S}{\mathrm{W}}_{\mathrm{n}}+\mathrm{m}\mathrm{i}\mathrm{n}\left(\mathrm{P}{\mathrm{V}}_{\mathrm{n}}-\mathrm{S}{\mathrm{W}}_{\mathrm{n}},I_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{m}\mathrm{a}\mathrm{t}}\left(1-\mathrm{S}{\mathrm{A}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}\right),P_{\mathrm{n}\mathrm{e}\mathrm{t}}\left(1-\mathrm{S}{\mathrm{A}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}\right)\right)$$For drainage (D) to the bottom layer, the spilling bucket approach is used plus a portion of water that goes straight through the macropores, calculated from the surface area of the macropores.

(8) $$D_{\mathrm{n}}=P_{\mathrm{n}\mathrm{e}\mathrm{t}}\mathrm{S}{\mathrm{A}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}+P_{\mathrm{n}\mathrm{e}\mathrm{t}}-\mathrm{m}\mathrm{i}\mathrm{n}\left(\mathrm{P}{\mathrm{V}}_{\mathrm{n}}-\mathrm{S}{\mathrm{W}}_{\mathrm{n}},I_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{m}\mathrm{a}\mathrm{t}}\left(1-\mathrm{S}{\mathrm{A}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}\right),P_{\mathrm{n}\mathrm{e}\mathrm{t}}\left(1-\mathrm{S}{\mathrm{A}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}\right)\right)$$For each pore size class the fraction water filled is calculated from the water content: so always one pore size is partially saturated and all others are either saturated or dry within one layer.

C flow

The KEYLINK model combines SOM modelling with soil food web modelling. The model conceptualized in Fig. 2 has 13 carbon pools (Table S1.2), visualised by boxes. Above and belowground litter is assumed to be provided from an external source (tree shoot in Fig. 2) not covered by this model. It could be given through experimental data or an external model, for example, a tree growth model that delivers the input of litter into the litter pool. All simulations presented here were made with constant C inputs (Table S2.6). Exudation is an input of organic carbon released from roots into the SOM pool. Every live pool has a respiration rate (r) and a turnover or death rate (d). On consuming a C pool, a fraction of this pool always becomes faeces and enters the SOM pool except for the microbial pools, that is, microbes and microbivores. SOM can be distributed in different fractions, particulate organic matter (POM) and dissolved organic matter (DOM), which can gain relevance in the addition and simulation of other nutrient cycles and processes as leaching. However, here, as a first version of the model, we present a simplification using SOM as a uniform pool. The growth (G, g C m⁻³ day⁻¹) of a biomass pool (B, g C m⁻³) is described according to Monod kinetic,

(9) $$G=\sum _{\mathrm{n}=1}^{\mathrm{N}}\left(g_{\mathrm{m}\mathrm{a}\mathrm{x}}{\left({\displaystyle \frac{\mathrm{S}{\mathrm{f}}_{\mathrm{a}}}{K_{\mathrm{s}}+S}}\right)}_{\mathrm{n}}\right)B$$where g_max (g C g C⁻¹ day⁻¹) is the maximal rate of growth, to which several modifiers are applied (see descriptions below). Substrate (S, g C m⁻³) is the consumable pool, litter, SOM or biomass of soil organism (n), that consumer pool (B) can use but corrected by its available fraction (f_a). All fluxes (N) of consumed C from each S are summed. K_s (g C m⁻³) is related to substrate quality, it gives the content required to get half the maximal growth. This is not related to the amount that will be consumed, because consumed C equals growth + faeces, but shows how dense the material needs to be ‘found’ by the consumer. Available fraction of a S to a consumer (as f_a) is calculated using the fraction from total porosity volume that is accessible for the consumer, by size, minus its fraction that is completely flooded or dry (see Eqs. (14)–16). This availability introduces the concept of physical recalcitrance, highlighting the role that soil structure plays affecting C fluxes in the soil, because SOM decomposition rates modelling use to rely on its chemical recalcitrance, from now on referred just as “recalcitrance”. But physical recalcitrance has proven to be also relevant for the calculation of SOM decomposition rates (Von Lützow et al., 2008), and soil matrix also affect other biotic interactions through the food web by this availability concept.

Rate of increase of a population of meso- or macrofauna depends on generation time (r, K strategies), age distribution of the population, different life stages. Models exist for only some soil fauna species (Osler & Sommerkorn, 2007; Chertov et al., 2017a). To offer a solution that can work for both the microbial biomass and the meso- and macrofauna, we use g_max as equal to the maximal rate of increase in biomass of any population, dB/dt = g_max when resources are non-limiting and assuming the population structure is stable and optimal, equal to what is often stated as the intrinsic growth rate of a species (Birch, 1948).

The net rate of change of a biomass pool is the sum of growth (G), respiration (R) and turnover (death, D_t), and possibly predation (P_d), all in g C m⁻³ day⁻¹: (10) $${\displaystyle \frac{\mathrm{d}\mathrm{B}}{\mathrm{d}\mathrm{t}}=G-R-D_{\mathrm{t}}-P_{\mathrm{d}}}$$R is a function of temperature, through respiration rate (r, g C g C⁻¹ day⁻¹), and biomass, assuming the same temperature sensitivity as growth; this is somewhat different to how it is seen in many models where a food source is turned over with a specific efficiency. From a more faunal point of view, this makes sense: a food source is “consumed”; the consumed material is partly excreted and partly assimilated and spent on respiration and growth (i.e., biomass formation).

(11) $$R=\mathrm{r}\mathrm{B}$$While the death rate d (g C g C⁻¹ day⁻¹) is constant.

(12) $$D_{\mathrm{t}}=\mathrm{d}\mathrm{B}$$Predation depends on biomass of predator or microbivore and is calculated from the growth of the predator (G_pred, from Eq. (9)) plus the fraction of the prey allocated to faeces (f_faec).

(13) $$P_{{\mathrm{d}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{y}}}=G_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}\left(1+f_{\mathrm{f}\mathrm{a}\mathrm{e}\mathrm{c}}\right)$$

Effect of H₂O

Drought or saturation of a pore leads to reduced availability of the C in the pore for its food web consumers. First, the overall effect of hydration is calculated as a modifier (m_H2Otot) in function of volumetric soil moisture (V) and pore volume (P_vol) (after Freytag & Luttich, 1985).

(14) $$m_{\mathrm{H}2\mathrm{O}\mathrm{t}\mathrm{o}\mathrm{t}}=\{\begin{array}{c}4{\displaystyle \frac{V}{P_{\mathrm{v}\mathrm{o}\mathrm{l}}}\left(1-{\displaystyle \frac{V}{P_{\mathrm{v}\mathrm{o}\mathrm{l}}}}\right)\mathrm{f}\mathrm{o}\mathrm{r}{\displaystyle \frac{V}{P_{\mathrm{v}\mathrm{o}\mathrm{l}}}<0.5}}\\ 1for{\displaystyle \frac{V}{P_{\mathrm{v}\mathrm{o}\mathrm{l}}}>0.5}\end{array}$$The activity is always in the pores that are not water-logged, therefore the pore size class that is partially filled with water, and the pore size above that is assumed not yet completely dry (after Kuka, Franko & Rühlmann, 2007).

(15) $$m_{{\mathrm{H}}_2\mathrm{O}}={\displaystyle \frac{P_{\mathrm{v}\mathrm{o}\mathrm{l}\mathrm{A}}}{P_{\mathrm{v}\mathrm{o}\mathrm{l}\mathrm{A}}+P_{\mathrm{v}\mathrm{o}\mathrm{l}\mathrm{W}}}{m}_{{\mathrm{H}}_2{\mathrm{O}}_{\mathrm{t}\mathrm{o}\mathrm{t}}}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{y}\mathrm{f}\mathrm{i}\mathrm{l}\mathrm{l}\mathrm{e}\mathrm{d},}$$

(16) $$m_{{\mathrm{H}}_2\mathrm{O}}={\displaystyle \frac{P_{\mathrm{v}\mathrm{o}\mathrm{l}\mathrm{W}}}{P_{\mathrm{v}\mathrm{o}\mathrm{l}\mathrm{A}}+P_{\mathrm{v}\mathrm{o}\mathrm{l}\mathrm{W}}}{m}_{{\mathrm{H}}_2{\mathrm{O}}_{\mathrm{t}\mathrm{o}\mathrm{t}}}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{o}\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{a}\mathrm{b}\mathrm{o}\mathrm{v}\mathrm{e},}$$where P_volW is the water filled pore volume and P_volA is the aerated pore volume. The availability (f_a) of a substrate to a consumer is defined by the inherent availability of the pore size to the consumer, multiplied with m_H2O. For surface litter these calculations are not possible since the surface litter is not in the soil matrix. However, on days without precipitation, litter humidity is assumed to be related to the soil humidity below, therefore the m_H2O calculated for the microbial biomass is used.

Simulating the variability in g_max

The maximum growth of biota is influenced by different environmental factors. Each one can lead to a modifier (m ∈ [0, 1]) on g_max. It is easy to change, add or turn off specific modifiers according to the soil studied. Here we present a modelling framework focused on abiotic controls of growth rates, but there is room for new add-ons as for example death rate modifiers as a density-dependent microbial turnover. While interaction processes affected by the demographic density of microbial communities (e.g., competition, space constraints) can play also a significant role controlling growth and decomposition rates and improve its modelling (Georgiou et al., 2017), our aim in this work is to link the key roles of fauna and soil structure in C cycle modelling, and together with the hydrology can simulate constraints in biotic interactions, which are also relevant controls in microbial growth and activity.

Simulating the effect of temperature (T)

To simulate the effect of T on growth rate through a temperature modifier (m_T), we use a Q10-shaped curve between maximum tolerable temperature (T_max) and minimum temperature for consumers activity (T_min), set as a default at 0 °C (Franko, 1989), but unlike many models, we assume a plateau above the optimal temperature (T_opt).

(17) $$m_{\mathrm{T}}=\{\begin{array}{cc}0,\hfill & T<T_{\mathrm{m}\mathrm{i}\mathrm{n}}\mathrm{o}\mathrm{r}T\ge T_{\mathrm{m}\mathrm{a}\mathrm{x}}\hfill \\ Q^{\left(T-T_{\mathrm{o}\mathrm{p}\mathrm{t}}\right)/10},\hfill & T_{\mathrm{m}\mathrm{i}\mathrm{n}}\le T<T_{\mathrm{o}\mathrm{p}\mathrm{t}}\hfill \\ 1,\hfill & T_{\mathrm{o}\mathrm{p}\mathrm{t}}\le T<T_{\mathrm{m}\mathrm{a}\mathrm{x}}\hfill \\ \end{array}$$However, temperature also increases respiration (R). To simulate this temperature effect, we assume the same Q10 function but without the plateau; in this way, when T is above the optimum, R increases while growth does not. At some point these lines will cross and cause a net reduction in biomass.

Effect of pH on growth

A good example of an optional effect is the effect of pH: for a system close to a threshold, simulating pH can be very important, assuming a good knowledge of the system. But for well-buffered systems, it is an unnecessary increase in complexity. g_max decreases at low pH for bacteria but increases for fungi (Rousk, Brookes & Bååth, 2009, 2010; Rousk et al., 2010). For this example, we put the thresholds at 8 for fungi and 3 for bacteria, with precision of one decimal, inducing a 10 fold reduction in g_max for a change in pH of 1.

(18) $$m_{\mathrm{p}\mathrm{H}}=1/\left(\left(\mathrm{p}\mathrm{H}–8\right)10\right)\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{g}\mathrm{i}\mathrm{i}\mathrm{f}\mathrm{p}\mathrm{H}\ge 8.1$$

(19) $$m_{\mathrm{p}\mathrm{H}}=1/\left(\left(3-\mathrm{p}\mathrm{H}\right)10\right)\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{i}\mathrm{f}\mathrm{p}\mathrm{H}\le 2.9$$In any other case for bacteria or fungi, m_pH = 1, and if m_pH goes above 1, it is replaced by 1. For engineer saprotrophs, their optimal g_max changes (becoming g_maxEng) with pH (Lavelle, Chauvel & Fragoso, 1995; Chertov et al., 2017b) according to the following equation: (20) $$g_{\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{E}\mathrm{n}\mathrm{g}}=\{\begin{array}{c}0,\mathrm{i}\mathrm{f}\mathrm{p}\mathrm{H}<3\\ \left({\displaystyle \frac{g_{\mathrm{m}\mathrm{a}\mathrm{x}}}{2}}\right)\left(\mathrm{p}\mathrm{H}-3\right),\mathrm{i}\mathrm{f}3\le \mathrm{p}\mathrm{H}<5\\ g_{\mathrm{m}\mathrm{a}\mathrm{x}},\mathrm{i}\mathrm{f}\mathrm{p}\mathrm{H}\ge 5\\ \end{array}$$

Effect of recalcitrance and C:N on g_max

Overall consumption of an organism that can consume different pools is computed by simply adding them up. However, litter is not necessarily as “palatable”, depending on its C:N ratio if not enough N, and on recalcitrance if low in energy, then it is needed to consume more litter, which is calculated through modifiers m_C:N and m_rec, respectively. This is simulated by changing g_max. The equation for m_rec is not necessary and only important if enough data on litter quality is available or the users are interested into looking into the effects of changes in litter quality. The litter pool can be consumed by both bacteria and fungi, and of course also detritivores. Depending on the C:N ratio, the competition between these two is different; this is simulated by the g_max of the bacteria being more variable with C:N ratio. The sensitivity is described by the parameter p_mC:N, between 0 and 1.

(21) $$\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{g}\mathrm{i}:m_{\mathrm{C}:{\mathrm{N}}_{\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{g}}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(1,{\left({\displaystyle \frac{\mathrm{C}:{\mathrm{N}}_{\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{g}}}{\mathrm{C}:{\mathrm{N}}_{\mathrm{l}\mathrm{i}\mathrm{t}}}}\right)}^{{{\mathrm{p}}_{\mathrm{m}\mathrm{C}:\mathrm{N}}}_{\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{g}}}\right)$$

(22) $$\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{a}:m_{\mathrm{C}:{\mathrm{N}}_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{t}}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(1,{\left({\displaystyle \frac{\mathrm{C}:{\mathrm{N}}_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{t}}}{\mathrm{C}:{\mathrm{N}}_{\mathrm{l}\mathrm{i}\mathrm{t}}}}\right)}^{{{\mathrm{p}}_{\mathrm{m}\mathrm{C}:\mathrm{N}}}_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{t}}}\right)$$For litter recalcitrance (Rec_lit), a linear equation instead of a power is chosen so that decay of the recalcitrant litter is 0 if p_mRec = 1 and is unaffected if p_mRec = 0. The reason for choosing a different equation than above is that constrain of labile litter decomposition by C:N ratio should not completely stop decomposition but adjust the decomposition rate, while recalcitrant fraction of litter could remain almost constant for long periods. The following equations determine if the recalcitrant fraction of litter remains stable or if it is affected by decomposers partially or even totally.

(23) $$\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{g}\mathrm{i}:m_{\mathrm{r}\mathrm{e}{\mathrm{c}}_{\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{g}}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(1,1-p_{\mathrm{m}\mathrm{R}\mathrm{e}{\mathrm{c}}_{\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{g}}}\mathrm{R}\mathrm{e}{\mathrm{c}}_{\mathrm{l}\mathrm{i}\mathrm{t}}\right)$$

(24) $$\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{a}:m_{\mathrm{r}\mathrm{e}{\mathrm{c}}_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{t}}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(1,1-p_{\mathrm{m}\mathrm{R}\mathrm{e}{\mathrm{c}}_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{t}}}\mathrm{R}\mathrm{e}{\mathrm{c}}_{\mathrm{l}\mathrm{i}\mathrm{t}}\right)$$

Adding up all these modifying effects on g_max

We assume a complete additivity of the effects, so the different modifiers on g_max are multiplied to get the overall effect, m_tot in Eq. (25). Another optional approach could be to use only the most limiting effect, setting m_tot equal to the lowest modifier and ignoring the rest.

(25) $$m_{\mathrm{t}\mathrm{o}\mathrm{t}}=m_{\mathrm{T}}{m}_{\mathrm{C}:\mathrm{N}}{m}_{\mathrm{r}\mathrm{e}\mathrm{c}}{m}_{\mathrm{p}\mathrm{H}}{m}_{{\mathrm{H}}_2\mathrm{O}}$$

Not assimilated C

The reduction in a substrate equals the growth of the consumer plus the C that goes to faeces (excrements) and to respiration. Fraction to excrement (f_faec) is a parameter of the consumer and assumed constant. However, one consumes more and a larger fraction becomes faeces at a lower substrate quality, for the meso-and macro fauna, because microbes do not produce excrements; the sensitivity of f_faec to C:N ratio is expressed by the modifier m_faec. This is however only relevant for the detritivores and engineers (Eq. (26)) that eat SOM and litter which can contain extremely variable amounts of nutrients; for the predators and herbivores we assume the variability is minimal. For the microbes, it was calculated as an effect on g_max.

(26) $$f_{\mathrm{f}\mathrm{a}\mathrm{e}\mathrm{c}\mathrm{E}\mathrm{f}\mathrm{f}}=f_{\mathrm{f}\mathrm{a}\mathrm{e}\mathrm{c}}+m_{\mathrm{f}\mathrm{a}\mathrm{e}\mathrm{c}}{\displaystyle \frac{\mathrm{C}:{\mathrm{N}}_{\mathrm{S}\mathrm{O}\mathrm{M}}-\mathrm{C}:{\mathrm{N}}_{\mathrm{e}\mathrm{n}\mathrm{g}}}{\mathrm{C}:{\mathrm{N}}_{\mathrm{S}\mathrm{O}\mathrm{M}}}{f}_{\mathrm{f}\mathrm{a}\mathrm{e}\mathrm{c}}}$$

Calculations regarding engineers

Soil changes made by engineer species depend on their body width, but in the model this is simplified using initial parameters for engineers’ effects that must be chosen based on an average width (see Table S2.5); the model then simulates their daily effects using their biomasses. Bioturbation is a function of engineer biomass (B_eng, g C m⁻³), which calculates organic matter moving to deeper layers: litter moving (g C_lit g C_eng⁻¹ day⁻¹) from litter layer to end of burrow, and SOM moving by mixing of soil between layers (g C_SOM g C_eng⁻¹ day⁻¹). In this first version of the model, with only one soil layer, bioturbation works as a C output flow, but in future versions with more layers it could be upgraded to C flows between them.

Burrow volume (PV_B, l m⁻³) is a function of engineer biomass (g C_eng m⁻³) and the ratio of pore volume to engineer biomass (VE_ratio, l g C_eng⁻¹) but towards a maximum (PV_Bmax, l m⁻³). On the other hand, burrow turnover (tPV_B) happens at a constant rate, with average burrow lifespan of 10 years; porosity decreases and burrows become mesopores.

(27) $$\mathrm{P}{\mathrm{V}}_{\mathrm{B}}=\mathrm{m}\mathrm{a}\mathrm{x}\left(\mathrm{m}\mathrm{i}\mathrm{n}\left(\mathrm{P}{\mathrm{V}}_{\mathrm{B}\mathrm{m}\mathrm{a}\mathrm{x}},\mathrm{V}{\mathrm{E}}_{\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}}{B}_{\mathrm{e}\mathrm{n}\mathrm{g}}\right),\left(\mathrm{P}{\mathrm{V}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}-\mathrm{P}{\mathrm{V}}_{\mathrm{t}\mathrm{e}\mathrm{x}\mathrm{t}\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}\right)\left(1-\mathrm{t}\mathrm{P}{\mathrm{V}}_{\mathrm{B}}\right)\right)$$where PV_macro is the current pore volume of the macropores and PV_textmacro is the textural porosity of the macropores (see in next section). This could be improved in future versions including perturbations as the possible effect of heavy rain.

Porosity calculations

The pore volume is distributed in five classes by pore size. Initial pore size distribution is given or measured as the total pore volume (PV, l m⁻³) in each class. The link between aggregation and porosity is hard to quantify. Regelink et al. (2015a) showed for different soils that overall soil porosity is the sum of the textural porosity determined by the proportion of clay, sand and silt fractions and aggregation porosity. They conclude that micropores, which they define <9 µm, are mainly situated within the aggregates, while mesopores are situated between dry-sieved aggregates. While Regelink et al. (2015a) have shown that total micro and mesoporosity (<1,000 µm) increases with total aggregate content, Grosbellet et al. (2011) have provided evidence that pores in the range 30–300 µm decrease with aggregation. Despite of the generally lower ranges for mesopores (9–1,000 µm) described for soil physics (Lal & Shukla, 2004; Regelink et al., 2015a), here mesopores are assumed to be physically accessible to mesofauna body size (ca. 100–2,000 µm), so we consider that mesopores ranging 30–1,500 µm are a reasonable compromise. Based on that, we decided to hypothesize that aggregation increases bacterial and micro- porosity while decreasing mesoporosity. However, we want to emphasize that further experimental studies are needed to establish robust relationships between aggregation and pore size distribution.

Aggregates are not calculated as a pool, but the effect of aggregation is included in the calculation of porosity as described below. The following three porosities contribute to total porosity:

• Textural porosity (PV_text): measured or calculated from % clay and sand.

• Additional aggregation porosity (PV_Ag): all porosity in surplus of textural, can be estimated, for example from PTF (PedoTransfer Function) or calculated empirically from SOM and fungal biomass, that is, mycorrhiza and other fungi, max 2% porosity extra (Eq. (29)). Aggregation (Ag) is the fraction (0–1) of the SOM aggregated calculated as (based on the data from Malamoud et al., 2009):

(28) $$\mathrm{A}\mathrm{g}=\mathrm{m}\mathrm{i}\mathrm{n}\left(1,{\displaystyle \frac{c\left(B_{\mathrm{f}\mathrm{u}\mathrm{n}\mathrm{g}}+B_{\mathrm{m}\mathrm{y}\mathrm{c}}\right)}{B_{\mathrm{S}\mathrm{O}\mathrm{M}}}}\right)$$with an empirical parameter c = 10. The aggregation porosity is then calculated as: (29) $$\mathrm{P}{\mathrm{V}}_{\mathrm{A}\mathrm{g}}=k\mathrm{A}\mathrm{g}{B}_{\mathrm{S}\mathrm{O}\mathrm{M}}$$with k = coefficient (2 l g C⁻¹ m⁻³) based on empirical data (Regelink et al., 2015a, 2015b).

• Bioporosity (PV_B): biopores created by engineers. Pore formation by engineers increases macroporosity, increasing soil layer density, but at the same time reduces mesoporosity as engineers push soil aside and produce casts that are denser than average soil. The relative importance of these two effects depends on the engineers’ activity patterns, and is reflected by the parameter f_PV ∈ (0, 1), which gives the fraction of the change in biopore volume that increases macroporosity. Therefore, the counterpart of the biopore volume (1 − f_PV) PV_B is ‘compensated’ by a decrease in mesoporosity.

Conceptually, the total soil porosity is then the sum of: (30) $$\mathrm{P}{\mathrm{V}}_{\mathrm{t}\mathrm{o}\mathrm{t}}=\mathrm{P}{\mathrm{V}}_{\mathrm{t}\mathrm{e}\mathrm{x}\mathrm{t}}+\mathrm{P}{\mathrm{V}}_{\mathrm{A}\mathrm{g}}+f_{\mathrm{P}\mathrm{V}}\mathrm{P}{\mathrm{V}}_{\mathrm{B}}$$

In the model, pore volume is calculated for each pore size separately.

The volume of micropores (PV_micro) and bacterial pores (PV_bact) increases with increasing aggregation. Apart from creating additional porosity depending on the total amount of aggregated SOM (Eq. (29)), aggregation also increases the relative micro- and bacterial pore volume at the expense of (textural) mesoporosity (PV_meso), therefore not increasing total porosity. This effect is controlled by available pore space between mineral particles (i.e., textural mesoporosity) and we assume that half of this mesoporosity can be affected by aggregation. In both cases, we assume that the increase in porosity due to aggregation is divided equally among micropores and bacterial pores. The pore volume in different size classes is calculated as: (31) $$\mathrm{P}{\mathrm{V}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}=\mathrm{P}{\mathrm{V}}_{\mathrm{t}\mathrm{e}\mathrm{x}\mathrm{t}\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}+\mathrm{P}{\mathrm{V}}_{\mathrm{B}}$$

(32) $$\mathrm{P}{\mathrm{V}}_{\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{o}}=\mathrm{P}{\mathrm{V}}_{\mathrm{t}\mathrm{e}\mathrm{x}\mathrm{t}\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{o}}-\left(1-f_{\mathrm{P}\mathrm{V}}\right)\mathrm{P}{\mathrm{V}}_{\mathrm{B}}-{\displaystyle \frac{\mathrm{A}\mathrm{g}}{2}\mathrm{P}{\mathrm{V}}_{\mathrm{t}\mathrm{e}\mathrm{x}\mathrm{t}\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{o}}}$$

(33) $$\mathrm{P}{\mathrm{V}}_{\mathrm{m}\mathrm{i}\mathrm{c}\mathrm{r}\mathrm{o}}=\mathrm{P}{\mathrm{V}}_{\mathrm{t}\mathrm{e}\mathrm{x}\mathrm{t}\mathrm{m}\mathrm{i}\mathrm{c}\mathrm{r}\mathrm{o}}+k{\displaystyle \frac{\mathrm{A}\mathrm{g}}{2}{B}_{\mathrm{S}\mathrm{O}\mathrm{M}}+{\displaystyle \frac{\mathrm{A}\mathrm{g}}{4}\mathrm{P}{\mathrm{V}}_{\mathrm{t}\mathrm{e}\mathrm{x}\mathrm{t}\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{o}}}}$$

(34) $$\mathrm{P}{\mathrm{V}}_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{t}}=\mathrm{P}{\mathrm{V}}_{\mathrm{t}\mathrm{e}\mathrm{x}\mathrm{t}\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{t}}+k{\displaystyle \frac{\mathrm{A}\mathrm{g}}{2}{B}_{\mathrm{S}\mathrm{O}\mathrm{M}}+{\displaystyle \frac{\mathrm{A}\mathrm{g}}{4}\mathrm{P}{\mathrm{V}}_{\mathrm{t}\mathrm{e}\mathrm{x}\mathrm{t}\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{o}}}}$$

Volume of inaccessible pores is assumed to be constant and equal to PV_textinac.

These changes are calculated daily to give a dynamic feedback to the hydrology and to the distribution of each C source among pore classes, affecting its availability.

Leaching

Water leaving one soil layer (n) is moved to the layer below (n + 1). Dissolved organic and inorganic compounds are a complex matter to simulate since they are strongly dependent on the pH and the mother-material, that is, clay and Ca rich or not. Nonetheless, in many systems simulating leaching of N and DOM is highly relevant. Unless better data are available, we suggest the following, semi-empirical method:

Dissolved organic matter can be simulated in relation to the CO₂ released as total respiration (R_tot, g C m⁻³ day⁻¹) based on the consideration that high “activity” in the soil is related to high R_tot. This is calculated as a fraction (f_DOM) of R_tot entering the DOM pool, similar to the concepts used in the LIDEL model, in addition to the directly exuded DOM (C_Exud, g C m⁻³ day⁻¹) which is an input (from data or a vegetation model). Assuming a short half-life of DOM and semi-empirically, because daily concentration is not “equal” to daily production (DOM_p, g C m⁻³ day⁻¹) but is linearly related to the daily production, we consider: (35) $$\mathrm{D}\mathrm{O}{\mathrm{M}}_{\mathrm{p}}={\mathrm{C}}_{\mathrm{E}\mathrm{x}\mathrm{u}\mathrm{d}}+f_{\mathrm{D}\mathrm{O}\mathrm{M}}{R}_{\mathrm{t}\mathrm{o}\mathrm{t}}$$

DOM has a short half-life but the dissolution is even faster (hours). We assume the daily concentration is in equilibrium between dissolved and adsorbed (DOM_ad) depending on adsorption coefficient K_D of the soil (m³ kg⁻¹ soil). Similar to the modelling in Orchidee-SOM (Camino-Serrano et al., 2018) we assume: (36) $$\mathrm{D}\mathrm{O}{\mathrm{M}}_{\mathrm{a}\mathrm{d}}=K_{\mathrm{D}}\mathrm{D}\mathrm{O}\mathrm{M}$$with K_D depending on the minerals and pH. More clay (fClay fraction) means less mobile DOM, and lower pH is also a cause of less mobile DOM.

(37) $$K_{\mathrm{D}}=a_{\mathrm{K}\mathrm{D}}-b_{\mathrm{K}\mathrm{D}}\mathrm{p}\mathrm{H}+c_{\mathrm{K}\mathrm{D}}f\mathrm{C}\mathrm{l}\mathrm{a}\mathrm{y}$$with values 0.001226, 0.000212 and 0.00374 respectively for a_KD, b_KD and c_KD, from Camino-Serrano et al. (2018).

DOM leaching is calculated as the volume of water leaching to a lower layer multiplied with the concentration of dissolved DOM.

Calculation order

The sequence of function sets used by the model to calculate all carbon fluxes and ecosystem changes is as follows:

1. Calculate the pore size fractions in 5 classes and the associated pore surface areas

2. Calculate the water volume of the relevant pore size

3. Use the precipitation leaching to calculate DOM leaching

4. Calculate for each biota group the accessibility of each of the pools it consumes

5. Calculate the g_max depending on temperature, H₂O, C:N, pH and recalcitrance

6. Solve the 12 differential equations for increase/decrease of all C pools

7. Update all C pools

8. Calculate the new C:N and recalcitrance of each pool

9. Calculate engineering effect

   • i1) Update macropores

   • i2) Update SOM from bioturbation

10. Calculate other changes in pore size distribution from weather or management

The KEYLINK core model consists of steps d to i; steps a, b, c and j are add-ons that could be replaced by other models (e.g., water flow model) coupled to KEYLINK. Steps a–c are used to calculate the distribution of porosity between the pore classes, the hydrology and daily soil water content (distributed among pore classes), and then step d calculates how that is affecting the availability of each C source to its consumers. That couples soil structure and hydrology with trophic interactions, allowing the resolution of differential equations for C fluxes.

Model coding and output

KEYLINK consists of a relatively limited, freely downloadable Python code (available at: https://github.com/Plant-Root-Soil-Interactions-Modelling/KEYLINK). Each of the modifiers on growth, that is, temperature, pH, H₂O, recalcitrance and C:N, as well as the primal shape of the growth equations can be adapted towards specific questions or ecosystems. The inputs in the current version are read from data-files but are easy to link to a mechanistic model. The output of the current version consists of all daily C pools as well as the main C fluxes. KEYLINK is also available as a stand-alone executable model, allowing it to be called from models in other languages. A single run of ten years could take less than one minute (depending on computing power). In this version, the average results over one hundred runs are calculated but also all daily outputs of each run are saved.

Model parameterization

The first version of KEYLINK model has been parameterized for a Scots pine forest stand situated in Brasschaat, in the Campine region in Belgium (51°18′ N and 4°31′ E) but without modelling the complete forest ecosystem (to simplify the interpretation of the results from KEYLINK). The goal of this parameterization was a model verification, not a model application for which a complete integration with an aboveground model or detailed above ground data would be necessary.

The soil of the Brasschaat forest is sandy but with high ground water table, so trees are generally not water-limited, but the topsoil is often dry. The soil is acidic (pH 3.5). The trees were planted around 1930 and formed a rather sparse vegetation in 1999, with leaf area index (LAI) ranging from 2.1 to 2.4.

For this model run, we used the following input data from the stand (Table 1). In this case, we did not use measured or modelled growing trees but constant input of aboveground and belowground litter (measured value). The top 90 cm of soil from the Brasschaat forest was analyzed in Janssens et al. (1999). Earthworm biomass, used in this case as an example of ecosystem engineers, is extremely low due to the low pH, it was not measured since 1993 by Muys (1993), but these data are used since there is no reason to expect there was a marked change.

Data availability on soil pools, biology and functioning is generally low, and it is currently not possible to find a dataset describing in detail, and with small error margins, the temporal evolution of all different soil biological compartments and SOM pools. Available data are often incomplete, or based on rough estimates, for example, from semiquantitative DNA analysis for microbial abundance in soils. To deal with this issue, a quite pragmatic approach combining different estimates from different sources is appropriate for most datasets where the soil is not the key focus, but a means to improve the simulation of an ecosystem.

The daily loss of water by evapotranspiration was calculated using an equation for potential evapotranspiration based on Thornthwait (1948) in this study.

Model calibration

Once the model is parameterized for an ecosystem, the next step is to optimize that model, calibrating the fit of its simulations to the ecosystem data. The optimization included in the KEYLINK model follows a Bayesian procedure, using the version of Markov Chain Monte Carlo (MCMC) method known as the Metropolis-Hastings random walk with reflection algorithm (Christian & Casella, 1999; Van Oijen, Rougier & Smith, 2005; Van Oijen, 2008).

A pragmatic assumption is that the starting values of the C pools (including the soil fauna initial biomasses) are at steady state for a given date (most often spring or summer, it would be unrealistic to keep the values constant through the year as they fluctuate with climate). The simplest calibration of any ecosystem can be done by assuming these 11 carbon pools (litter, SOM and the nine functional groups in food web) need to be stable over the simulated years, for example, for 9 years that gives us 99 data points by taking the same value for each C pool every year (Table 2). Initial litter, SOM and biomasses of bacteria, fungi and engineers were taken from the references cited in Table 1. For other C pools, data were estimated using measured data for previous C pools and similar proportions between C pools as in the Swedish pine forest in Persson et al. (1980); predator biomass was assumed to be the 20% of all biomass in their consumed C pools. Errors were assumed as a percentage of biomass, 10% for predators, 12.5% for litter and SOM and 20% for the rest C pools.

It is common to apply a correction (“burn-in”) deleting part of the posterior, for example, the first half of the runs, to avoid the effect of the starting distribution (Gelman & Shirley, 2011). For this calibration, a sample of the last one hundred accepted parameter vectors was taken from the posterior distribution, and it was used for all further model runs, so every run was performed with 100 different parameter sets.

Input parameters of species

The KEYLINK model framework is conceptualized as an adaptable framework. Each user needs to determine for their specific site and questions the main drivers and pools required. Depending on the dataset, it is in general better to use less pools and equations if sparse data are available. Moreover, KEYLINK is not a soil fauna model and was not designed to simulate specific soil fauna species in detail. The soil fauna groups used consist of a wide range of species, for which average data are used. For a description of the species categories, we refer to the review on the KEYLINK concepts (Deckmyn et al., 2020).

Microbes and meso-macro fauna have a temperature curve using an optimum, minimum and maximum temperature. Each soil biota group also has its own maximum growth rate, C:N ratio, respiration rate and size. Death rate (d) is the inverse of turnover, mostly given in days. In Supplemental File S1 we briefly review main input parameters. We propose setting K_s, the concentration of the food source at which growth rate is half the maximum, equal to the existing concentrations for all meso- and macro fauna, so assuming growth could double at unlimited food source. But for microbial biomass the difference between growth of bacteria on a petri-dish unlimited in nutrients compared to field data of soil microbes clearly indicates that g_max in the soil is not comparable to laboratory data; if such data of g_max are used, the K_s should be increased considerably.

Calibration for Brasschaat pine forest

We show here the results from a calibration towards data measured and assumed, using proportions between fauna groups in Persson et al. (1980), in the Brasschaat Scots pine stand in Belgium. This forest stand is relatively well described in many publications concerning the trees and the total ecosystem fluxes, but less concerning the soil and very little was measured on soil fauna. Therefore, the partially assumed data refers to a hypothetical ecosystem that does not fully fit with reality in Brasschaat forest. We use this forest as an example of how the KEYLINK model can be used to improve our understanding of the system even when detailed soil faunal data are limited.

The parameters g_max, K_s and r are linked (increasing g_max has a similar effect to decreasing K_s or r). However, g_max or r ranges can be found in literature relatively easily. Therefore, we use fixed values for K_s (see Supplemental File S2) and parameterize g_max within the known limits. In this way, the number of parameters to be calibrated is 9, which is a reasonable number for most cases where limited data to calibrate towards are available. Of course, a user could decide to optimize more parameters if more data are available. A useful “rule of the thumb” is limiting the number of parameters to the square root of the number of calibration data available (Jörgensen, 2009), which means we can get a reasonable result for nine parameters assuming 81 data points.

In our case, no measurements of growth rates were available and information in the literature was scant. Therefore, we deliberately defined wide ranges for the prior values of each parameter to cover all the possible values found in the literature (Chuine, 2000; Linkosalo, Lappalainen & Hari, 2008). For species for which no prior parameter information was available, we assumed parameter values equal to the mean value of the range. The initial uncertainty of each parameter is quantified in terms of a prior probability distribution with lower and upper bounds. Because of lack of detailed knowledge, we assumed the distribution as uniform and non-correlated.

The g_max values were optimized using the prior range for g_max (Table 3). The data used to calibrate against were chosen to give a ‘standard’ procedure, so limited to biomass of the different C pools. Including all available data s.a. soil respiration, soil humidity could improve the run for Brasschaat, but would not be a representative run for the model. Other parameter settings, for example, sensitivity to C:N and recalcitrance, were based on model runs of the Brasschaat site by Deckmyn et al. (2011).

We ran the model for the time period 1999–2008, because this was the period in which the forest was still clearly dominated by Scots pine; since then a transition to more deciduous trees has been taking place. We calibrated towards stable C pools over the ten years for all C pools, with an allowed error margin of 20% for all faunal pools, except 10% for predators, and 12.5% for litter and SOM. Daily climate data (temperature and precipitation) were used. The full range of input data can be found in Supplemental File S2, except climate data, which can be downloaded with the model. Choosing to calibrate towards one or more pools can yield different results, and it depends on the end-user’s goal which calibration is preferred.

Model evaluation

Although coupling KEYLINK to real or simulated data of the aboveground ecosystem would yield more realistic results, in this exercise we used KEYLINK as a stand-alone model with quite constant input (e.g., litter, plant water uptake) to minimize the feedback effects and give a clear view on the model behaviour. This is a model evaluation, not a full model validation.

After calibration to the Brasschaat dataset, a set of scenarios was performed to evaluate the model: (I) Basic results; (II) Sensitivity to initial soil structure; (III) Changing initial litter C:N ratio; (IV) Changing initial litter recalcitrance; (V) Changing soil pH; (VI) Excluding predators.

Scenario I was done with the reference input parameters (Supplemental File S2), and used as a basal one to be compared with the other five alternative scenarios: scenario II with higher clay content in the soil (clay 15%); scenario III with lower litter C:N ratio (40); scenario IV with lower litter recalcitrance (20%); scenario V with higher pH (5.9); and scenario VI without predators by setting its initial biomass to 0 (B_pred = 0).

In each one of the five alternative scenarios, input parameters were the same than in the basal scenario, except for the parameter changed to generate the new scenario (see Supplemental File S2). All the six scenarios were run 100 times using a sample of 100 parameter vectors from the posterior distribution of the calibration, consisting each run in a simulation of 10 years at a daily time-step (3,653 days). Then, averages of biomasses were calculated for each C pool among the 100 simulations of 10 years, for each scenario, comparing the effects of disturbances on average values.

Calibration

The model was run ca. 100,000 times with different parameter settings sampled from the prior parameter distribution. A sample of the posterior distribution was taken with the last 100 accepted parameter vectors for g_max (Table 4).

The optimization showed the link between the different groups of soil biota, for example, a high g_max for bacteria was coupled to a high g_max for bacterivores. The alternative five scenarios compared to the basal one can show very different results concerning specific C pools (Table 5). Running the model 100 times using the sample of the g_max values resulted in predictions with a quite wide range (Table 6).

Basic results

Mycorrhiza, herbivores and detritivores are relatively uncoupled, though influenced by predators, and follow the yearly climate curves. The bacterial and fungal biomasses are very strongly linked. The high g_max of bacteria allows steep peaks, which are generally followed by peaks in bacterivore biomass. As we used constant litter input into the soil and used a calculated constant fraction of potential evapotranspiration as water uptake from the soil, it cannot be expected that these results follow the normal annual trends in fluctuations of those fluxes. This can, at least partially, explain the relatively low bacterial biomass found in our results, since the bacteria would profit most from a rapidly changing environment, but under some unrealistic simulated conditions fungi could be displacing bacteria by competitive exclusion. For more realistic simulations the model can be coupled with other models that give that information as outputs, or with measured datasets.

All C pools tend to reach some stability after the first years, suggesting the model is well-balanced; however, stability values seem to be more sensitive to changes in g_max parameters for some pools (e.g., SOM). The set-up of the model, where we only calibrate the faunal g_max, does not allow calibration towards different ratio of litter and SOM decay. This depends on the uncalibrated parameter fragmentation, the sensitivity to recalcitrance, but also the temperature used for the litter and SOM. Here we used the same temperature, while in reality it would be expected that there would be certain differences in mean temperatures or their variability at different depths.

There was a high variability within the 100 simulations of each scenario, for example, basal scenario (Fig. 3), calling into question the reliability of predictions. We suggest here some theoretical predictions based on this example, but it is clear that for a realistic application it will be necessary to improve the calibration, by using more detailed data or by linking the model to a vegetation model. This highlights the relevance of developing databases including enough details for the key parameters of the different parts of the soil system (i.e., soil structure, hydrology, food webs).

An overview of all C pools under the different simulation scenarios (Fig. 4) shows how changing one input parameter at a time influences the results. It must be taken into account that KEYLINK was run as a stand-alone model, which can explain why some of the resulting outputs seem not very realistic; linking it to a model or more detailed data of the aboveground ecosystem would greatly influence the results, but would not allow clear interpretation of the model functioning due to feedbacks. Since our goal here was to introduce the belowground model itself rather than a realistic application to a particular case, we chose to avoid those feedbacks with other parts of the ecosystem that are required for more realistic simulations. To further elucidate these effects and to show some of the potential outputs the model can give, we show a few of the most interesting fluxes (Table 7).

The simulated scenarios showed that increasing clay content (i.e., changing initial soil structure) resulted in an increase in water content (Fig. 5) and a decrease in litter and SOM decay while fungal/bacterial ratio decreased. In fact, this scenario caused the largest change in soil water content, showing the sensitivity of the system to initial soil structure and its crucial role for soil hydrology. On the other hand, lowering litter recalcitrance or C:N ratio resulted in an increase in microbial biomass, mainly fungi, which caused an increase in litter decay, while SOM decay did not show a clear change. The scenario with higher pH allowed engineer species (and predators with them) to increase in biomass, altering the soil structure with an increase in macroporosity, which caused a clear decrease in SOM stabilization. Finally, the exclusion of predators totally changed the soil porosity and the trophic interactions along the food web, causing the largest increase in the decay of SOM and litter, which highlights the crucial role of predators in the regulation of the soil C cycle.