Modeling coupled nitrification–denitrification in soil with an organic hotspot

The emission of nitrous oxide (N₂O) from agricultural soils to the atmosphere is a significant contributor to anthropogenic greenhouse gas emissions. The recycling of organic nitrogen (N) in manure and crop residues may result in spatiotemporal variability in N₂O production and soil efflux which is difficult to capture by process-based models. We propose a multi-species, reactive transport model to provide detailed insight into the spatiotemporal variability in nitrogen (N) transformations around such N₂O hotspots, which consists of kinetic reactions of soil respiration, nitrification, nitrifier denitrification, and denitrification represented by a system of coupled partial differential equations. The model was tested with results from an incubation experiment at two different soil moisture levels (−30 and −100 hPa) and was shown to reproduce the recorded N₂O and dinitrogen (N₂) emissions and the dynamics of important carbon (C) and N components in soil reasonably well. The simulation indicated that the four different microbial populations developed in closely connected but separate layers, with denitrifying bacteria growing within the manure-dominated zone and nitrifying bacteria in the well-aerated soil outside the manure zone and with time also within the manure layer. The modeled N₂O production within the manure zone was greatly enhanced by the combined effect of oxygen deficit, abundant carbon source, and supply of nitrogenous substrates. In the wetter soil treatment with a water potential of −30 hPa, the diffusive flux of nitrate (NO${}_{\mathrm{3}}^{-}$) across the manure–soil interface was the main source of NO${}_{\mathrm{3}}^{-}$ for denitrification in the manure zone, while at a soil water potential of −100 hPa, diffusion became less dominant and overtaken by the co-occurrence of nitrification and denitrification in the manure zone. Scenarios were analyzed where the diffusive transport of dissolved organic carbon or different mineral N species was switched off, and they showed that the simultaneous diffusion of NO${}_{\mathrm{3}}^{-}$, ammonium (NH${}_{\mathrm{4}}^{+}$), and nitrite (NO${}_{\mathrm{2}}^{-}$) was crucial to simulate the dynamics of N transformations and N₂O emissions in the model. Without considering solute diffusion in process-based N₂O models, the rapid turnover of C and N associated with organic hotspots can not be accounted for, and it may result in the underestimation of N₂O emissions from soil after manure application. The model and its parameters allow for new detailed insights into the interactions between transport and microbial transformations associated with N₂O emissions in heterogeneous soil environments.

1 Introduction

Nitrous oxide (N₂O) is a long-lived greenhouse gas (LLGHG) that accumulates in the atmosphere, accounting for about 7 % of the radiative forcing by LLGHGs (World Meteorological Organization, 2021). Globally, N₂O emissions increased from 10–12 Tg N yr⁻¹ before the industrial era (Davidson, 2009; Syakila and Kroeze, 2011) to an average of ca. 17 Tg N yr⁻¹ in the last decade (Thompson et al., 2019). Agriculture is the dominant contributor to this change, with emissions having increased from 0.3–1.0 Tg N yr⁻¹ in 1850 to 3.9–5.3 Tg N yr⁻¹ in 2010 (Davidson, 2009; Syakila and Kroeze, 2011; Thompson et al., 2019; Tian et al., 2020). Manure from animal production systems is responsible for as much as 30 %–50 % of the global N₂O emissions from agriculture (Oenema et al., 2005). The application of manure to arable land is a widely recommended practice to recycle nitrogen (N) and other nutrients for crop production. However, in a wet temperate climate it is also a large and highly variable source of N₂O emissions, of which the extent is determined by manure and soil properties and field management, with liquid manure having the greatest risk for emissions (Charles et al., 2017).

Manure has the potential to stimulate two key biochemical processes governing N₂O emissions, nitrification and denitrification. They are both regulated by multiple factors such as temperature, acidity, and availability of electron donors and acceptors, among which the interactions are highly nonlinear and difficult to predict with simple approaches (Tian et al., 2020). Such interactions may create hotspots and hot moments of N₂O emissions at small scales and have implications at landscape scales (Butterbach-Bahl et al., 2013; Groffman et al., 2009; Wagner-Riddle et al., 2020). Short-lived pulses of N₂O emissions can be induced by precipitation if anoxic soil conditions develop owing to impeded oxygen (O₂) supply from the atmosphere, provided that mineral nitrogen is present (Christensen et al., 1990a; Sexstone et al., 1985). In agricultural soils, however, manure and crop residues rich in degradable organic matter can also develop anoxic conditions by acting as a temporary sink for O₂, leading to local anoxia even in well-drained soil (Christensen et al., 1990b; Kravchenko et al., 2017) with a temporal stability that allows for microbial growth (Petersen et al., 1992, 1996). Accordingly, high spatial and temporal variations in nitrification and denitrification activity, and N₂O emissions, have been reported in manure- and plant-residue-amended soils (Kravchenko et al., 2017; Petersen et al., 1992; Taghizadeh-Toosi et al., 2021). When modeling nitrification and denitrification activity in soil, it is important to be able to include the effects of such hotspot environments.

Liquid manure (slurry) containing degradable organic carbon (C) and water is particularly prone to create anoxic hotspots upon field application. A part of the slurry infiltrates the surrounding soil in response to the soil water potential gradient (Olesen et al., 1997a; Petersen et al., 2003), but particulate matter is immobile and suspended organic particles carried with liquid manure may be trapped in the soil matrix. The extent of slurry redistribution depends on the application method, i.e., determining the extent of the manure–soil contact, and on the water retention properties of manure solids and soil (Petersen et al., 2003). Manure-saturated soil can retain a higher water content than the surrounding bulk soil for a long period of time (Olesen et al., 1997a, b). The elevated water content in conjunction with intensified O₂ consumption rates will result in gradients in the distribution of oxygen (Petersen et al., 1996; Zhu et al., 2015), which have implications for N transformations.

Ammonium (NH${}_{\mathrm{4}}^{+}$) in liquid manure infiltrating the soil will likely be adsorbed by soil particles (Olesen et al., 1997a), and the growth of nitrifying bacteria can therefore be greatly stimulated at short distance from manure-saturated volumes, where O₂ and NH${}_{\mathrm{4}}^{+}$ are both non-limiting factors (Petersen et al., 1992). Meanwhile, the lack of oxygen and higher availability of degradable carbon inside the manure-saturated zone can stimulate the activity of heterotrophic denitrifying bacteria provided that NO${}_{\mathrm{3}}^{-}$ is available. Controlled experiments (e.g., Nielsen et al., 1996; Nielsen and Revsbech, 1994) showed that denitrification rates in active organic hotspots were promoted by NO${}_{\mathrm{3}}^{-}$ from the soil, as well as NO${}_{\mathrm{3}}^{-}$ newly produced through nitrification activity, and that coupled nitrification–denitrification around oxic–anoxic interfaces can account for a large proportion of total denitrification in manure-amended soil (Meyer et al., 2002; Nielsen et al., 1996; Nielsen and Revsbech, 1994; Zhu et al., 2015). In a soil without convective water transport, NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ ions are transported by diffusion only, and the supply of soil-borne NO${}_{\mathrm{3}}^{-}$ for denitrification in manure hotspots will decline over time (Nielsen et al., 1996; Petersen et al., 1996), whereas the availability of NO${}_{\mathrm{3}}^{-}$ produced via nitrification will increase, and diffusion rates between nitrifying and denitrifying niches will be especially high around soil–manure interfaces where steep concentration gradients can develop (Petersen et al., 1992).

A close association between nitrification and denitrification activity greatly complicates the description of N transformations and N₂O emissions in models. There have been attempts to describe N₂O formation processes in soil at millimeter scale through modeling and experiments focusing on soil aggregates, where the effective denitrification rate is governed by the physical constraints on the transport of dissolved O₂ associated with aggregate size, external O₂ content, and soil respiration (Kremen et al., 2005; Schlüter et al., 2018; Smith, 1980). For example, Kremen et al. (2005) found that, with increasing aggregate radius, anaerobic conditions developed inside the aggregates, and NO${}_{\mathrm{3}}^{-}$ availability gradually became the limiting factor for denitrification. However, the implication of solute diffusion for denitrification in models of soil with macroscale heterogeneity has not been widely studied. There are numerous models to describe nitrification and denitrification processes, as well as the production of N₂O and dinitrogen (N₂). Simplified models (e.g., Conen et al., 2000; Sozanska et al., 2002) are available to use at the field or regional scale, which do not describe solute and gas movement and associated microbial processes. In process-based models (e.g., Jansson and Moon, 2001; Li et al., 2000), a relatively complete suite of biochemical processes is generally embedded to describe cycles of water, C, and N for target ecosystems. These models often include the transport of C and N species due to convective water flow but not always the diffusion process driven by concentration gradients of solutes also in the static water phase. DAISY (Hansen et al., 2012) and APSIM (Holzworth et al., 2014) are examples of models that account for diffusion, but both models do not account for the sequence of oxidation and reduction processes of mineral N, nor do they simulate O₂ directly. More commonly, process models include only C and N (or only NO${}_{\mathrm{3}}^{-}$) solute transport in soil, together with water movement, without considering the diffusion process. Example models are DayCent (Parton et al., 1998), DNDC (Li et al., 1992), DSSAT (Jones et al., 1998), PaSim (Riedo et al., 1998), and STICS (Brisson et al., 2003), i.e., models which are widely used for simulating biogeochemical cycles and related outputs (Brilli et al., 2017). This common model design for solute transport can be expected to explain field conditions where water movement (i.e., rainfall and drainage) controls the distribution of O₂, degradable C, and mineral N. However, it can lead to difficulties in reflecting the turnover of N in soil with organic hotspots where active transport of N species is important for nitrification and denitrification but uncoupled from water flow, as will be the case for extended periods after spring fertilization. Besides the solute transport module, spatial extrapolations are in most modeling approaches made using average site parameters for soil moisture and other key drivers, and the interactions among these drivers leading to hotspots either are not included or are represented at insufficient resolution (Groffman et al., 2009), making those models unsuitable to describe organic hotspots in soils.

In this work, we propose a depth- and time-varying model of N₂O emissions from an organic hotspot that builds on a system of partial differential equations (PDEs) describing the transformation of several organic or inorganic components of interest. The model presented in this study simultaneously accounts for diffusional transport of gases (i.e., O₂, CO₂, N₂O, and N₂) and solutes (i.e., dissolved C and N components), as well as biochemical processes including soil respiration, nitrification, nitrifier denitrification, and denitrification. Our main focus in this study was not to fit the model rigorously to reproduce observed N₂O emissions and associated components from manure-amended soil but rather to account for the dynamics of relevant biochemical processes around manure hotspots in soil in a way that is consistent with experimental evidence. We hypothesized that diffusion constraints on substrates and O₂ play an important role in regulating microbial activity in general and denitrification in particular. More specifically, we modeled a one-dimensional laboratory system, with a manure hotspot embedded within a repacked soil core, during a 4-week incubation period. We aimed to (1) characterize the multi-species temporal and spatial dynamics controlling C and N transformations using a modeling approach and (2) investigate the extent to which solute diffusion is important for simulating N₂O fluxes by developing and testing a series of scenarios with the model.

2 Methods and materials

2.1 Conceptual model and governing equations

The system investigated in the current study was a repacked soil core with a stagnant water phase and constant soil temperature. To model the transport and reactions of C and N components, the mass conservation equations were employed to account for time- and depth-varying concentrations of multiple components. The fate of individual components in water-filled pores was governed by biochemical reactions and diffusion, while convection–dispersion was assumed to be negligible. Gas diffusion in air, being about 4 orders of magnitude faster than that in water, was considered to be the only gas transport scheme in the model.

The mass balance equation for every component γ in the model can be described by a PDE:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\theta }(z,t)C_{\mathit{\gamma }}(z,t)}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(D_{\mathit{\gamma }}(z,t){\displaystyle \frac{\partial C_{\mathit{\gamma }}(z,t)}{\partial z}}\right)+\sum S_{\mathit{\gamma }}(z,t),\end{array}$$

where γ is dissolved organic carbon (DOC, aq), NO${}_{\mathrm{3}}^{-}$ (aq), NO${}_{\mathrm{2}}^{-}$ (aq), NH${}_{\mathrm{4}}^{+}$ (aq), carbon dioxide (CO₂, g), O₂(g), N₂O (g), N₂ (g), aerobic heterotrophs (s), ammonia-oxidizing bacteria (s), denitrifiers (s), and nitrite-oxidizing bacteria (s), in which the letters aq, g, and s indicate whether component γ is associated with the aqueous, the gaseous, or the solid phase in the model; z and t represent soil depth and time, respectively; C_γ is the concentration of substrate γ; θ represents either the water-filled (θ_aq) or the air-filled porosity (θ_g) depending on the phase of C_γ; and D_γ is the effective diffusion coefficient for component γ (see Sect. S7.1 for details). For the bacterial populations, the diffusion term is excluded in the equation. S_γ is the source or sink term of a γ component caused by biochemical reactions, which will be depicted in Sect. 2.2. As a component may participate in several reactions, the 12 PDEs for the above components are coupled by the reaction terms and must be solved simultaneously.

For the cationic species NH${}_{\mathrm{4}}^{+}$ in the model, the adsorbed and dissolved forms were assumed to be at equilibrium and follow the Freundlich model as defined below (Olesen et al., 1999):

$$\begin{array}{}\text{(2)}& C_{\mathrm{s}}=K_{\mathrm{F}}\left(C_{\mathrm{w}}^N\right),\end{array}$$

where C_s is the concentration of adsorbed NH${}_{\mathrm{4}}^{+}$, C_w is the concentration of dissolved NH${}_{\mathrm{4}}^{+}$, K_F is the Freundlich distribution coefficient, and N is the dimensionless Freundlich isotherm exponent.

A retardation factor $R_{{\mathrm{NH}}_{\mathrm{4}}^{+}}$, describing the effect of cation adsorption to soil particles causing a transport delay, was multiplied with the left-hand side of Eq. (1):

$$\begin{array}{}\text{(3)}& R_{{\mathrm{NH}}_{\mathrm{4}}^{+}}=\mathrm{1}+\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{b}}}{{\mathit{\theta }}_{\mathrm{w}}}}\right)K_{\mathrm{d}}=\mathrm{1}+\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{b}}}{{\mathit{\theta }}_{\mathrm{w}}}}\right)K_{\mathrm{F}}N{C}_{\mathrm{w}}^{N-\mathrm{1}},\end{array}$$

where ρ_b is the soil bulk density, and K_d is the distribution coefficient of NH${}_{\mathrm{4}}^{+}$ between soil solids and water.

Figure 1 Schematic diagrams of (a) the simulated 10 cm soil core with slurry application and two surfaces connected to the ambient air, (b) conceptualized initial profiles of some key drivers of N cycling in the system, and (c) biochemical processes included in the model where the notation (−) and (+) indicate the negative or positive response of the process to O₂ or DOC. ND: nitrifier denitrification; AER: aerobic heterotrophs; AOB: ammonia-oxidizing bacteria; NOB: nitrite-oxidizing bacteria; and DEN: denitrifiers.

2.2 Reaction processes

As presented in Fig. 1, the mathematical model developed in this study integrated relevant functional groups of microorganisms in a soil–manure system including aerobic heterotrophs, ammonia-oxidizing bacteria (AOB), nitrite-oxidizing bacteria (NOB), and denitrifying bacteria. The model considered the following processes: aerobic respiration, nitrification, nitrifier denitrification, and denitrification.

In aerobic decomposition, living heterotrophic microorganisms, in the presence of oxygen, feed upon the organic carbon which serves as a source of energy and is respired as CO₂:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{CO}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\text{(5)}& {\displaystyle }\begin{array}{rl}{S}_{{\mathrm{CO}}_{\mathrm{2}}+,\mathrm{r}}& ={\mathit{\rho }}_{\mathrm{b}}{\mathit{\mu }}_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{r}}{B}_{\mathrm{AER}}{\displaystyle \frac{\left[\mathrm{C}\right]}{\left[\mathrm{C}\right]+k{M}_{\mathrm{C}-{\mathrm{CO}}_{\mathrm{2}}}}}\\ & \times {\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+k{M}_{{\mathrm{O}}_{\mathrm{2}}-{\mathrm{CO}}_{\mathrm{2}}}}},\end{array}\end{array}$$

where ${\mathit{\mu }}_{{\mathrm{CO}}_{\mathrm{2}},\mathrm{r}}$ (mmol CO₂ g⁻¹ biomass d⁻¹) is the maximum reaction rate regarding microbial biomass, B_AER (g biomass g⁻¹ dw – dry weight) is the total biomass of aerobic heterotrophs in soil, ρ_b (g L⁻¹) is the bulk soil density, and kM (mmol L⁻¹) is the half-saturation constant of substrates. Further, [C] and [O₂] represent available concentrations of carbon and oxygen at the reactive sites of the enzyme, which will be described in the following. In the subscript of the reaction velocity, S, we use the sign “+” or “−” to indicate production or consumption of the component, and we use the letters “r”, “n”, “nd”, and “dn” to represent the aerobic respiration, nitrification, nitrifier denitrification, and denitrification processes; the same applies to the equations below.

Three pools of organic carbon were considered: immobile C associated with soil organic matter (immobile SOC), immobile C associated with particulate manure solids (manure POC), and DOC. We assumed that only DOC was the substrate for aerobic respiration, and first-order kinetics were used to describe the conversion from immobile SOC and manure POC to DOC:

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{SOC}-}={\displaystyle \frac{\partial \mathrm{SOC}(z,t)}{\partial t}}=-{\mathit{\alpha }}_{\mathrm{SOC}}\mathrm{SOC},\text{(7)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{POC}-}={\displaystyle \frac{\partial \mathrm{POC}(z,t)}{\partial t}}=-{\mathit{\alpha }}_{\mathrm{POC}}\mathrm{POC},\end{array}$$

where S_SOC− and S_POC− (g C g⁻¹ dw d⁻¹) are the rates of SOC and POC conversion to DOC, and α (d⁻¹) is the conversion rate.

Other main biochemical reactions in the model, including respiration, nitrification, nitrifier denitrification, and denitrification, were all assumed to follow Michaelis–Menten kinetics. In the nitrification process, the model assumes that NH${}_{\mathrm{4}}^{+}$ is oxidized directly to NO${}_{\mathrm{2}}^{-}$ by AOB and subsequently to NO${}_{\mathrm{3}}^{-}$ by NOB (Eqs. 8 and 9). This setup thus does not include the intermediates NH₂OH and NO of ammonia oxidation, in contrast to the explicit description found in Chang et al. (2022) and Chen et al. (2019). N₂O is generated as a by-product of incomplete oxidation of NH${}_{\mathrm{4}}^{+}$ to NO${}_{\mathrm{2}}^{-}$ by AOB (Eq. 10, adapted from Eqs. 1–3 in Chang et al., 2022). The consumption of O₂ by nitrifiers is included in the source and sink terms. AOB is also responsible for nitrifier denitrification where NO${}_{\mathrm{2}}^{-}$ is reduced to N₂O (Eq. 11, adapted from Eqs. 1 and 6 in Chang et al., 2022). A two-sided effect of O₂, both promotion and inhibition (Eq. 15), is included on this pathway, wherein O₂ is required to support NH${}_{\mathrm{4}}^{+}$ oxidation while also inhibiting the reduction of NO${}_{\mathrm{2}}^{-}$ to N₂O (Wrage et al., 2001).

The governing equations and reaction velocities of nitrification and nitrifier denitrification are

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathrm{NH}}_{\mathrm{4}}^{+}+\mathrm{1.5}{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{NO}}_{\mathrm{2}}^{-}+\mathrm{2}{\mathrm{H}}^{+}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\text{(9)}& {\displaystyle }{\mathrm{NO}}_{\mathrm{2}}^{-}+\mathrm{0.5}{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{NO}}_{\mathrm{3}}^{-},\text{(10)}& {\displaystyle }\mathrm{2.5}{\mathrm{NH}}_{\mathrm{4}}^{+}+\mathrm{2.75}{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{N}}_{\mathrm{2}}\mathrm{O}+\mathrm{0.5}{\mathrm{NO}}_{\mathrm{2}}^{-}+\mathrm{3}{\mathrm{H}}^{+}+\mathrm{3.5}{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\text{(11)}& {\displaystyle }{\mathrm{NO}}_{\mathrm{2}}^{-}+{\mathrm{NH}}_{\mathrm{4}}^{+}+\mathrm{0.5}{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{N}}_{\mathrm{2}}\mathrm{O}+\mathrm{2}{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\text{(12)}& {\displaystyle }\begin{array}{rl}{S}_{{\mathrm{NO}}_{\mathrm{2}}^{-}+,\mathrm{n}}& ={\mathit{\rho }}_{\mathrm{b}}{\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{2}}^{-},\mathrm{n}}{B}_{\mathrm{AOB}}{\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]}{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]+k{M}_{{\mathrm{NH}}_{\mathrm{4}}^{+}-{\mathrm{NO}}_{\mathrm{2}}^{-}}}}\\ & \times {\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+k{M}_{{\mathrm{O}}_{\mathrm{2}}-{\mathrm{NO}}_{\mathrm{2}}^{-}}}},\end{array}\text{(13)}& {\displaystyle }\begin{array}{rl}{S}_{{\mathrm{NO}}_{\mathrm{3}}^{-}+,\mathrm{n}}& ={\mathit{\rho }}_{\mathrm{b}}{\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}^{-},\mathrm{n}}{B}_{\mathrm{NOB}}{\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{2}}^{-}\right]}{\left[{\mathrm{NO}}_{\mathrm{2}}^{-}\right]+k{M}_{{\mathrm{NO}}_{\mathrm{2}}^{-}-{\mathrm{NO}}_{\mathrm{3}}^{-}}}}\\ & \times {\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+k{M}_{{\mathrm{O}}_{\mathrm{2}}-{\mathrm{NO}}_{\mathrm{3}}^{-}}}},\end{array}\text{(14)}& {\displaystyle }\begin{array}{rl}{S}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}+,\mathrm{n}}& ={\mathit{\rho }}_{\mathrm{b}}{\mathit{\mu }}_{{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}^{+},\mathrm{n}}{B}_{\mathrm{AOB}}{\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]}{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]+k{M}_{{\mathrm{NH}}_{\mathrm{4}}^{+}-{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}}\\ & \times {\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+k{M}_{{\mathrm{O}}_{\mathrm{2}}-{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\mathrm{n}}}},\end{array}\text{(15)}& {\displaystyle }\begin{array}{rl}{S}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}+,\mathrm{nd}}& ={\mathit{\rho }}_{\mathrm{b}}{\mathit{\mu }}_{{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}^{+},\mathrm{nd}}{B}_{\mathrm{AOB}}{\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{2}}^{-}\right]}{\left[N{O}_{\mathrm{2}}^{-}\right]+k{M}_{{\mathrm{NO}}_{\mathrm{2}}^{-}-{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}}\\ & \times {\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]}{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]+k{M}_{{{\mathrm{NH}}_{\mathrm{4}}}^{+}-{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}}\\ & \times {\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+k{M}_{{\mathrm{O}}_{\mathrm{2}}-{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\mathrm{nd}}}}\times {\displaystyle \frac{k{I}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+k{I}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\mathrm{nd}}}},\end{array}\end{array}$$

where μ (mmol g⁻¹ biomass d⁻¹) is the maximum reaction rate of microbial biomass in individual steps, B_AOB and B_NOB (g biomass g⁻¹ dw) are the total microbial biomass of AOB and NOB in soil, and kI (mmol L⁻¹) is the inhibition constant of O₂.

In the modeling of denitrification, denitrifiers use the carbon source (DOC) to gain energy and reduce NO${}_{\mathrm{3}}^{-}$ stepwise to NO${}_{\mathrm{2}}^{-}$, N₂O, and N₂ (Eqs. 16–18). The stimulation of denitrification by DOC, as well as the inhibition by O₂, is considered in the model. In each modeled step of denitrification, the reduction in nitrogenous oxides is accompanied by anaerobic respiration, whereby CO₂ is produced.

The governing equations and reaction velocities of denitrification are

$$\begin{array}{}\text{(16)}& {\displaystyle }{\mathrm{NO}}_{\mathrm{3}}^{-}+\mathrm{0.5}{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}\to {\mathrm{NO}}_{\mathrm{2}}^{-}+\mathrm{0.5}{\mathrm{CO}}_{\mathrm{2}}+\mathrm{0.5}{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\text{(17)}& {\displaystyle }\mathrm{2}{\mathrm{NO}}_{\mathrm{2}}^{-}+{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}+\mathrm{2}{\mathrm{H}}^{+}\to {\mathrm{N}}_{\mathrm{2}}\mathrm{O}+{\mathrm{CO}}_{\mathrm{2}}+\mathrm{2}{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\text{(18)}& {\displaystyle }{\mathrm{N}}_{\mathrm{2}}\mathrm{O}+\mathrm{0.5}{\mathrm{CH}}_{\mathrm{2}}\mathrm{O}\to {\mathrm{N}}_{\mathrm{2}}+\mathrm{0.5}{\mathrm{CO}}_{\mathrm{2}}+\mathrm{0.5}{\mathrm{H}}_{\mathrm{2}}\mathrm{O},\text{(19)}& {\displaystyle }\begin{array}{rl}{S}_{{\mathrm{NO}}_{\mathrm{2}}^{-}+,\mathrm{dn}}& ={\mathit{\rho }}_{\mathrm{b}}{\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{2}}^{-},\mathrm{dn}}{B}_{\mathrm{DEN}}{\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]}{\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]+k{M}_{{\mathrm{NO}}_{\mathrm{3}}^{-}-{\mathrm{NO}}_{\mathrm{2}}^{-}}}}\\ & \times {\displaystyle \frac{\left[\mathrm{C}\right]}{\left[\mathrm{C}\right]+k{M}_{\mathrm{C}-{\mathrm{NO}}_{\mathrm{2}}^{-}}}}\times {\displaystyle \frac{k{I}_{{\mathrm{NO}}_{\mathrm{2}}^{-}}}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+k{I}_{{\mathrm{NO}}_{\mathrm{2}}^{-},\mathrm{dn}}}},\end{array}\text{(20)}& {\displaystyle }\begin{array}{rl}{S}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}+,\mathrm{dn}}& ={\mathit{\rho }}_{\mathrm{b}}{\mathit{\mu }}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\mathrm{dn}}{B}_{\mathrm{DEN}}{\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{2}}^{-}\right]}{\left[{\mathrm{NO}}_{\mathrm{2}}^{-}\right]+k{M}_{{\mathrm{NO}}_{\mathrm{2}}^{-}-{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}}\\ & \times {\displaystyle \frac{\left[\mathrm{C}\right]}{\left[\mathrm{C}\right]+k{M}_{\mathrm{C}-{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}}\times {\displaystyle \frac{k{I}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+k{I}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\mathrm{dn}}}},\end{array}\text{(21)}& {\displaystyle }\begin{array}{rl}{S}_{{\mathrm{N}}_{\mathrm{2}}+,\mathrm{dn}}& ={\mathit{\rho }}_{\mathrm{b}}{\mathit{\mu }}_{{\mathrm{N}}_{\mathrm{2}},\mathrm{dn}}{B}_{\mathrm{DEN}}{\displaystyle \frac{\left[{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]}{\left[{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]+k{M}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}}\\ & \times {\displaystyle \frac{\left[\mathrm{C}\right]}{\left[\mathrm{C}\right]+k{M}_{\mathrm{C}-{\mathrm{N}}_{\mathrm{2}}}}}\times {\displaystyle \frac{k{I}_{{\mathrm{N}}_{\mathrm{2}}}}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+k{I}_{{\mathrm{N}}_{\mathrm{2}},\mathrm{dn}}}},\end{array}\text{(22)}& {\displaystyle }{S}_{{\mathrm{CO}}_{\mathrm{2}}+,\mathrm{dn}}=\mathrm{0.5}{S}_{{\mathrm{NO}}_{\mathrm{2}-}+,\mathrm{dn}}+S_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}+,\mathrm{dn}}+\mathrm{0.5}{S}_{{\mathrm{N}}_{\mathrm{2}}+,\mathrm{dn}}.\end{array}$$

We assumed that substrate concentrations in the reaction velocity equations represent the available substrates at the reactive sites of enzymes and that the availability of substrates is affected by the diffusion of substrates through water films, following the work by Davidson et al. (2012). Therefore, since a decline in soil moisture is accompanied by more disconnected water films, the concentrations of dissolved substrate at the site of enzyme reaction will decline and the available concentration of gaseous reactants will increase. The available concentrations [C_aq] in the aqueous phase at the enzyme reaction site are calculated based on a dimensionless diffusivity D_aq and water content, and similarly gaseous concentrations D_g and air content are used for calculating the available concentration [C_g] in the gas phase:

$$\begin{array}{}\text{(23)}& {\displaystyle }& {\displaystyle }\left[C_{\mathrm{aq}}\right]=C_{\mathrm{aq}}\times D_{\mathrm{aq}}\times {\mathit{\theta }}_{\mathrm{aq}}^{\mathrm{3}},\text{(24)}& {\displaystyle }& {\displaystyle }\left[C_{\mathrm{g}}\right]=C_{\mathrm{g}}\times D_{\mathrm{g}}\times {\mathit{\theta }}_{\mathrm{g}}^{\mathrm{4}/\mathrm{3}},\end{array}$$

where C_aq and C_g indicate the actual concentrations of dissolved substrates in water and concentrations in air. D_aq and D_g are unitless diffusion coefficients of solute in water and air, respectively (Davidson et al., 2012). The value of D_aq was determined by assuming the extreme condition that [C_aq]=C_aq for saturated soil; i.e., all of the soluble substrate is available at the reaction site under this condition (Papendick and Campbell, 2015). The value of D_g is determined by another assumed extreme condition that all of the gas is available at the reaction site in completely dry soil (Millington, 1959). The effective diffusivity used in Fick's law and in Eq. (1), originated from the same sources we used for the dimensionless diffusivities. See the Supplement (Sects. S6.2 and S7.1) for more details.

In the defined model, the total biomass of each microbe, B_microbe, is composed of two parts: an existing biomass associated with the bulk soil (B_base,microbe) and new biomass produced during the simulation (B_new,microbe). The growth of microbial biomass in the model is proportional to the rate of substrate consumption following the Monod equation and the microbial decay following first-order kinetics:

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\left[B_{\mathrm{new},\mathrm{AER}}\right]}{\mathrm{d}t}}=-y_{\mathrm{AER}}/{\mathit{\rho }}_{\mathrm{b}}{f}_{\mathrm{Cbio}}{S}_{\mathrm{DOC}-,\mathrm{r}}-a_{\mathrm{AER}}{B}_{\mathrm{new},\mathrm{AER}},\text{(26)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\left[B_{\mathrm{new},\mathrm{AOB}}\right]}{\mathrm{d}t}}=-y_{\mathrm{AOB}}/{\mathit{\rho }}_{\mathrm{b}}{f}_{\mathrm{Nbio}}{S}_{{\mathrm{NH}}_{\mathrm{4}}^{+}-,\mathrm{n}}-a_{\mathrm{AOB}}{B}_{\mathrm{new},\mathrm{AOB}},\text{(27)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\left[B_{\mathrm{new},\mathrm{NOB}}\right]}{\mathrm{d}t}}=-y_{\mathrm{NOB}}/{\mathit{\rho }}_{\mathrm{b}}{f}_{\mathrm{Nbio}}{S}_{{\mathrm{NO}}_{\mathrm{2}}^{-}-,\mathrm{n}}-a_{\mathrm{NOB}}{B}_{\mathrm{new},\mathrm{NOB}},\text{(28)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\left[B_{new,DEN}\right]}{\mathrm{d}t}}=-y_{\mathrm{DEN}}/{\mathit{\rho }}_{\mathrm{b}}{f}_{\mathrm{Cbio}}{S}_{\mathrm{DOC}-,\mathrm{dn}}-a_{\mathrm{DEN}}{B}_{\mathrm{new},\mathrm{DEN}},\end{array}$$

where y is the yield coefficient (g C g⁻¹ C or g N g⁻¹ N), f_Cbio (g C g⁻¹ biomass) is the C content of microbial biomass, f_Nbio (g N g⁻¹ biomass) is the N content of the microbial biomass, S (g C L⁻¹ soil d⁻¹ or g N L⁻¹ soil d⁻¹) is the rate of substrate consumption, and a (d⁻¹) is the microbial decay rate.

2.3 Incubation experiment

The proposed model was first benchmarked against a lab incubation experiment (Kolstad et al., 2023) investigating the emission of N₂O, the total denitrification products, and the spatial distribution of inorganic N after manure incorporation under contrasting levels of soil water content.

In brief, the incubation experiment was set up as follows. Acrylic cylinders (height 5 cm, diameter 8.4 cm) were packed with partially dried and sieved (<6 mm) sandy loam soil collected in November from the plow layer (0–25 cm) of a long-term field experiment at the Foulumgaard research facility, Denmark. On day 0 of the experiment, the injection of liquid cattle manure was simulated by applying the manure to the surface of two uniformly packed soil cores which, after the liquid had infiltrated the two cores, were combined with the manure-saturated zone embedded at the center of the now 10 cm high soil cores. The manure application rate was 0.40 g cm⁻² (corresponding to 40 t ha⁻¹), and basic properties of the manure are listed in Table S1 in the Supplement. Two soil–manure treatments having a bulk density of 1.4 g cm⁻³ and final soil moisture levels at either −30 or −100 hPa (corresponding to a water-filled pore space of 85 % and 70 %, respectively) were used to benchmark the model. For each water potential there was a control treatment, where the soil was packed as described above to the same moisture levels, but slurry was not added. The samples were prepared in triplicate and incubated for 28 d, with gas sampling on days 1, 3, 7, 10, 14, 21, and 28 to determine fluxes of N₂O and CO₂ using gas chromatography (GC). Separate samples were sectioned on days 1, 14, and 28 to determine the distribution of NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$, NO${}_{\mathrm{2}}^{-}$, and water, while loss on ignition (LOI) was analyzed as a measure of organic matter content on day 1 only. The samples were sliced at 0, 2, 6, 10, 14, 18, and 30 mm distance from the center to either side of the manure layer. The samples for sectioning and analysis of soil variables were prepared in duplicate. The soil NO${}_{\mathrm{3}}^{-}$ pool was labeled with ¹⁵N to allow for both N₂O and N₂ fluxes to be determined by isotope ratio mass spectrometry (IRMS). Two different solutions of K¹⁵NO₃ were prepared for addition during packing to ensure the same final concentration at both water potentials. Extra gas samples were collected on days 1, 3, 7, and 14 for IRMS analysis.

The experimental results from the two treatments were used for comparison with model outputs, including the daily fluxes of three gases (N₂O, N₂, and CO₂) and measurements of LOI, NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$, and NO${}_{\mathrm{2}}^{-}$. The N₂ fluxes were calculated from N₂O fluxes determined by GC and ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{N}}_{\mathrm{2}}$ ratios determined by IRMS. When N₂O or N₂ fluxes were below detection limits, ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{N}}_{\mathrm{2}}$ were considered invalid, and these sampling values were then treated as missing data.

2.4 Model implementation

2.4.1 Initial conditions

The coupled PDEs, as described by Eq. (1) and associated text, had to obey relevant initial and boundary conditions to complete the model setup for the realistic system investigated. The buildup of the soil–manure system started with the application of manure and redistribution of liquid manure into a larger soil volume until the final water potential in the core reached equilibrium at −30 or −100 hPa, a process with a duration of a few hours. In conceptualizing the experimental setup in a model system, we took conditions from when the slurry redistribution had ceased as the starting point of the model. The soil water distribution was largely constant between day 1 and day 28 (Fig. S4), and therefore Eq. (1), where water convection was excluded, was suitable to use in the model system.

Given the lack of information about the soil–manure system prior to the first soil sampling on day 1, we made a few assumptions to set the initial conditions of the model system:

• i.

  Liquid manure was redistributed within a limited “manure zone” and displaced soil pore water resulting in elevated water content compared to before application.

• ii.

  During slurry redistribution, the fate of dissolved components including DOC, NH${}_{\mathrm{4}}^{+}$, and NO${}_{\mathrm{3}}^{-}$ in the slurry was dominated by water flow, while biochemical reactions and sorption were neglected, entailing that the components were uniformly concentrated in the limited zone with a constant total amount.

• iii.

  When slurry redistribution was complete, the concentration of oxygen in the 2 mm center of the limited manure zone, wherein particulate manure solids were concentrated, was zero.

The soil core in the model had a depth of 10 cm from the top (z=0 cm) to the bottom surface (z=10 cm), and the center of manure application was at 5 cm depth. The soil core in the model had the same length as in the experiment with the manure at its center. The vertical extent of the manure zone was determined to be 8 mm, from 4.6 to 5.4 cm depth, in which the amount of water on day 1 corresponded to 20.2 and 17.5 g in the −30 and −100 hPa treatments, respectively, close to the 20.9 g of water in the applied slurry. We used the recorded profiles of water content on day 1 as the constant water condition in the modeled soil core.

Dissolved and particulate C and N

Contributions of DOC, NH${}_{\mathrm{4}}^{+}$, and NO${}_{\mathrm{3}}^{-}$ from both soil and manure were considered in the initial conditions. Soil properties of the control soil on day 1 in the two layers closest to the surfaces, i.e., 0–0.01 m and 0.09–0.1 m, were averaged (Table S2). Soil organic carbon (SOC) was estimated from the relationship $\mathrm{SOC}=\mathrm{0.39}\cdot \mathrm{LOI}-\mathrm{0.28}$ (Jensen et al., 2018). DOC concentrations in the Foulum soil in winter are fairly constant at 20–25 mg C L⁻¹ (Gjettermann et al., 2008), and a conversion factor of $\mathrm{3.5}\times {\mathrm{10}}^{-\mathrm{4}}$ between DOC and SOC was used to make estimates of DOC consistent with this reported range. Concentrations of dissolved NH${}_{\mathrm{4}}^{+}$ were calculated from the total content of NH${}_{\mathrm{4}}^{+}$ and the abovementioned Freundlich model. NO${}_{\mathrm{3}}^{-}$ was assumed to only exist in the water phase. See the Supplement (Sect. S6.1 and S6.3) for details.

Regarding the contributions from the manure, we considered that the total amount of each dissolved component in the applied manure was uniformly distributed in the manure zone of 0.046–0.054 m. Based on the manure application rate and content of volatile solids (VSs) in the manure, we estimated that the total organic carbon (TOC) accounted for a fraction of 0.42 of VSs (Petersen et al., 2016). We also assumed that the fraction of DOC in TOC was 0.5 (Petersen et al., 1996, 2016), and hence the amount of DOC in the applied manure was estimated at 30.9 g C m⁻². Particulate organic carbon (POC) in the manure was then defined as (TOC − DOC); assuming that manure POC was not mobile, this fraction was considered to concentrate within 1 mm from the depth of application, i.e., from 0.049 to 0.051 m depth. Dissolved NH${}_{\mathrm{4}}^{+}$ from manure was estimated from total NH${}_{\mathrm{4}}^{+}$ as above. As there was no detected NO${}_{\mathrm{3}}^{-}$ in the manure, the initial NO${}_{\mathrm{3}}^{-}$ content in the 8 mm manure zone was set to be zero. The initial NO${}_{\mathrm{2}}^{-}$ content was considered to be zero in the entire soil–manure profile. We did not include the mineralization of organic N from POC in the model, as previous studies have associated high $\mathrm{C}/\mathrm{N}$ with reduced N mineralization in slurry (Barbarika et al., 1985), whereas we considered that mineral N could be immobilized by microbes, i.e., AOB and NOB, and become part of the microbial biomass pool.

The initial content of each component in the 14 sampling layers was thus obtained by summing the contributions of manure and soil.

Gases

We assumed that the initial concentrations of N₂, O₂, CO₂, and N₂O corresponded to ambient air, i.e., 0.78, 0.21, $\mathrm{4.1}\times {\mathrm{10}}^{-\mathrm{4}}$, and $\mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{7}}$ atm, respectively (World Meteorological Organization, 2021); according to the ideal gas equation, their concentrations in soil air were calculated to be 33.0, 8.9, 0.017, and $\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{5}}$ mmol L⁻¹, respectively. One exception was O₂ at the center of the manure zone (i.e., 0.049–0.051 m) where anoxic conditions were assumed and the concentration was zero.

Microbial populations

The initial biomasses of denitrifiers, AOB, and NOB in the bulk soil were all set to be 12.7 mg kg⁻¹ soil, referring to the measurement in an arable soil by Khalil et al. (2005) using chloroform fumigation–extraction and microbial enumeration, while the initial biomass of aerobic heterotrophs was assumed to be 10 times higher. The new biomasses of aerobic heterotrophs, AOB, and NOB initially added with manure were assumed to be zero, as we considered the concentrations of these aerobic organisms to be negligible in the anaerobic liquid manure. In contrast, we expected denitrifiers to survive better in the manure, and the denitrifier biomass at the center of the applied manure (0.049–0.051 m) was assumed (including the contribution from soil) to be 10 times the base value in soil and thus similar to the base value of aerobic heterotrophs.

For all of the above components to be solved with Eq. (1), their initial conditions over the modeled soil core were constructed as continuous profiles by interpolating the middle points of all sampling layers and extrapolating to the two boundaries (see Figs. S1 and S2 in the Supplement).

2.4.2 Boundary conditions

For components dissolved in the water phase and microbial populations, we specified no flux at the two boundaries:

$${\displaystyle \frac{\partial C_{\mathit{\gamma }}(z,t)}{\partial z}}=\mathrm{0},z=\mathrm{0}\mathrm{and}z=\mathrm{0.1}.$$

For the concentration of each gaseous component, we specified a constant value for the two boundaries, which was the gas concentration in the ambient air as described in Sect. 2.4.1:

$$C_{\mathit{\gamma }}(z,t)=C_{\mathrm{atm}}\left(t\right),z=\mathrm{0}\mathrm{and}z=\mathrm{0.1}.$$

2.4.3 Computation

The system of coupled equations for eight chemical components and four microbial populations was numerically solved using the partial differential equations solver pdepe in MATLAB. The first-order kinetics of SOC and POC were analytically solved and incorporated into the source terms of DOC production to reduce the complexity of the model system. The solver internally uses a dynamic numerical time stepping, and the output time step was set to 6 h. A non-uniform mesh was applied for the discretization of the 0.1 m one-dimensional model domain. The mesh size was 1 mm in the area away from the manure-applied zone (z<0.04 m, z>0.06 m) and 0.1 mm in the interval of 0.04–0.06 m to obtain stable numerical solutions at different parameter values in preliminary tests.

For each component, Eq. (1) was valid when the corresponding θ was above zero. When θ was equal to zero, there is an inequality between the zero left-hand side and the non-zero right-hand side with the reaction term. In simulating the −30 hPa treatment, a narrow 2 mm saturated zone existed at the center of the core, and to avoid the invalid model domain and not to disturb the rest of model system, a small air fraction of 0.01 was allowed for in the saturated zone, corresponding to a few drained macropores. The relevant discussion on model uncertainty is included in Sect. 4.3.

The parameter fitting was performed by a combination of algorithm and manual adjustment. We used the surrogate algorithm as the optimizer, as it is suitable for evaluating time-consuming models such as the complex PDE system. The objective function used for optimization was the sum of relative root mean squared errors (rRMSEs) of simulated variables from the two treatments:

$$\begin{array}{}\text{(29)}& \mathrm{obj}=\sum _{j=\mathrm{1}}^N{w}_j{\displaystyle \frac{\sqrt{\frac{\mathrm{1}}{n_j}\sum _{i=\mathrm{1}}^{n_j}{\left(y_{i,j}-{\widehat{y}}_{i,j}\right)}^{\mathrm{2}}}}{\frac{\mathrm{1}}{n_j}\sum _{i=\mathrm{1}}^{n_j}{y}_{i,j}}},\end{array}$$

where obj is the value of objective function, y_i,j is the ith observation for variable j, ${\widehat{y}}_{i,j}$ is the ith simulation for variable j, n_j is the number of observations for variable j, w_j is the weight of variable j, and N is the number of variables.

Parameter values of the C module, k_{C_CO2_r}, k_{O2_CO2_r}, and u_{C_CO2_r}, were optimized by fitting model simulations with measured CO₂ fluxes in the −30 and −100 hPa treatments ($w_{{\mathrm{CO}}_{\mathrm{2}}}=\mathrm{1}$ for both treatments), while the rest of the parameters were fixed to the starting values, as we considered the majority of CO₂ fluxes to come from aerobic respiration. The values of parameters in the N module were then optimized by fitting the N-related variables with measured data. The N-related variables here included N₂O and N₂ fluxes, as well as NO${}_{\mathrm{3}}^{-}$ and NO${}_{\mathrm{2}}^{-}$ content. DOC and NH${}_{\mathrm{4}}^{+}$ were not included as the model-data error was not much affected by parameter fitting in preliminary tests. When presenting the data used to compare with model results in the Results section, we used the labels “measured” and “estimated” to distinguish between the variables obtained directly from measurements (CO₂ fluxes, N₂O fluxes, NO${}_{\mathrm{3}}^{-}$ content, and NO${}_{\mathrm{2}}^{-}$ content, where the latter two were considered to only exist in dissolved form) and the variables (N₂ fluxes, DOC, and dissolved NH${}_{\mathrm{4}}^{+}$) estimated from relevant experimental data and assumptions stated earlier.

To ensure good agreements with the temporal and spatial dynamics of multiple variables in the experiment, the selection of precise and accurate weights in objective functions was often difficult as small perturbations of the weights can lead to quite different solutions (Konak et al., 2006). By trial and error, we assigned the weights of 0.2, 0.05, 0.012, and 0.001 to the rRMSEs regarding N₂O, N₂, NO${}_{\mathrm{2}}^{-}$, and NO${}_{\mathrm{3}}^{-}$ in the −30 hPa treatment, as well as weights of 0.05, 0.05, 0.001, and 0.001 to the corresponding rRMSEs in the −100 hPa treatment. Model errors for simulated mineral N were prone to be much higher than for gas emissions, and we therefore assigned them smaller weights to avoid overfitting mineral N accompanied by the poor prediction of gas fluxes while minimizing the objective function. Since the peak N₂O flux in the −30 hPa treatment was of particular interest, we included also a term representing peak flux error, (${\mathrm{peak}}_{\mathrm{obs}}-{\mathrm{peak}}_{\mathrm{sim}})/{\mathrm{peak}}_{\mathrm{obs}}$, in the objective function to ensure this interest was met. Here, peak_obs and peak_sim indicate the maximum daily N₂O flux found in the experiment and the model, respectively. After checking the optimized model outputs, which in general aligned well with the measurements, we manually adjusted two parameters (μ_{NO2_dn} and μ_{N2O_nd}) by trials to make the temporal trend of NO${}_{\mathrm{2}}^{-}$ profiles more consistent with the measurements. Model input parameters are described in Sect. S7.

2.4.4 Flux calculation

The efflux of gases from the soil cores was calculated according to Fick's first law, where the concentration gradient between the ambient air and the outmost soil layers at the top and bottom was the driving force of gas flux:

$$\begin{array}{}\text{(30)}& \begin{array}{rl}J\left(t\right)& ={\displaystyle \frac{D_{\mathrm{eff},\mathrm{g}}\left(z=z_{\mathrm{1}},t\right)}{\mathrm{\Delta }z}}\left(C\left(z=z_{\mathrm{1}},t\right)-C_{\mathrm{atm}}\right)\\ & +{\displaystyle \frac{D_{\mathrm{eff},\mathrm{g}}\left(z=z_{\mathrm{2}},t\right)}{\mathrm{\Delta }z}}\left(C\left(z=z_{\mathrm{2}},t\right)-C_{\mathrm{atm}}\right),\end{array}\end{array}$$

where J is the gas fluxes (mol m⁻² d⁻¹), D_eff,g (m² d⁻¹) is the effective gas diffusion coefficients at the two depths closest to the borders in the discretization of the model where z₁=0.001 m and z₂=0.099 m, and Δz is the depth interval of 0.001 m.

2.4.5 Scenario tests

We took the original setup of the model where the diffusion process of all the solutes was included as the baseline scenario. Based on the calibrated model, we simulated four scenarios in which solute diffusion in the soil–manure system was eliminated to different degrees. Diffusive fluxes of different solutes were selectively turned off in Scenarios 1–4 as shown in Table 1, where Scenario 1 did not allow any solute diffusion; Scenario 2 allowed only NO${}_{\mathrm{3}}^{-}$ diffusion; Scenario 3 allowed only NH${}_{\mathrm{4}}^{+}$ diffusion; and Scenario 4 allowed diffusion of NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$, and NO${}_{\mathrm{2}}^{-}$ but not DOC. In Eq. (1), turning off the diffusive flux of component γ meant that the effective diffusion coefficient D_γ was set to zero.

Table 1 Tested scenarios of the influence of solute diffusion.

3 Results

3.1 Testing model predictive capability with experimental data

Emissions of three gases, i.e., N₂O, N₂, and CO₂, during 4 weeks in the experiment were compared with model simulations (Fig. 2). In the −30 hPa treatment with manure application, the model showed a peak N₂O flux on day 3 which was comparable with the peak value of 2880 µg N m⁻² h⁻¹ observed on day 1 but with a delay of 2 d. In the −100 hPa treatment, the simulated N₂O fluxes in the first 3 d were between 500 and 1500 µg N m⁻² h⁻¹, which was much higher than the recorded fluxes with a maximum of 130 µg N m⁻² h⁻¹. After the first week the simulation results corresponded well with the tail of fluxes recorded in the experiment, and small negative values were obtained in the model for the −100 hPa treatment. The simulated N₂ fluxes agreed with the magnitude of N₂ fluxes estimated from ¹⁵N labeling in the experiment and were approximately 1 order of magnitude higher than the observed N₂O fluxes in both treatments. Specifically, the temporal trend of N₂ fluxes in the −30 hPa treatment was simulated well despite the limited experimental data available. The simulated CO₂ fluxes were in general comparable to those observed, especially in the −100 hPa treatment with an early CO₂ peak around day 1 in both model and experiment. In the experiment, CO₂ fluxes in the −30 hPa treatment were generally lower than in the −100 hPa treatment, and model results showed the same pattern. A secondary increase in CO₂ flux, observed during the second week in the −30 hPa treatment, was not shown in the simulation.

Figure 2 Measured and modeled N₂O, N₂, and CO₂ fluxes: (a) N₂O fluxes at −30 hPa, (b) N₂O fluxes at −100 hPa, (c) N₂ fluxes at −30 hPa, (d) N₂ fluxes at −100 hPa, (e) CO₂ fluxes at −30 hPa, and (f) CO₂ fluxes at −100 hPa. The N₂ fluxes used for comparison with model results were estimated from the observed N₂O fluxes and ¹⁵N-labeled ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}/{\mathrm{N}}_{\mathrm{2}}$ ratios. Error bars represent standard deviations (n=3).

Mineral N and DOC concentrations in the soil profile were also compared with model simulations (Figs. 3 and 4). The dissolved NH${}_{\mathrm{4}}^{+}$ in Figs. 3a and 4a, used to compare with model results, was estimated according to a Freundlich adsorption isotherm as described in Sect. 2.4.1. In the manure zone, the estimated concentrations of dissolved NH${}_{\mathrm{4}}^{+}$ by day 1 had dropped from the initial value of 23.1 to 4.1 mmol L⁻¹ in the −30 hPa treatment and from 24.1 to 4.8 mmol L⁻¹ in the −100 hPa treatment. The total NH${}_{\mathrm{4}}^{+}$ content in the manure zone was calculated to be 435 mg N kg⁻¹ initially if distributed within the 8 mm layer, which was approximately 4 times the content actually observed on day 1 of 112 and 123 mg N kg⁻¹ in the −30 and −100 hPa treatments, respectively (Figs. S7a and S8a). In the model, neither the dissolved nor the total NH${}_{\mathrm{4}}^{+}$ content decreased as much as that observed in the experiment by day 1, and it took 4 d before the majority of NH${}_{\mathrm{4}}^{+}$ was reduced to the level observed by day 1 in the experiment. By day 14, the dissolved NH${}_{\mathrm{4}}^{+}$ content in the manure hotspot was reduced to the level in the bulk soil both in the experiment and in the model.

Figure 3 Model results (b, d, f, h) of soil variables (DOC and dissolved N species) compared with direct measurements or estimations (a, c, e, g) in the −30 hPa treatment: (a) estimated NH${}_{\mathrm{4}}^{+}$ concentrations from the experiment, (b) modeled NH${}_{\mathrm{4}}^{+}$ concentrations, (c) measured NO${}_{\mathrm{3}}^{-}$ concentrations from the experiment, (d) modeled NO${}_{\mathrm{3}}^{-}$ concentrations, (e) estimated DOC concentrations from the experiment, (f) modeled DOC concentrations, (g) measured NO${}_{\mathrm{2}}^{-}$ concentrations from the experiment, and (h) modeled NO${}_{\mathrm{2}}^{-}$ concentrations. In the legend, IC indicates the initial condition applied in the model to represent the lab experiment with slurry application. Error bars represent standard deviations (n=2).

Figure 4 Model results (b, d, f, h) of soil variables (DOC and dissolved N species) compared with direct measurements or estimations (a, c, e, g) in the −100 hPa treatment: (a) estimated NH${}_{\mathrm{4}}^{+}$ concentrations from the experiment, (b) modeled NH${}_{\mathrm{4}}^{+}$ concentrations, (c) measured NO${}_{\mathrm{3}}^{-}$ concentrations from the experiment, (d) modeled NO${}_{\mathrm{3}}^{-}$ concentrations, (e) estimated DOC concentrations from the experiment, (f) modeled DOC concentrations, (g) measured NO${}_{\mathrm{2}}^{-}$ concentrations from the experiment, and (h) modeled NO${}_{\mathrm{2}}^{-}$ concentrations. In the legend, IC indicates the initial condition applied in the model to represent the lab experiment with slurry application. Error bars represent standard deviations (n=2).

In the experiment the NO${}_{\mathrm{3}}^{-}$ content in the manure zone was negligible 1 d after manure addition in both treatments. In the −30 hPa treatment, NO${}_{\mathrm{3}}^{-}$ content then increased and approached the background level in the unamended control soil after 4 weeks. In the simulation, NO${}_{\mathrm{3}}^{-}$ content was comparable with the measurements on day 1 and day 14 but showed a net removal of NO${}_{\mathrm{3}}^{-}$ in the vicinity of the manure zone on day 28, which was not observed in the experiment. In the −100 hPa treatment of the experiment, NO${}_{\mathrm{3}}^{-}$ content showed an overall increase within the soil profile after 2 and 4 weeks compared to day 1 (Fig. 4c). The NO${}_{\mathrm{3}}^{-}$ content on day 28 in the soil outside the manure zone was relatively higher than that on day 14, whereas the NO${}_{\mathrm{3}}^{-}$ content at the center of the manure zone was lower by day 28 compared to day 14, and this trend was well captured by the model (Fig. 4d).

The DOC content was estimated from measured values of LOI and manure VSs for comparison with model results. The estimated DOC profile on day 1 showed total amounts of 20.9 and 23.9 mmol in the −30 and −100 hPa treatments, respectively. The modeled DOC by day 1 showed comparable total amounts of 27.4 and 26.8 mmol in the −30 and −100 hPa treatments, respectively, but with a wider distribution in the soil than that estimated in the experiment (Figs. 3e, f and 4e, f). According to the model simulation, the DOC in the manure zone was mostly removed within 2 weeks in both treatments, but no further organic matter analyses were done in the experiment to confirm this.

For both treatments the modeled NO${}_{\mathrm{2}}^{-}$ content showed a similar temporal trend to the experimental data. The total NO${}_{\mathrm{2}}^{-}$ content in the soil was higher on day 1 than after 14 and 28 d of incubation (Figs. 3g, h and 4g, h). The simulation values were generally larger than those recorded, but considering the analytical uncertainty about the low NO${}_{\mathrm{2}}^{-}$ levels in the experiment, this discrepancy was not unacceptable. In contrast to the bell-shaped observed NO${}_{\mathrm{2}}^{-}$ profile on day 1, the simulated NO${}_{\mathrm{2}}^{-}$ profile showed a trough at the center of the manure-saturated zone where NO${}_{\mathrm{2}}^{-}$ concentration developed more in the following days.

3.2 Spatiotemporal features of simulated N₂O and other variables

Figure 5 shows the modeled process rates related to N₂O transformations in the −30 hPa treatment, including N₂O production via nitrification, nitrifier denitrification, and denitrification; N₂O reduction via denitrification; and rates of N₂O change via diffusion. According to the model, N₂O production mainly occurred at the depth interval 0.04–0.06 m encompassing the manure zone from 0.046 to 0.054 m (Fig. 5a–d). At the center of the manure zone, around the 2 mm saturated area, denitrification was the major pathway of N₂O production, with the highest rate of 6 mmol N₂O L⁻¹ soil d⁻¹ on day 3 (Fig. 5d). Nitrification contributed to N₂O production in the zone outside the 2 mm saturated area and peaked around the soil–manure interface (Fig. 5b). N₂O production via nitrification mainly occurred during the first week, and the reaction rates up to 0.07 mmol N₂O L⁻¹ soil d⁻¹ were markedly lower than the simultaneous denitrification rates according to the model. The modeled nitrifier denitrification took place close to the saturated zone where the reaction rates were much higher than in other layers (Fig. 5c). Similar to nitrification and denitrification, N₂O production from nitrifier denitrification mainly occurred during the first week, and the highest rate was approx. 0.1 mmol N₂O L⁻¹ soil d⁻¹ on day 3.

Figure 5 Simulated rates of N₂O production, consumption, and transport in the −30 hPa treatment: (a) net N₂O production rate, (b) N₂O production rate by nitrification, (c) N₂O production rate by nitrifier denitrification, (d) N₂O production rate by denitrification, (e) N₂O consumption rate by denitrification, (f) N₂O diffusive rate, (g) the net rate of N₂O changes by reactions and transport, and (h) N₂O diffusive flux, where the negative sign represents the downward movement towards the lower soil–air interface (z=0.1 m) and the positive sign the flow towards the upper soil–air interface (z=0 m).

Microbial reduction in N₂O to N₂ showed a similar spatiotemporal pattern to N₂O production via denitrification (Fig. 5e). In accordance with this, the net reaction, i.e., the sum of production and consumption rates, peaked at the center of the manure zone and had a moderate intensity in the soil volumes mainly associated with nitrification (Fig. 5a). Negative rates around the central manure-saturated zone were paralleled by N₂ production (Fig. 5a).

The rates in Fig. 5f, i.e., the first term on the right-hand side of Eq. (1), indicate that N₂O concentration change is caused by gas diffusion. The negative values in the soil profile represent the decrease in N₂O concentrations owing to efflux to the surrounding soil which mainly occurred in the manure zone. Positive values represent a gain in N₂O concentrations due to gas diffusion and mainly occurred in the bulk soil where the generated N₂O accumulated before escaping to the atmosphere. The net change in N₂O concentrations in the soil showed an overall increase on day 1 and then a decrease on day 3 that was followed by minor changes on the following sampling days according to the simulation (Fig. 5g). The direction and magnitude of the internal gas transport of N₂O are illustrated in Fig. 5h.

In the −100 hPa treatment, the intensities of N₂O production from the three processes were different from the −30 hPa treatment, whereas the spatial distributions of the individual processes were similar to the patterns observed in the −30 hPa treatment (Fig. S10a–d). Nitrification rates were generally higher than those of the −30 hPa treatment, especially at the center of the manure zone within the first 3 d (Figs. S10b and 5b). Specifically, in the −30 hPa treatment the nitrification rate in the central manure-saturated zone was approx. 0.02 mmol N₂O L⁻¹ soil d⁻¹ between day 1 and day 3, which was ca. 50 % lower than the rates at the manure–soil interface (Fig. 5b). In the −100 hPa treatment, the corresponding rate at the manure center increased from 0.04 to 0.06 mmol N₂O L⁻¹ soil d⁻¹ between day 1 and day 3 in the model simulation, and a higher rate was also found at the manure–soil interface (Fig. S10b). Denitrification and nitrifier denitrification were less intensive in the −100 hPa treatment compared to the −30 hPa treatment. N₂O produced from nitrifier denitrification was up to 0.04 mmol N₂O L⁻¹ soil d⁻¹ on day 3 in the −100 hPa treatment compared to a rate of 0.1 mmol N₂O L⁻¹ soil d⁻¹ in the −30 hPa treatment (Figs. S10c and 5c). Denitrification-derived N₂O showed the highest rate with up to 0.8 mmol N₂O L⁻¹ soil d⁻¹ on day 1, markedly lower than the peak value of ca. 6 mmol N₂O L⁻¹ soil d⁻¹ on day 3 in the −30 hPa treatment (Figs. S10d and 5d). The net changes in N₂O concentrations showed an increase during the first 3 d and slight changes afterwards (Fig. S10g).

The fate of NO${}_{\mathrm{3}}^{-}$ in the model was determined by three processes: production via nitrification, consumption via denitrification, and diffusive exchange. According to the model, NO${}_{\mathrm{3}}^{-}$ production mainly occurred outside the saturated zone and peaked around the borders of the manure-saturated zone, where the rate was up to 1.4 mmol NO${}_{\mathrm{3}}^{-}$ L⁻¹ soil d⁻¹ by day 1, and the peak values continued to increase to reach 2.2 mmol NO${}_{\mathrm{3}}^{-}$ L⁻¹ soil d⁻¹ by day 3 (Fig. 6a). The increase in NO${}_{\mathrm{3}}^{-}$ production rate continued until day 7, followed by a decrease to the end of the 4-week simulation. The consumption of NO${}_{\mathrm{3}}^{-}$ via denitrification dominated at the center of the manure zone (Fig. 6b) where the soil was saturated and the O₂ availability was low (Figs. S4 and S9). The NO${}_{\mathrm{3}}^{-}$ consumption rate greatly increased during the first week with a maximum value of around 7 mmol NO${}_{\mathrm{3}}^{-}$ L⁻¹ soil d⁻¹ on day 3 (Fig. 6b), greater than the simultaneous production rate. A net consumption rate within the manure-saturated zone was therefore simulated (Fig. 6c) based on these two microbial pathways. Besides biochemical reactions, the supply of NO${}_{\mathrm{3}}^{-}$ used for denitrification also came from the process of NO${}_{\mathrm{3}}^{-}$ diffusion, as depicted in Fig. 6d. Driven by the concentration gradient between the manure zone and the bulk soil, nitrification-derived NO${}_{\mathrm{3}}^{-}$ migrated to the center of the manure hotspot where NO${}_{\mathrm{3}}^{-}$ was effectively consumed, as indicated in Fig. 6f. Between day 1 and day 3, the NO${}_{\mathrm{3}}^{-}$ diffusion rate within the saturated zone doubled from ca. 3 mmol NO${}_{\mathrm{3}}^{-}$ L⁻¹ soil d⁻¹ to the maximum value of ca. 6 mmol NO${}_{\mathrm{3}}^{-}$ L⁻¹ soil d⁻¹, compensating for the relatively lower nitrification rate to sustain denitrification (Fig. 6d).

Figure 6 Simulated rates of NO${}_{\mathrm{3}}^{-}$ production, consumption, and transport in the −30 hPa treatment: (a) NO${}_{\mathrm{3}}^{-}$ production rate by nitrification, (b) NO${}_{\mathrm{3}}^{-}$ production rate by denitrification, (c) NO${}_{\mathrm{3}}^{-}$ net reaction rate, (d) NO${}_{\mathrm{3}}^{-}$ diffusive rate, (e) the net rate of NO${}_{\mathrm{3}}^{-}$ changes by reactions and transport, and (f) NO${}_{\mathrm{3}}^{-}$ diffusive flux. In panel (f), the negative sign represents the downward movement towards the lower soil–air interface (z=0.1 m) and the positive sign the flow towards the upper soil–air interface (z=0 m).

In the −100 hPa treatment, the spatial patterns of NO${}_{\mathrm{3}}^{-}$ transformations were different from the −30 hPa treatment, and here high rates of NO${}_{\mathrm{3}}^{-}$ production also occurred at the center of the manure zone (Fig. S11) where the water content was relatively higher than in the surrounding soil but lower than in the −30 hPa treatment. The diffusion rate ranged from 0.8 to 1.6 mmol NO${}_{\mathrm{3}}^{-}$ L⁻¹ soil d⁻¹ at the center of the manure zone (Fig. S11d), and the production rate ranged from 1.8 to 2.8 mmol NO${}_{\mathrm{3}}^{-}$ L⁻¹ soil d⁻¹ (Fig. S11a). Figure S11b shows that NO${}_{\mathrm{3}}^{-}$ reduction during denitrification mainly occurred on day 1 at around 1.8 mmol NO${}_{\mathrm{3}}^{-}$ L⁻¹ soil d⁻¹, which declined greatly to 0.4 mmol NO${}_{\mathrm{3}}^{-}$ L⁻¹ soil d⁻¹ during the following 2 d and thereafter continued to show a declining trend. A net positive change in NO${}_{\mathrm{3}}^{-}$ concentration was simulated during the first 3 d, and then a net NO${}_{\mathrm{3}}^{-}$ depletion was predicted in the manure zone after 1 week (Fig. S11e).

The biomass of nitrifiers (i.e., AOB and NOB) and denitrifiers increased rapidly during the first week in the simulations (Figs. S5 and S6) and dominated in spatially distinct niches. Thus, the growth of nitrifying bacteria mainly occurred inside and around the manure-affected area, with a visible increase in the soil at 6–8 mm from the initial borders of manure application. The biomasses of AOB and NOB peaked by day 7 close to the soil–manure interface and were by then, respectively, up to 4 times as high and twice as high as the background level in both treatments. Microbial decay then dominated with the result that the nitrifier biomass decreased but still remained above the background level by day 28. The denitrifier biomass was initially unevenly distributed within the soil where the 2 mm center of the manure zone had a denitrifier biomass 10 times that of the bulk soil to account for denitrifiers existing in the liquid manure, as described in Sect. 2.4.1. The growth of denitrifiers mainly occurred inside the manure zone and peaked at up to 30 times the background level in bulk soil by day 10 in the −30 hPa treatment, which was followed by a temporally stable phase until day 21 and then a declining trend until the end of the experiment.

In the −30 hPa treatment, the O₂ consumption rate was well aligned with the CO₂ production rate through aerobic respiration (Figs. S14b and S16b). By day 1 and day 3, O₂ consumption by aerobic respiration ($S_{{\mathrm{O}}_{\mathrm{2}},\mathrm{r}}$) showed peak values of −32 and −83 mmol O₂ L⁻¹ soil d⁻¹, respectively, around the soil–manure interface, and simultaneously the O₂ consumption via nitrification ($S_{{\mathrm{O}}_{\mathrm{2}},\mathrm{n}}$) had peak values of −4.5 and −7.5 mmol O₂ L⁻¹ soil d⁻¹ in the same position. At the center of the manure zone, CO₂ was mainly produced by respiration of denitrifiers during the simulated 4-week period (Fig. S14c), with up to 12 mmol CO₂ L⁻¹ soil d⁻¹ on day 1, in contrast to the aerobic respiration that dominated in most of the core volume. In the −100 hPa treatment, oxygen consumption was mainly owing to aerobic respiration by day 1, where the peak values of $S_{{\mathrm{O}}_{\mathrm{2}},\mathrm{r}}$ and $S_{{\mathrm{O}}_{\mathrm{2}},\mathrm{n}}$ were −73 and −5 mmol O₂ L⁻¹ soil d⁻¹. On day 3 nitrifier denitrification dominated the use of O₂ within the manure-application zone (peak: −8.6 mmol O₂ L⁻¹ soil d⁻¹), whereas aerobic respiration has greatly declined and was only elevated at the center (peak: −7.4 mmol O₂ L⁻¹ soil d⁻¹) (Figs. S15 and S17). CO₂ produced by denitrification peaked at the center of the manure zone but was overall less than aerobic respiration during the 4 weeks in the −100 hPa treatment (Fig. S15).

The oxygen diffusion rates closely followed the consumption rates (Figs. S16a, d and S17a, d). This indicates that the overall air-filled porosity within the 10 cm soil core, including a small fraction of air-filled porosity (i.e., 0.01) in the saturated zone, was able to ensure adequate soil aeration at the oxygen consumption rates occurring in the experiment. The simulated oxygen content in the pore air began to increase within 1 d (Figs. S5 and S6). However, the O₂ availability to enzyme reaction sites in the model depends not only on the oxygen content in soil pore air but also on the soil tortuosity in relation to the soil water content. As shown in Fig. S9, the simulated O₂ availability showed dramatic heterogeneity across the soil profile, with a steep downward slope at the center of the manure zone, but no evident temporal changes were simulated after day 1 in the two treatments. The −100 hPa treatment in general showed a higher level of O₂ availability than the −30 hPa treatment, but the O₂ available in the near-saturated center (i.e., 0.049–0.051 m) was still lower than that in the −30 hPa bulk soil and thus facilitated local denitrification.

3.3 Scenario tests: the importance of solute diffusion in the model

Figure 7 shows the results of four scenarios to investigate the importance of solute diffusion for gas emissions in the baseline model scenarios and, presumably, measurement results presented above. Switching off the diffusion of all solutes greatly reduced N₂O. N₂ fluxes during the 4-week simulation of the −30 hPa treatment were compared to the baseline simulation and the observations, and the emission peaks during the first 3 d disappeared (Fig. 7a and c). By allowing the diffusion of NO${}_{\mathrm{3}}^{-}$ but not DOC, NH${}_{\mathrm{4}}^{+}$, and NO${}_{\mathrm{2}}^{-}$ (Scenario 2), N₂O and N₂ fluxes were enhanced to some degree but still could not explain the temporal patterns in the observation. Allowing NH${}_{\mathrm{4}}^{+}$ diffusion only (Scenario 3) increased the simulated N₂O and N₂ fluxes at the start of the simulation compared to Scenario 1, although not as much as the change caused by allowing NO${}_{\mathrm{3}}^{-}$ diffusion only. The simulated N gas fluxes dropped to negligible levels after 1 week in Scenario 3. In Scenario 4, allowing the diffusion of three N solutes but not DOC, the simulated N₂O and N₂ fluxes were greatly increased compared to Scenarios 1–3, and the temporal pattern, characterized by an early flux peak and a long tail, was similar to the baseline. However, the timing of peaks came later than in the measurement, and the simulated peak was even higher when preventing DOC flux between manure layer and soil. Simulated CO₂ emission showed a peak around day 2 in the baseline scenario, which was delayed by approx. 1 d and declined from ca. 5×10⁵ to ca. 2×10⁵ µg C m⁻² h⁻¹ in the four scenarios (Fig. 7e). Not much difference was found in the simulated CO₂ emissions between the four scenarios.

Figure 7 Comparison of different scenarios with respect to N₂O, N₂, and CO₂ fluxes: (a) N₂O fluxes at −30 hPa, (b) N₂O fluxes at −100 hPa, (c) N₂ fluxes at −30 hPa, (d) N₂ fluxes at −100 hPa, (e) CO₂ fluxes at −30 hPa, and (f) CO₂ fluxes at −100 hPa. In each panel, there are measured data, a baseline simulation where the diffusion of all solutes is included, and four scenarios (1–4), where Scenario 1 does not allow any solute diffusion; Scenario 2 allows only NO${}_{\mathrm{3}}^{-}$ diffusion; Scenario 3 allows only NH${}_{\mathrm{4}}^{+}$ diffusion; and Scenario 4 allows diffusion of NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$, and NO${}_{\mathrm{2}}^{-}$ but not DOC. Error bars represent standard deviations (n=3).

In the −100 hPa treatment, Scenario 1 and Scenario 2 showed similar results in that the early N₂O and N₂ fluxes within the first 3 d were lower than that in the baseline scenario but still higher than the observations (Fig. 7b and d). The differences between the baseline and the two scenarios were much smaller than that in the −30 hPa treatment. Whereas the inclusion of NO${}_{\mathrm{3}}^{-}$ diffusion (Scenario 2) did not show any stimulation of N₂O emissions compared to Scenario 1, including NH${}_{\mathrm{4}}^{+}$ diffusion resulted in a marked increase in simulated gas fluxes at the early stage which was identical with the result of Scenario 4. The measured gas fluxes were overestimated in Scenario 3 and Scenario 4, similar to the tendency in the −30 hPa treatment. The simulated CO₂ peak flux was lower than the baseline simulation and was delayed by approx. 1 d (Fig. 7f).

4 Discussion

4.1 Model performance

The model generally reflected the magnitudes of individual N₂O and N₂ emissions measured in the experiment where two moisture levels were tested, with the larger model errors in simulating the −100 hPa treatment as shown in Table S7. In the −30 hPa treatment, the predicted cumulative N₂O and N₂ fluxes were close to the measured data with the relative mean error (ME) of 0.19 and −0.29, respectively, while in the −100 hPa treatment, the predicted cumulative fluxes were approx. double or 4 times those of the observed data (relative ME: 0.93–3.54). This is attributed to the overestimation of gas fluxes in the early stage of the −100 hPa treatment, and the total errors indicated by rRMSE are more sensitive to such extreme daily fluxes. Previous studies have also reported high variability in model performance of predicting daily N₂O fluxes (e.g., Grosz et al., 2021; Zimmermann et al., 2018). Since the model captured well the higher fluxes of N₂O and N₂ in the −30 hPa treatment and the temporal trend of observed fluxes with R² values ranging from 0.43 to 0.71, model predictions were considered acceptable. Also, the main purpose of this study was not to accurately reproduce emissions but to investigate C and N transformations via major microbial pathways in a soil environment with liquid manure representing a type of organic hotspot that is characteristic of intensive agriculture. The different model-data errors under contrasting moisture levels could be related to the use of relationships between soil moisture and diffusivity factors to calculate the availability of substrates at the sites of enzymatic reactions rather than empirical moisture response functions, as process-based models often do. Since the distribution of soil water content in the soil cores remained almost constant during the 4-week incubation (Fig. S4), temporal variations in soil moisture responses and gas diffusivity associated with water movement were less important than solute diffusion in the experimental setup used in this study.

The simulated profiles of solutes, including mineral N and DOC, showed generally the same spatiotemporal patterns as the measured or estimated data from the experiment, albeit with some discrepancies. The overestimation of NH${}_{\mathrm{4}}^{+}$ on day 1 (Figs. 3a, b and 4a, b) may have several reasons. First, there was some uncertainty about the NH${}_{\mathrm{4}}^{+}$ initial condition set in the model due to a potential for ammonia volatilization from the cattle slurry while infiltrating into the soil during sample preparation. Next, we assumed that a few hours after manure application the NH${}_{\mathrm{4}}^{+}$ content was uniformly distributed within an 8 mm depth interval, but the liquid manure could have spread to a wider range and hence become more diluted. These two mechanisms would lead to overestimation of the initial NH${}_{\mathrm{4}}^{+}$ concentrations in the modeled manure zone compared to the actual level on day 0 in the experiment. Uncertainty about the NH${}_{\mathrm{4}}^{+}$ adsorption isotherm (Hunt and Adamsen, 1985; Sieczka and Koda, 2016) could be another reason for the high model estimates of NH${}_{\mathrm{4}}^{+}$ in the manure zone if the modeled NH${}_{\mathrm{4}}^{+}$ mobility was lower than the actual mobility in the loamy sandy soil. Ammonium uptake by bacteria as part of osmoregulation could also have contributed to the rapid disappearance of NH${}_{\mathrm{4}}^{+}$ in the experiment, since the osmotic pressure of solutes in livestock slurry can exceed −1000 kPa (Petersen and Andersen, 1996), and hyperosmotic conditions induce bacteria to accumulate amino acids (Li et al., 2011). Also, the immobilization of NH${}_{\mathrm{4}}^{+}$ would occur in connection with the bacterial growth on the basis of DOC (Recous et al., 1990). The processing of NH${}_{\mathrm{4}}^{+}$ by AOB was also included in the model and partly accounted for this process, while the related parameter, yield coefficient, was fixed by taking values from the literature where uncertainty existed.

The spatiotemporal pattern of modeled NO${}_{\mathrm{3}}^{-}$ aligned well with measurements in the −100 hPa treatment, while in the −30 hPa treatment a depletion of NO${}_{\mathrm{3}}^{-}$ was simulated around the manure zone by day 28, which was absent from the NO${}_{\mathrm{3}}^{-}$ profile observed. The rebound of NO${}_{\mathrm{3}}^{-}$ in the experiment could be related to the mineralization of labile organic N in the manure followed by nitrification, a mechanism that is currently not included in the model and a potential source of model structure uncertainty. Møller et al. (2004) reported an average protein content of 150 g kg⁻¹ of VSs in cattle slurry, part of which could have been degraded within the first days of incubation. The modeled denitrifier biomass in the −30 hPa treatment remained much higher (>23 times) in the manure zone compared to the background level on day 28, whereas the biomass of AOB and NOB had declined to less than twice the background level (Fig. S5). In an experimental study by Petersen et al. (1992), the potential for denitrification in the manure-saturated zone decreased in the late stage of the incubation, which was explained by the decay of microorganisms or inhibition of enzymatic activities caused by NO${}_{\mathrm{2}}^{-}$ toxicity, specifically the unionized nitrous acid species (Abeling and Seyfried, 1992). Microbial decay in the model followed first-order kinetics with a constant decay rate, and it may be possible in future work to improve relevant parameters and consider the response to substrates on microbial decay and maintenance based on previous studies (Ni et al., 2011).

The spatial stratification of nitrifier and denitrifier growth in the model (Figs. S5 and S6) was consistent with experimental results (e.g., Nielsen et al., 1996; Petersen et al., 1992) showing a very high potential for nitrification developing within a few millimeters from the manure–soil interface and the potential for denitrification growing within the manure-saturated zone. In both treatments, after significant growth of microbial populations, the reduction in biomass of aerobic heterotrophs, AOB, NOB, and denitrifiers indicated a predominance of decay over growth which had been gradually limited by the depletion of substrates, i.e., DOC for heterotrophs and denitrifiers, ammonium for AOB, and nitrite for NOB, as the simulated nitrogen removal from the soil continued. The spatial distributions of aerobic heterotrophs and nitrifying bacteria were similar to and distinctly different from that of denitrifying bacteria.

Variation in the availability of O₂ at enzymatic sites was the main factor accounting for the stratification of microbial communities in the model, where the growth of denitrifying bacteria was promoted at the center of the manure zone with O₂ limitation, while aerobic bacteria preferentially developed outside the saturated or near-saturated zone with better access to O₂. The simulated O₂ concentration in soil pore air indicated that O₂ re-entered the center of the manure zone within the first 24 h (Figs. S5 and S6). If the anoxic volume had been larger, as observed in experiments with centimeter-scale hotspots (Markfoged et al., 2011; Petersen et al., 1996), the O₂ depletion zone would probably have been larger and temporally more stable. Thus, micro-sensor measurements showed an oxygen penetration into a manure hotspot of less than 1 mm during the first week (Petersen et al., 1996). The simulated rapid recovery of O₂ in the air-filled pores in the manure zone predicted by the model was related to the small extent of the anoxic zone that was predicted for the experimental system, which was as narrow as 2 mm. However, the limited gas diffusivity within the manure zone still led to a low level of O₂ available for enzymatic reactions throughout the 4-week simulation (Fig. S9).

4.2 Effects of solute diffusion

Switching off solute diffusion between layers greatly reduced the modeled N₂O emissions in the −30 hPa treatment (Fig. 7), whereas at −100 hPa soil water potential the effect was much less and indicated that solute diffusion between manure and soil layers was already low. The contrasting effects that we found in the scenario tests with the complete or partial elimination of solute diffusion were considered to reflect the importance of different mechanisms contributing to N₂O emissions at the two moisture levels.

Dissimilatory NO${}_{\mathrm{3}}^{-}$ reduction by denitrifying bacteria initiates the chain of reactions resulting in a stepwise reduction in nitrogenous oxides to N₂ and emission to the ambient air. In the −30 hPa treatment, a hotspot of NO${}_{\mathrm{3}}^{-}$ consumption was present at the center of the manure zone where, by day 1, the diffusive supply was ca. 103 % of the NO${}_{\mathrm{3}}^{-}$ demand by denitrification and by day 3 ca. 99 % (Fig. S18). After day 7, NO${}_{\mathrm{3}}^{-}$ generated from nitrification accounted for a greater proportion of the demand for denitrification and was 12 %–17 % in the second week. The NO${}_{\mathrm{3}}^{-}$ diffusing from the soil–manure interface towards the anoxic center was due to nitrification outside the hotspot, as well as transport of soil NO${}_{\mathrm{3}}^{-}$. Eliminating the diffusion of NO${}_{\mathrm{3}}^{-}$, the initially low NO${}_{\mathrm{3}}^{-}$ at the anoxic center could not be replenished by either source, thereby hindering the potential for denitrification in the simulation. The lack of NO${}_{\mathrm{3}}^{-}$ transport could account for the dramatic decrease in N₂O flux in the Scenario 1 test for the −30 hPa treatment. By day 1, the NO${}_{\mathrm{3}}^{-}$ produced from nitrification, which could be estimated from Fig. 6a and soil moisture, was still lower than the initial NO${}_{\mathrm{3}}^{-}$ content in the bulk soil, which was 4.5 mmol L⁻¹ in the −30 hPa soil and 6.2 mmol L⁻¹ in the −100 hPa soil (Figs. S5 and S6). Nitrification-derived NO${}_{\mathrm{3}}^{-}$ gradually became comparable to, and even exceeded, the soil NO${}_{\mathrm{3}}^{-}$ level after day 3. However, the transport of NO${}_{\mathrm{3}}^{-}$ from the bulk soil, initially present at z>0.046 m and z>0.054 m, to the anoxic zone (0.049–0.051 m) became less efficient than the nitrification-derived NO${}_{\mathrm{3}}^{-}$ which was produced immediately outside the anoxic zone, thus increasing the relative importance of nitrification-derived NO${}_{\mathrm{3}}^{-}$ over time, as shown previously in controlled experiments with ¹⁵N enrichment of the external NO${}_{\mathrm{3}}^{-}$ supply (Nielsen et al., 1996). Allowing NO${}_{\mathrm{3}}^{-}$ diffusion only (Scenario 2) did not greatly improve the N₂O emission compared to Scenario 1, which reflects the limited role of NO${}_{\mathrm{3}}^{-}$ supply from the bulk soil owing to the increasing transport distance and absence of nitrification without access to NH${}_{\mathrm{4}}^{+}$. However, allowing for NO${}_{\mathrm{3}}^{-}$ diffusion in the −30 hPa soil was still better than allowing NH${}_{\mathrm{4}}^{+}$ diffusion only in this treatment (Scenario 3), and hence the co-occurrence of nitrification and denitrification within the manure zone was insufficient to account for the observed N₂O emissions. These limitations could be alleviated by allowing the diffusion of all nitrogenous substrates from the simulation (Scenario 4), allowing for the coupling between nitrification outside the anoxic zone and denitrification inside this zone.

In the −100 hPa treatment, in contrast, the spatial stratification of nitrification and denitrification was less evident because the production of NO${}_{\mathrm{3}}^{-}$ within the denitrification zone was only moderately lower than the production around the interface and could cover the demand for denitrification in this scenario (Fig. S11). Specifically, the nitrification rate and diffusion rate within the 0.049–0.051 m interval were, respectively, 103 % and 89 % of the denitrification rate on day 1, and the nitrification rate continued to increase by day 3 (Fig. S19), showing that the NO${}_{\mathrm{3}}^{-}$ consumption within the denitrification center in the −100 hPa treatment was much lower than the NO${}_{\mathrm{3}}^{-}$ supply, with nitrification rather than solute diffusion as the dominating source. Compared with the −30 hPa treatment, the larger air-filled porosity (Fig. S3) and higher availability of O₂ (Fig. S9) in the −100 hPa treatment thus supported a higher nitrification activity which was sufficient to sustain the denitrification activity at the center of the manure zone. The limited solute diffusion in the −100 hPa treatment was due to the much lower effective diffusion coefficients in the aqueous phase, and with more extensive nitrification, gradients of NO${}_{\mathrm{3}}^{-}$ concentrations between the soil and the manure zone could be smaller relative to that in the −30 hPa treatment. The simulated N₂O emission was not influenced by cutting off the solute diffusion as much as that of the −30 hPa treatment, indicating that tightly coupled nitrification and denitrification within the manure zone, rather than solute diffusion from surrounding soil, sustained denitrification. This is in accordance with the results of Manzoni and Katul (2014), who found that hydrological connectivity may continue to exist at microscale even when macroscale solute diffusivity is interrupted. Hence, the model predicted that when manure is applied to drier soil, the redistribution of water may support an even closer coupling between nitrification and denitrification than is the case in wetter soil with much less redistribution of liquid manure (Petersen et al., 2003).

Some experimental studies (Meyer et al., 2002; Nielsen et al., 1996; Petersen et al., 1991, 1992) have reported that in manure-amended soil there may be a distance of only a few millimeters between active nitrification and denitrification zones, facilitating coupled nitrification–denitrification. These observations are consistent with our model predictions for the −30 hPa treatment, where a very pronounced stratification of microbial populations was simulated. We also found support for an even closer, though short-lived, coupling of nitrification and denitrification in the simulation of the −100 hPa treatment. The manure hotspots investigated here were relatively small, with only a 2 mm layer dominated by manure, and field application of liquid cattle manure is likely to form, depending on application method, centimeter-scale hotspots with a greater potential for developing coupled nitrification–denitrification. Based on the discussion above, we suggest that solute diffusion should be considered a key control of denitrification and N₂O emissions in soil–manure systems. The extent will depend on factors such as soil water content and texture controlling water film continuity, as well as manure application rate and application method determining soil–manure contact and stability of anoxic volumes that can support denitrification.

4.3 Model uncertainties and future applications

The model presented was able to qualitatively reflect gaseous emissions and changes in carbon and nitrogen pools in and around an organic hotspot in soil, but several aspects of model uncertainty remain. First, the initial condition used in the model was a simplification of actual conditions due to the lack of detailed soil analysis within the first 24 h of incubation. We considered the stabilized slurry redistribution as a “box”; i.e., liquid manure was redistributed by piston flow away from the initial placement of the added manure. However, in reality the joint effects of dispersion and diffusion could result in a more continuous distribution of dissolved slurry components over a wider zone. Also, some reactive C and N components contained in the manure may have been lost between the time of manure sampling for chemical analyses and the initiation of incubations. These two factors could result in the estimated initial concentrations of manure components in the assumed manure zone being overestimated. We suggest that if incubation experiments are designed to support modeling, more frequent sampling early in the experiment can provide an opportunity to better characterize the behavior of the system.

In the model, we assumed the gases generated from biochemical reactions existed in air-filled pores instead of modeling the water–air exchange. Therefore, the current model is suitable for simulating processes in unsaturated rather than saturated soil where N₂O and CO₂ can only exist in dissolved form. However, when saturated volumes only account for a small fraction of the total samples, small adaptations can be made in model application. In this study, we included a small fraction of air-filled porosity (i.e., 0.01) in the 2 mm saturated zone when modeling the −30 hPa treatment to ensure the validity of Eq. (1) over the entire model domain. The value of the air-filled porosity did influence the absolute values of simulated gas fluxes. However, when the value of the air-filled porosity was kept within a reasonable range, specifically well below the air-filled porosity (i.e., 0.042) at the center of the −100 hPa treatment samples, the dominating processes associated with the emissions were not influenced much. Thus, when keeping the zero air-filled porosity in the saturated zone or changing the added air-filled porosity to be 0.001, the modeled N₂O peak, compared to the baseline simulation and observations, showed an increase in both cases in the −30 hPa treatment, as anoxia dominated more at the center of soil cores, whereas the qualitative results regarding solute diffusion and the spatiotemporal patterns of reaction rates remained the same (Figs. S22 and S23). Based on the above two air-filled porosity settings in the −30 hPa treatment, we also attempted to re-optimize the model parameters in order to specifically reduce the model-data errors, and the main conclusions we obtained were not changed (data not shown).

In describing the microbial activities in the model, the possible substrate inhibition of nitrifiers and denitrifiers was neglected. For example, denitrification may be inhibited by high concentrations of NO${}_{\mathrm{3}}^{-}$ (Francis and Mankin, 1977), which has been attributed to toxicity in situations where NO${}_{\mathrm{2}}^{-}$ accumulates (Abeling and Seyfried, 1992; Glass et al., 1997). The mineralization of organic N, including the transfer of microbial organic N to mineral pools via microbial death, was also not included in the model at this time and is a topic for future research. However, in the manure hotspot environment the N available for nitrifiers is dominated by ammonium derived from urea in excretal returns or organic matter decomposed during storage (Jost et al., 2013; Webb et al., 2013). Following the study by Chen et al. (2019), only AOB and NOB were capable of N incorporation in the model, and heterotrophic growth was linked to C incorporation. It should be noted that the focus of the study was on the coupled nitrification–denitrification and the resulting spatiotemporal dynamics of multiple C and N species, and the consideration of additional microbial processes, e.g., detailed mineralization and immobilization, would further complicate the system.

The model was able to depict the spatiotemporal patterns of the main N₂O forming processes in soil including nitrification, nitrifier denitrification, and denitrification. However, uncertainty exists in the relative contributions of these processes, especially the nitrifier denitrification pathway, as relevant parameters were not well constrained by the available measurements, and the simulated gas fluxes were less sensitive to the parameters regarding nitrifier denitrification than to the parameters regarding the other two processes (Fig. S24). Detailed quantification of N₂O sources from different processes was beyond the scope of our study, but isotopic signatures of N₂O are a promising tool in this respect, as shown in the modeling study by Chang et al. (2022).

For a complex biochemical model such as the one presented, the detailed parameter design gave physical meaning to each parameter, allowing great flexibility in describing the processes, but there is a potential for parameter equifinality. When analyzing parameter correlations in the posterior ensemble (i.e., 1 % best runs), we found that parameter correlations in the ensemble were mostly between −0.6 and 0.6, although some parameters, such as the maximum reaction rates in the stepwise denitrification process, showed relatively higher correlation coefficients (Fig. S25). The parameter values in the ensemble were overall well-constrained and normally distributed (Fig. S25). Therefore, we consider that parameter equifinality was not a severe issue here, although a formal analysis is required to make more concrete conclusions. Besides, when optimizing parameters, manually adjusting the variable weights in the objective function may have led to biased model-data errors in some variable estimates (Konak et al., 2006). Exploring better designs for multi-objective optimization (e.g., Cheng et al., 2002; Konak et al., 2006; Li et al., 2010; Nguyen et al., 2019), which involve not only temporal variation in gas emissions (e.g., CO₂, N₂O, and N₂) but also spatiotemporal variation in substrates (e.g., DOC and inorganic N) in the soil, is a topic for further studies to improve N₂O model performance.

The contributions of different biochemical pathways to N₂O emissions have been elaborately studied in wastewater systems (e.g., Chen et al., 2019; Ni et al., 2011; Peng et al., 2014), and the role of solute diffusion has long been a topic in the experimental and modeling field of soil aggregates (e.g., Kremen et al., 2005; Schlüter et al., 2018; Smith, 1980) and permeable sediments (e.g., Kessler et al., 2013). These processes, however, have not been explored well in meso- and large-scale soil modeling. The experiment investigated in this study represented a period immediately after the incorporation of liquid manure with no rainfall causing infiltration or leaching around the surfaces. With this constraint, the proposed diffusion–reaction model is capable of simulating soil conditions for constant soil moisture levels below or at the water holding capacity (i.e., no leaching). Water convection, e.g., during and after rainfall or irrigation, will require an extension of the model concept by adding hydrological processes.

“Hotspot” effects depend not only on application rate but also on the application method which defines the contact area between soil and manure and hence the exchange of gases and solutes. With a high surface-to-volume ratio (i.e., smaller hotspots), more oxygen is available for C turnover, and as a consequence the duration of anoxic conditions is too short for the significant growth of nitrifiers and coupled nitrification–denitrification. For example, Petersen et al. (1992) incubated cattle slurry as 16 or 8 mm hotspots surrounded by soil in cylinders, and a slurry corresponding to a 12 mm layer was fully mixed with soil. The N recovered as denitrification in the samples with 16 and 8 mm hotspots and the mixed sample corresponded to 40 %, 33 %, and 0.3 % of total ammoniacal N in the cattle slurry, respectively. Future applications of the model could include predicting effects of surface-to-volume ratios of manure-amended soil corresponding to different application methods (e.g., incorporated by plowing or injected), which may then be examined in experiments under field conditions. Also, the influence of soil moisture associated with modeled diffusional constraints on soil N₂O emissions could be re-assessed based on properly designed experiments.

5 Conclusion

Understanding the spatiotemporal distribution of nitrification and denitrification in agroecosystems, and particularly that associated with field-applied manure, may be crucial for reducing the uncertainty about agricultural N₂O emission estimates in process-based soil biogeochemistry models. The present modeling study, supported by an incubation experiment, gives new insights into the importance of biochemical processes in heterogeneous soil environments. The simulation results confirmed previous experimental work suggesting that nitrifier and denitrifier communities develop in closely coupled but separate niches in soils amended with manure. Denitrifiers mainly developed in the predominantly anoxic zone, whereas nitrifiers and aerobic heterotrophs proliferated around the interface between bulk soil and manure-amended soil. The stratification was affected by soil moisture and became tighter with decreasing soil water content in response to better aeration and constraints on solute diffusion. In accordance with the spatial features of bacterial communities, the saturated or near-saturated center of the manure zone was a hotspot of N₂O production by denitrification and possibly nitrifier denitrification, whereas nitrification dominated the N₂O production at the better-aerated interface according to the model. In the manure-amended soil volume where NO${}_{\mathrm{3}}^{-}$ was initially depleted, the NO${}_{\mathrm{3}}^{-}$ demand of denitrification was largely maintained by nitrification, demonstrating a strong coupling of nitrification and denitrification either in the same layer or separated by a short distance depending on soil water content. Breaking down the sources and sinks of NO${}_{\mathrm{3}}^{-}$ in the model clearly revealed the important contribution of solute diffusion to the supply of NO${}_{\mathrm{3}}^{-}$ for denitrification, and neglecting this process will significantly deteriorate the accuracy of N₂O emission estimates, as demonstrated with scenarios without the concurrent transport and reactions of NO${}_{\mathrm{2}}^{-}$, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ in the vicinity of organic hotspots. The implications of this study are of importance for soil N₂O modeling in general by suggesting that if solute transport is not included in process-based models when simulating stagnant soil conditions and if organic hotspots are present, then model-estimated N₂O emissions could be much lower than the actual emissions.