Land use changes and the intensification of agricultural production have dramatically altered the flow of nutrients resulting in unprecedented transcontinental inter-watershed and intra-watershed transfers of nitrogen (N), phosphorus (P), and other nutrients with fertilizer, harvest product, and pollutant flow (Galloway et al., 2008). Nutrient cycling in soils plays a major role in the control of these flows. Carbon (C), N, and P cycling are intimately linked through soil, plant and microbial processes. These processes affect the level of inorganic N and P and the C:N and C:P ratios of SOM. To realistically represent these C, N, and P transfers in river-basin scale models such as the Soil Water Assessment Tool (SWAT+, Arnold et al., 1998). a comprehensive integration of the cycling of these nutrients through soil organic matter (SOM) is required.

Most conceptual and quantitative SOM cycling models compartmentalize soil C and N in pools with different, yet stable, turnover rates and C:N ratios (Paul et al., 2006, McGill et al., 1981; Parton et al., 1988; Verberne et al., 1990). Incubation experiments also suggest the existence of pools with varying turnover rates (e.g. Collins et al, 2000); however, Six et al. (2002) indicated after an extensive literature review that the success at matching measurable and modelable SOM pools has been minimal. Furthermore, the division of SOM in pools has been criticized on mathematical grounds, as continuous turnover rates distributions can be artificially represented by discrete pools (Bruun and Luxhoi, 2006). These criticisms do not deny the existence of pools but rather emphasize the difficulty in establishing generalized methods to measure or predict their size and turnover rate. This unpredictability can limit the applicability of multi-pool SOM sub-models if the parameterization for different agricultural soils, pasturelands, forestlands and organic soils is uncertain or requires intensive calibration.

In this chapter, we describe the one-pool SOM sub-model implemented in SWAT+. This sub-model is conceptually based on the model described in Kemanian and Stockle (2010), and was adapted to the SWAT+ algorithms and integrated to the cycling of N and P.

The previous versions of SWAT+ have separate pools for soil organic N (two pools, active and stable) and organic P (one pool) and do not include explicitly $S_C$, $R_C$, or $M_C$ pools. Nitrogen and P decompose separately, excepting the common effect of soil moisture and temperature on the actual decomposition rate and the indirect control of decomposition through $R_{CN}$ and $R_{CP}$. The higher the ratios the lower the decomposition rate, but no immobilization of mineral N or P occurs. Of the N and P decomposed from residues, 20% is allocated to the active soil organic N and the soil organic P pools, respectively, and 80% is allocated to the nitrate ($NO_3$) and P in the soil solution. The 20% allocated to the N and P organic pools is roughly similar to the humification of C described in the new sub-model.

In the new sub-model, mineralization of N from organic pools feeds the $NH_4$ pool instead of the $NO_3$ pool. This will possibly increase volatilization of ammonium because in SWAT the latter is tightly linked to nitrification. Since mineralization and immobilization are explicitly calculated, the fluctuations in mineral $NO_3$, should be more realistically simulated, which has obvious implications for $NO_3$ transport and other processes across the watershed. Tillage accelerates SOM turnover and mixes layers according to the tillage depth. In the previous version tillage had some effects on the surface properties (residue cover) and mixing only involved two layers. No effect of tillage on the organic N or P decomposition rate was previously simulated. The new sub-model also enhances the capacity of SWAT+ to simulate the impact of tillage on nutrient cycling.

The sub-model described here pertains to soil processes at the hydrologic unit response level. The new sub-model has one pool for soil organic C, N, and P ($S_C,S_N,S_P$ respectively, kg m$^{-2}$), and separate pools for residue and manure C, N, and P($R_C,R_N,R_P$ and $M_C,M_N,M_P$, respectively, kg m$^{-2}$). The pools are not separated in active and stable pools. Microbial activity on SOM, manure, and residues decomposes simultaneously the organic C, N, and P. Decomposition of residues may add to (mineralization) or subtract from (immobilization) the inorganic N and P pools, depending on the C:N and C:P ratio of the decomposing pools ($R_{CN},R_{CP},M_{CN}$, and $M_{CP}$, respectively) and those of the SOM or destiny pool ($S_{CN},S_{CP}$), and on the humification rate (h). The humification rate represents the fraction of C in the decomposing residues that are effectively incorporated in the $S_C$ pool, the reminder being respired as $CO_2$. The balance equations for each soil layer for the organic residue, manure, and soil pools are as follows.

Residue and manure pools (represented as R and M) gain mass through additions ($I_R$ and $I_M$ for residue and manure, respectively) and lose mass from decomposition following first order kinetics:

$\frac{dR_c}{dt}=I_{RC}-f_Ek_RR_C,$ 3:5.1.1a

$\frac{dM_C}{dt}=I_{MC}-f_Ek_MM_C,$ 3:5.1.1b

$f_E=(f_Tf_Wf_0)^{fp},$ 3:5.1.1c

where $f_E$ is the combined effect of the soil factors temperature ($f_T$), moisture ($f_W$), and aeration ($f_0$), $fp$ (0.67) is a power regulating the multiplicative effect of the three environmental factors, and $k_R$ (0.05 day$^{-1}$) and $k_M$ (0.025 day$^{-1}$) are the optimum decomposition rate (day$^{-1}$) for residues and manure, respectively. The addition of residues (or manures) from different sources are bulked with the existing pool and not tracked separately, with rates shown above assumed to apply to all sources. If these decomposition rates are made residue- or manure-specific, every time there is an addition of residues or manure the effective decomposition rates would have to be calculated as weighted averages of the rate of the already decomposing pools and that attributed to the newly added materials, with a weighting function that should change over time, or each pool tracked separately, both rather impractical propositions.

For the soil pools, the differential equations are as follows:

$\frac{dS_C}{dt}=h_Rf_Ek_RR_C+h_Mf_Ek_MM_C-k_SS_C,$ 3:5.1.2a

$\frac{dS_N}{dt}=\frac{h_R f_E k_R R_C+h_Mf_Ek_MM_C}{S_{CN}}-k_SS_N,$ 3:5.1.2b

where $h_R$ and $h_M$ are the residue and manure humification rates (kg kg$^{-1}$) and $k_S$ is the apparent organic matter decomposition rate (day$^{-1}$). The humification rates depend on the current $S_C$ and a reference $S_C$ ($S_{CC}$):

$h_R=h_x(1-(\frac{S_C}{S_{CC}})^{\alpha}),$ 3:5.1.3a

$h_x=0.09(2-e^{-5.5clay}),$

$h_M=1.6h_R,$ 3:5.1.3b

$S_{CC}=S_{BD}Z_l(0.021+0.38clay)$ 3:5.1.4

Clay is the soil layer clay fraction (kg clay kg$^{-1}$ dry soil), $S_{BD}$ is the soil layer bulk density (kg m$^{-3}$), and $Z_l$ is the soil layer thickness (m). The constant $\alpha$ (default $\alpha=6$) modulates the response of the humification the current $S_C$. The maximum attainable residue humification is approximately 0.18 kg kg$^{-1}$. The same humification is used for below and aboveground residues. The humification of C from manure was assumed to be 60% higher than that of fresh residues (approximately 0.29 kg kg$^{-1}$). The reference $S_C$ or $S_{CC}$ depends linearly on the soil layer clay fraction (Eq. 3:5.1.4) as proposed by Hassink and Whitmore (1997). When $S_C=S_{CC}$ then the humification is 0 and no accumulation of $S_C$ above $S_{CC}$ can occur (Kemanian and Stockle, 2010). In addition, this approach assumes a dependence of humification on the clay fraction, a dependence represented in many different ways in other models (e.g. Jenkinson, 1990; Bradbury et al., 1993). The control of $h_R$ and $h_M$ in Eq. [3.5.1.3] with $\alpha=6$ implies a non-linear response of $h_R$ and $h_M$ to $S_C$. It is likely that a better formulation is needed for these functional equations to apply in organic horizons, tropical soils, or soils with a high proportion of volcanic ashes.

The C:N ratio of the newly formed organic matter or $S_{CN}$ ranges from 8.5 to 14. The $R_{CN}$ and $M_{CN}$ determine how high $S_{CN}$ can be when there is no mineral N available that can be an N source for the soil microbes. This estimate of $S_{CN}$ is reduced depending on the mineral N availability, as follows:

$S_{CN}=8.5+2.7(1-\frac{1}{1+(\frac{R_{CN}}{110})^3})(1+\frac{1}{1+(\frac{N_{min}}{8})^3})$ 3:5.1.5

where $N_{min}$ is the mineral N in the layer (mg N kg$^{-1}$ soil in this equation). The first term within brackets represents the control of the residues and the second term represents the control of mineral N on $S_{CN}$ so that the higher the ratio of $R_{CN}$ the higher the resulting $S_{CN}$, and the higher the amount of $N_{min}$ the lower $S_{CN}$, with $S_{CN}$ ranging from 8.5 to 14 kg C kg$^{-1}$ N in SOM. For manures, $M_{CN}$ substitutes for $R_{CN}$ and 55 substitutes for the constant 110. The ratio $S_{CP}$ is derived from $S_{CN}$ by assuming that the newly formed SOM will conserve the N:P ratio of the decomposing residue and manure. Thus, the C:N and C:P ratios of SOM are not constant but fluctuates according to these equations during the simulation.

Residues and manure may not supply sufficient N and P to satisfy the $S_{CN}$ (Eq. 3:5.1.2b) and $S_{CP}$ (Eq. 3:5.1.2c), in which case N and P will be mined from the inorganic soluble pools. If the inorganic pools cannot supply N and P for decomposition to proceed, then the decomposition rate of residue and manure is reduced. Therefore, lack of mineral N and P in solution may slow down decomposition. The net mineralization is obtained from:

$MIN_{RN}=\frac{dR_C}{dt}(\frac{1}{R_{CN}}-\frac{h_R}{S_{CN}}),$ 3:5.1.6a

$MIN_{MN}=\frac{dM_C}{dt}(\frac{1}{M_{CN}}-\frac{h_M}{S_{CN}}),$ 3:5.1.6b

where $MIN_{RN}$ and $MIN_{MN}$ are the net mineralization rates (kg m$^{-2}$ day$^{-1}$) from decomposing residues (Eq. 3:5.1.6a) and manure (Eq. 3:5.1.6b), respectively, with negative values indicating immobilization and positive values net mineralization. The same equation applies for P with appropriate C:P ratios. The SOM decomposition rate ($k_S$) is calculated from:

$k_S=k_xf_{tool}f_E(\frac{S_C}{S_{CC}})^{\beta}.$ 3:5.1.7

The apparent $S_C$ turnover rate ($k_S$) is scaled down from an optimum of 4.5% yr$^{-1}$ for undisturbed soils ($k_x$ = 0.000123 day$^{-1}$ so that $k_x$ × 365 = 0.045) based on environmental conditions and $S_C$. The power $\beta$ modulates $k_S$ (default $\beta$= 0.5) so that when $S_C$ is low, the turnover rate slows sharply. Substituting Eq. 3:5.1.7 for $k_S$ in Eq. 3:5.1.2a with $\beta$= 0.5 renders a kinetics of order 3/2 for $S_C$ decomposition, as opposed to the typical first order kinetics commonly use in SOM decomposition models (see Kemanian and Stockle, 2009). Tillage can enhance $k_S$ through the factor $f_{tool}$. This factor is calculated independently for each soil layer and depends on the tillage tool mixing factor ($f_{mix}$, range 0 to 1) and the soil texture. The mixing factor and the tillage depth determine the fraction of the soil layer that is mixed by a tillage operation. The $f_{tool}$ basal value is 1 and it is enhanced immediately after a tillage event based on the estimated cumulative $f_{mix}$ (or $f_{cm}$):

$f_{tool}=1+(3+5e^{-5.5clay})(\frac{f_{cm}}{f_{cm}+e^{1-2f_{cm}}})$ 3:5.1.8

The factor ($f_{tool}$) is reduced on a daily basis based on soil moisture to simulate soil settling. If $f_{tool}$ > 1 and a tillage operation is executed, the corresponding $f_{mix}$ has to be added to the current $f_{cm}$. This requires solving for $f_{cm}$ by inverting Eq 3:5.1.6 before recalculating $f_{tool}$. Since Eq 3:5.1.6 is non-linear on $f_{cm}$ the solution has to be obtained iteratively. In the SWAT+ code the solution for $f_{cm}$ has been approximated by a functional equation to prevent recurrent iterations that can be computationally expensive. The SOM decomposition always causes net mineralization to be positive because humification is assumed to be zero. In actuality, the humification is not zero but since no explicit microbial pool is considered, the decomposition rate is an $apparent$ decomposition rate that represents the net loss of C from the SC and underestimates the true turnover rate of organic C, N, and P (Jenkinson and Parry, 1989).

The balanced equation for $S_C$ in a single layer (Eq. [3:5.1.2a]) can be combined with Eq. [3:5.1.3a] and Eq. [3:5.1.7] and with a few assumptions allow for an analytical solution to the differential equation. The soil clay content and the C input rate were assumed to be constant and the powers $\alpha$ and $\beta$ set to 1, so that the balance equation can be solved explicitly for $S_C$. The time step was considered to be a year and we assumed that residues and manure fully decompose in that time frame so that $k_R=k_M=1$. Since the influence of residues and manure on the $S_C$ balance is similar we assumed that C inputs are only through residues. With these substitutions Eq. [3:5.1.2a] becomes:

$\frac{dS_C}{dt}=h_xR_C-\frac{h_xR_C}{S_x}S_C-\frac{k}{S_x}S^2_C$ 3:5.3.1

The constant k substitutes for $k_xf_Ef_{tool}$. This differential equation can be solved analytically, with integration rendering the following solution:

$S_C(t)=\frac{h_xR_C}{2k}[(\frac{\phi e^{\frac{\gamma_t}{S_x}}-1}{\phi e^{\frac{\gamma _t}{S_x}}+1})\sqrt{1+\frac{4kS_x}{h_xR_C}}-1]$ 3:5.3.2

$\gamma=h_xR_C\sqrt{1+\frac{4kS_x}{h_xR_C}},$

$\phi=-\frac{2kS_C(t=0)+h_xR_C+\gamma}{2kS_C(t=0)+h_xR_C-\gamma}$

The integration constant $\phi$ depends on the initial $S_C$. The steady state solution for Eq. 3:5.3.2 is:

$S_C=\frac{h_xR_C}{2k}\sqrt{1+\frac{4kS_x}{h_xR_C}}-1]$ 3:5.3.3

The ratio $\frac{h_xR_C}{k}$ is the equilibrium $S_C$ that would be obtained if neither $h_C$ nor $k_S$ had a dependence on $S_C$. As $R_C$ increases, the value of the fraction $\frac{4kS_x}{h_xR_C}$ will get smaller. Therefore, the squared root term in Eq. 3:5.3.3 can be approximated as $1+\frac{1}{2}\frac{4kS_x}{h_xR_C}$ by preserving just the first two terms of the binomial expansion, from which Eq. 3:5.3.4 can be re-written as:

$S_C=\frac{h_xR_C}{2k}(1+(\frac{1}{2})\frac{4kS_x}{h_xR_C}-1)=S_x$ 3:5.3.4

Therefore, as $R_C$ increases $S_C$ approaches $S_x$ asymptotically (Figure 1). Taking as a reference a soil layer with $S_x$ = 25 Mg C ha$^{-1}$, $h_x$ = 0.2, and $k=0.01$ $y^{-1}$, it can be seen in Figure 1 that doubling $h_x$ and $k$ have a similar effect but of opposite sign such that the equilibrium $S_C$ increases with increasing $h_x$ and decreases with increasing $k$. In both cases the increase and the decrease in $S_C$ are less than proportional to the increase in these two parameters. The equilibrium $S_C$, however, is very sensitive to changes in $S_x$, which makes this variable critical for a correct representation of $S_C$ dynamics. This formulation is a mathematical representation of the concept of $S_C$ saturation (Hassink and Whitmore, 1997; Six et al., 2002), enhanced with a control of the decomposition rate by $S_C$.

The transient trajectory of $S_C$ is controlled by the quotient of the two exponential terms in Eq. [3:5.3.5].

$d=1-\frac{\sum^N_{i=1}(O_i-S_i)^2}{\sum^N_{i=1}(|S_i-\overline O|+|O_i-\overline O|)^2}$ [3.5.3.5]

For a given $S_x$ and initial $S_C$, increasing the inputs ($R_C$) changes the steady state $S_C$ with decreasing marginal increments as $S_C$ approaches $S_x$, yet the steady state condition is approached faster with higher inputs (Figure 3:5-1). For a given $R_C$, changing $S_x$ has a substantial impact on the rate of change of $S_C$ when the inputs are medium to high (Figure 3:5-1) but a minor effect if inputs are too low. This formulation strongly suggests that soils with higher carbon storage capacity (higher $S_x$) that are currently depleted of $S_C$ should be the primary targets for storing $S_C$, or that soils with low $S_x$ may store carbon quickly for a few years but the rate of gains will decrease earlier than in soils with higher $S_x$.

The conditions for which the $S_C$ can be modeled analytically as shown here are very restrictive. The numerical solution implemented in the model is more flexible as the constants $\alpha$ and $\beta$ are allowed to differ from 1. The model can be expanded to accommodate saturation of different SOM pools, instead of just one uniform pool, as strongly suggested by the results and analysis of Stewart et al. (2008). Yet, this will require a level of parameterization for which we consider there is simply not sufficient information for a realistic implementation in numerical models.

Figure 3:5-1. Equilibrium soil organic carbon (SC, Mg C ha$^{-1}$) for the steady state condition (Eq. [3:5.3.3]) with different values for humification ($h_x$, kg kg$^{-1}$), SOM apparent turnover rate ($k$, yr$^{-1}$), and saturation soil organic carbon ($S_{CC}$, Mg C ha$^{-1}$). The line without a symbol in both panels was arbitrarily chosen as a reference. The linear, no asymptotic line in Panel A shows the equilibrium for the case in which $h_x$ and $k$ do not depend on $S_C$ so that $S_C$ at equilibrium = $h_xR_C/k$, where $R_C$ is the residue carbon input rate (Mg C ha$^{-1}$ yr$^{-1}$). Panel A shows the equilibrium when $S_{CC}$ is doubled and Panel B shows the equilibrium $S_C$ when $h_x$ or $k$ is doubled.