Computational Geosciences
https://doi.org/10.1007/s10596-019-09862-3

ORIGINAL PAPER

A benchmark for soil organic matter degradation under variably
saturated flow conditions
M. Jia 1 & D. Jacques 2 & F. Gérard 3 & D. Su 1 & K. U. Mayer 1 & J. Šimůnek 4
Received: 12 January 2019 / Accepted: 18 July 2019
# Springer Nature Switzerland AG 2019

Abstract
Soil contains the largest terrestrial pool of organic matter, and the cycling of organic carbon in soils plays a crucial role incontrolling atmospheric carbon dioxide (CO2) and global climate change. Although considerable progress has been made inprevious modeling studies on the fate of soil organic matter (SOM), only a few models used a process-based approach forinvestigating these strongly coupled and complex soil systems, which involve SOM oxidation, transient water flow, and masstransport processes in aqueous and gaseous phases. Typically, physically based models for water flow, as well as solute and gas
transport, are not coupled with state-of-the-art SOM degradation models. Reactive transport models (RTMs) provide a flexible
framework for implementing different SOM degradation concepts and integrating biogeochemical processes with water flow andmass transport. Given the complex nature of carbon cycling in soils coupled with flow and mass transport, code intercomparison
using well-defined benchmarks is in many cases the only practical method of model verification. The benchmark presented in thismanuscript focuses on SOM oxidation under variably saturated flow conditions. The benchmark consists of three problems
characterized by increasing complexity. The problems were solved using two different reactive transport codes, namely HP1 and
MIN3P-THCm. The first supporting problem introduces a batch-type simulation to assess kinetic networks of SOM degradation.
In the second supporting problem, transient water flow, solute transport, gas generation, and diffusive gas transport are considered. The principal problem combines the kinetic networks of SOM degradation with reactive transport under variably saturated
flow conditions, including CO2 transport from soils to the atmosphere. Simulation results for the benchmark problems demonstrate an overall excellent agreement between the two codes, building confidence in the ability of RTMs to simulate complex Ccycling in dynamic environments.
Keywords Reactive transport modeling . Model intercomparison . Soil organic matter . Benchmarking . HP1 . MIN3P-THCm

1 Introduction
Electronic supplementary material The online version of this article(https://doi.org/10.1007/s10596-019-09862-3) contains supplementarymaterial, which is available to authorized users.
* M. Jia
mjia@eoas.ubc.ca* D. Jacques
djacques@sckcen.be1

Department of Earth, Ocean and Atmospheric Sciences, The
University of British Columbia, Vancouver, BC, Canada

2

Belgian Nuclear Research Centre SCK.CEN, B-2400 Mol, Belgium

3

Supagro, Univ. Montpellier, CIRAD, INRA, IRD, Eco&Sols,
Montpellier, France

4

Department of Environmental Sciences, University of California
Riverside, Riverside, CA 92521, USA

At the global scale, soil organic matter (SOM) contains threeor four times more carbon than is found either in the atmosphere or terrestrial vegetation [1]. SOM turnover is an essential process regulating atmospheric gases, particularly the production of CO2, and it thus has a significant influence on
global climate change [2]. SOM is a highly complex material
with extreme variations in physical and chemical properties,
spatial distribution, and temporal evolution at different scales[3]. Understanding the role of SOM dynamics using direct
measurements alone is challenging, largely due to the complexity of biotic and abiotic interactions, and the heterogeneous and evolving nature of SOM. However, models canbe useful tools for integrating complex data sets and providing
their quantitative interpretations. They have thus become
promising tools for the assessment of complex carbon

Comput Geosci

turnover processes in soil [4]. Models simulating flow andtransport processes in soils are essential for assessing ecosystem services provided by the soil at different scales [5].
Accurate representation of SOM dynamics is thus an important component of such models. Abiotic factors such as watercontent and temperature can strongly affect SOM degradation
rates, whereas mass transport processes govern both the migration of dissolved organic matter (DOM) in the aqueous
phase and the CO2 transport in the gaseous phase [6]. As aresult, the SOM turnover models allow us to improve ourassessment of soil ecosystem services and to develop a better
understanding of intricate biogeochemical interactions in soils
and their effect on climate change.
During the last 80 years, different approaches have beenused to describe SOM dynamics across a broad range of temporal and spatial scales [7]. At present, an increasingly
expanding range of SOM models is being used in sustainability research and decision-making [4]. A generally accepted
concept for modeling of SOM degradation is to define a number of compartments or SOM pools, where each pool is composed of SOM with varying chemical composition and degradation characteristics [8–10]. These SOM models differ in
the number of compartments and their connectivity by having
parallel compartments, compartments in series, and feedback
between compartments. Most of the existing SOM compartment degradation models can be represented by a generalized
system of linear equations, as discussed by Sierra [11] andexamples therein. Besides connectivity, degradation formulations representing the effect of substrate concentration (i.e.,
SOM pools) and decomposers (biota) vary between models,
ranging from simple first-order kinetics without representationof decomposers to nonlinear representations including dynamic interactions with decomposers [12].
Many SOM degradation models have a fixed mathematicalstructure with a rigid, static, and problem-specific character thatlacks the flexibility and extensibility to evaluate different SOM
dynamics within the same numerical framework. However, theanalysis of different SOM dynamics (e.g., for the definition of ahierarchy of models, for model abstraction, for multi-model ensemble predictions, or for scenario and sensitivity analyses)
would benefit from a flexible modeling framework that allows
implementations of different model structures within the samemathematical framework. Only a few models combine SOM
degradation with transient water flow and relevant solute andgas transport processes, while at the same time accounting forthe influence of dynamic soil water contents and temperatures on
SOM degradation processes [13, 14]. On the other hand, typicalwater and mass transport models for soils commonly lack thestate-of-the-art formalisms to cope with SOM degradation. For
example, while the DAYCENT model focuses on assessing
multi-pool SOM degradation, it has no explicit description ofgas transport and gas exchange, which limits its ability to simulate the temporal evolution of SOM turnover processes [15, 16].

Braakhekke et al. [17] included in their model the transport ofdissolved SOM using a constant rate of advection for all solutes.
However, this model cannot be applied under all soil conditions(e.g., to grassland soils). Models developed by Fang andMoncrieff [18] described in detail the gas transport and exchange, while greatly simplifying SOM degradation. Jassal
et al. [19] developed a complex model for simulating CO2 transport in the forest floor that considered diffusion in the gaseous
phase as well as diffusion and dispersion in the aqueous phase.
However, this model only considered two SOM pools with simplified first-order degradation rates for biological CO2 production. Šimůnek and Suarez [20] also developed a CO2 production
and transport model that accounted for transient variably

Fig. 1 Schematics for organic matter degradation models: a model
without leachable soil organic matter (SOM-L) and b model withleachable soil organic matter (SOM-L). Orange lines indicate FOM
degradation to easily decomposable soil organic matter (SOM-P1),
brown lines indicate degradation to slowly decomposable soil organic
matter (SOM-P2), green lines indicate assimilation into microbial
biomass (Biomass), blue lines indicate respiration and production of
inorganic carbon species (CO2), and red lines indicate degradation to
SOM-L. The dashed black line represents kinetic exchange between
SOM-L and dissolved organic matter (DOM)

Comput Geosci
Table 1
R

Stoichiometry of SOM degradation reactions for networks 1 and 2

Reaction

FOM

SOMP1

SOMP2

CO32
−

H+

Biomass

SOML

DOM

General reactions
1

FOM ➔ SOM-P1

−1

1

0

0

0

0

0

0

2
3ǂ

Biomass ➔ SOM-P1
SOM-L ➔ DOM

0
0

1
0

0
0

0
0

0
0

−1
0

0
−1

0
1

0
0
0
0
0

−1
0
−1
0
0

0.1
−1
0.1
−1
0

0.6
0.6
0.55
0.5
0.9

1.2
1.2
1.1
1.0
1.8

0.3
0.4
0.3
0.4
0.1

0
0
0.05
0.1
−1

0
0
0
0
0

Non-explicit biomass/linear interaction/exponential interaction
4
SOM-P1 ➔ SOM-P2 + CO32− + H+ + biomass
5
SOM-P2 ➔ CO32− + H+ + biomass
6ǂ
SOM-P1 ➔ SOM-P2 + CO32− + SOM-L + H+ + biomass
7ǂ
SOM-P2 ➔ CO32− + H+ + SOM-L + biomass
8ǂ
SOM-L ➔ CO32− + H+ + biomass
Stoichiometric coefficients for products > 0 and for reactants < 0ǂOnly for network 2

saturated water flow, advective-dispersive solute transport, andCO2 production in soil, but the model structure was fixed (hardcoded). Moreover, environmental parameters (such as water content or temperature) that affect degradation rates are often calculated separately from the SOM degradation model using a flow
and transport model (e.g., Pansu et al., 2010 [21]; Batlle-Aguilar
et al., 2011 [22]) or are directly derived from observations [23].
However, this approach is not suitable for simulating the diffusion of gaseous CO2 originating from respiration of SOM orleaching and sorption of DOM because of the close coupling
between flow, transport, and SOM degradation.
In our study, we illustrate that existing generic reactive
transport models (RTMs) represent a flexible and extensible
framework for implementing different SOM degradation
model structures while providing tight integration of SOM
degradation with flow and mass transport processes. Several

Table 2

RTMs can integrate transient variably saturated flow and masstransport equations with a biogeochemical reaction network
involving equilibrium and kinetic processes. Some recent
studies have included SOM degradation dynamics into
RTMs. For example, Thayssen et al. [24] implemented themodel of Šimůnek and Suarez [20] into the RTM HP1 [25].
Within the same RTM, the SOM degradation network directly
accounted for the effect of environmental variables such as thespatiotemporal evolution of moisture and temperature and theeffect of bio-diffusion of SOM pools [26]. However, the moresophisticated nature of RTMs leads to challenges during theirimplementation and code development [27]. Additional complexities arise when model structures of SOM degradation are
defined by users via input files. To improve the reliability ofsimulation tools, model verification is essential [28]. As aresult, benchmarking of different RTMs by implementing

Overview of benchmark cases

BenchmarkType

Supporting problem 1
C1a1
B
C1a2
B
C1b1
B
C1b2
B
C1c1
B
C1c2
B
Supporting problem 2
C2a
T
C2b
T
Principal problems
P1
B-T
P2
B-T

Description

Degradation network 1 with non-explicit biomass
Degradation network 2 with non-explicit biomass
Degradation network 1 with a linear substrate-biomass term
Degradation network 2 with a linear substrate-biomass term
Degradation network 1 with an exponential substrate-biomass term
Degradation network 2 with an exponential substrate-biomass term
Transient water flow and conservative solute transport
Transient water flow, solute transport, and diffusion in the gas phase with a kinetic production term
Transient variably saturated flow, solute transport, and organic matter degradation
Transient variably saturated flow, solute transport, and organic matter degradation with a kinetic adsorption-desorption termComput Geosci
Table 3

Parameters of SOM degradation reactions for networks 1 and 2

SymbolParameter description

Values

Unit

ρb
iF
rB
Non-explicit dependency

Soil bulk density
Addition rate of FOM
Decay rate of biomass pool

1.43
1.40 × 10−89.80 × 10−8

kg dm−3 bulk
g dm−3 bulk s−1
s−1

kα,SOM-P1
kα,SOM-P2
kα,SOM-L*
Linear dependency

Rate constant for SOM-P1 degradation
Rate constant for SOM-P2 degradation
Rate constant for SOM-L degradation

3.50 × 10−7
5.80 × 10−8
3.50 × 10−8

s−1
s−1
s−1

kβ,SOM-P1
kβ,SOM-P2
kβ,SOM-L*
Exponential dependency

Rate constant for SOM-P1 degradation
Rate constant for SOM-P2 degradation
Rate constant for SOM-L degradation

3.50 × 10−6
5.60 × 10−7
4.20 × 10−7

L H2O mol−1 s−1
L H2O mol−1 s−1
L H2O mol−1 s−1

kγ,SOM-P1
kγ,SOM-P2
kγ,SOM-L*cRate constant for SOM-P1 degradation
Rate constant for SOM-P2 degradation
Rate constant for SOM-L degradation
Exponential coefficient3.50 × 10−7
5.80 × 10−8
3.50 × 10−8
12.01

s−1
s−1
s−1
L H2O mol−1

General parameters*Only for network 2

the same formulation for SOM degradation is an important
prerequisite to increasing confidence in using such models
for simulating SOM degradation processes in soil systems.
This study is aimed at presenting a benchmark problem
on C-cycling in soils that includes SOM degradation, transient water flow, and mass transport processes. We illustratesome of the issues mentioned above by accounting for (i)two different SOM degradation models, (ii) different equations describing the role of decomposers, and (iii) waterflow and mass transport processes including the transportof DOM and total inorganic carbon (TIC) and diffusive CO2transport in soil air under variably saturated flow conditions. The benchmark provides a relatively high degree ofTable 4
Symbolcomplexity involving nonlinear coupling and facilitates adirect intercomparison of different implementations of thegoverning equations in the two participating codes (i.e.,
HP1 [25] and MIN3P-THCm [29]). Three individual benchmark problems of varying complexity were solved in orderto evaluate formulations and implementations of soil carbon
turnover, and flow and transport processes. Results obtainedwith both reactive transport codes are presented and compared. Within the scope of the current benchmark, the effectof selected environmental factors (e.g., oxygen availability
and temperature dependency) on SOM degradation was notconsidered. However, RTMs allow the inclusion of suchprocesses with relative ease [26, 30].

Parameters for kinetic adsorption-desorption of SOM-L and DOM
Parameter description

Value

Unit

5.80 × 10−7
1.93 × 10−7
1.67 × 10−1
5.35 × 10−2*

s−1
mol C dm−3 bulk s−1
mol L−1 H2O
mol L−1 H2O

5.8 × 10−7
0.40
6.01
18.68*s−1
g C dm−3 bulk
L H2O mol−1
L H2O mol−1

MIN3P-THCm formalism (Eq. 15)
k1
k2
K

Rate coefficient for fractional term
Rate coefficient for the Michaelis-Menten term
Half saturation constant for the Michaelis-Menten term

HP1 formalism (Eqs. 13 and 14)First-order kinetic sorption term
ksor
Smax
Binding capacity of DOM
Binding affinity of DOM
KL*Only for the principal problem (P2)Comput Geosci

2 Mathematical formulation of SOM
degradation and model formulations

Table 5 Boundary and initial conditions—supporting problem C1, totalcomponent concentrations
Components

InitialUnit

The mathematical framework for multicomponent reactive
transport models including transient variably saturated water
flow, advective-dispersive solute and gas transport, equilibrium chemistry, and kinetic reaction networks is provided in[31] and is therefore not repeated here. The model description
given here focuses on specific equations of the benchmark.

H+a
Ca2+
CO32−
FOM

5.00
1.00 × 10−10
1.00 × 10−10

mol L−1 H2O
mol L−1 H2O
mol L−1 H2Okg C kg−1 bulk

SOM-P11.75 × 10−4

2.1 Organic matter degradation schemes

SOM-P2
Biomass
SOM-Lb
DOMb

In the benchmark, the SOM degradation model structure follows the classical compartment-type models with a number offunctionally different SOM pools. The first degradation
scheme includes four carbon pools with fresh organic matter(FOM), soil organic residues (SOM-P), microbial biomass,
and inorganic carbon (Fig. 1a). The soil organic residues are
divided between a fast-decomposing pool (SOM-P1) and aslow-decomposing pool (SOM-P2). Both aboveground plant
litter and the decay of fine roots provide a source for FOM,
which can act as a source for energy-rich substrates, increasing
the rate of SOM mineralization [17]. For the sake of model
simplicity, we make the following assumptions: (i) no distinction in degradation rates is made between aboveground litterand root litter; (ii) the input of FOM to the soil is continuous
and unlimited; and (iii) the FOM degradation rate (iF) can be
described by a zero-order kinetics rate equation and FOM
inputs contribute only to SOM-P1. The slowly decomposing
SOM pool (i.e., SOM-P2) is derived from rapidly
decomposing SOM (SOM-P1) via humification, which is inagreement with the findings of Jenkinson [8]. Although different functional groups of microbial populations coexist, onlya single pool is considered to account for microbial biomass.
In many cases, quantitative information on microbial communities is limited [9], justifying this approach. SOM decomposition involved enzymatic oxidation which produces carbon
dioxide (CO2), which is, via transport through the soil gas
phase, partly returned to the atmosphere (soil respiration).
The partitioning of CO2 between the aqueous phase and soil
gas is also included via Henry’s law. The four carbon pools
provide a simplified, yet general model structure, accounting
for organic carbon fluxes between the main components of thesystem. A set of four mass balance equations defines the degradation scheme:

a

As pH

b

Only for network 2

1.40 × 10−3

kg C kg−1 bulk
kg C kg−1 bulk
kg C kg−1 bulk
kg C kg−1 bulk
mol L−1 H2O−57.00 × 10
7.00 × 10−5
1.33 × 10−5
8.33 × 10−3

where Rx [mol L−1 H2O s−1] defines the degradation reaction rates, subscript x denotes different SOM pools (B = biomass, TIC = inorganic carbon), and rx denotes the degradation
rates described in the next section. Parameter va→b represents
a fraction of degradation products of pool a going to pool b.
Fractions va→b for a given pool a sum up to one, ensuring thatmass balance is closed. Table 1 provides the stoichiometry of
the reactions used in network 1.
The second degradation scheme (Fig. 1b), named network
2, is a variant of network 1. In addition to network 1, network
2 includes a leachable SOM pool (SOM-L), representing
SOM that is reversibly adsorbed onto soil constituents and isthus the most important source of SOM deeper in a soil profile[17]. SOM-L partitions between the solid and aqueous phases
according to a kinetic adsorption/desorption isotherm. The
desorbed fraction corresponds to DOM that is transported byadvection and dispersion throughout the soil profile.
Compared to SOM and biomass concentrations, DOM concentrations are relatively small, although DOM transport plays

Table 6 Boundary and initial conditions—supporting problem C2, totalcomponent concentrations in units of [mol L−1 H2O]
Components

RSOM−P1 ¼ i F þ rB −rSOM−P1

ð1Þ

RSOM−P2 ¼ vSOM−P1→P2 rSOM−P1 −rSOM−P2

ð2Þ

Supporting problem C2a
H+
CO32− (as a conservative tracer)
Supporting problem C2b
H+a
CO32−aRB ¼ vSOM−P1→B rSOM−P1 þ vSOM−P2→B rSOM−P2 −rB

ð3Þ

H+b
CO32−b

RTIC ¼ vSOM−P1→TIC rSOM−P1 þ vSOM−P2→TIC rSOM−P2

ð4Þ

aFor shallow soilb

For deep soil

Recharge

Initial2.00 × 10−3
1.00 × 10−3

2.00 × 10−10
1.00 × 10−10

3.70 × 10−3
1.90 × 10−3

3.70 × 10−3
1.90 × 10−3––

3.70 × 10−3
1.90 × 10−3

Comput Geosci
Table 7 Physical properties forthe flow and transport model

SymbolParameter description

Value

Unit

θr*
θs*
Sw*

Residual water content
Saturated water content
Residual saturation

0.08
0.45
0.178

m3 m−3m3 m−3–α*n*Shape parameter
Shape parameter

1.0
1.4

m−1
–

Ks
l*

Saturated hydraulic conductivity
Pore connectivity parameter
Longitudinal dispersivity
Free phase diffusion coefficient in the aqueous phase
Free phase diffusion coefficient in the gas phase

2.894 × 10−60.5
0.1
2.315 × 10−92.315 × 10−5

m s−1
–
m
m2 s−1m2 s−1DL
Dw
Da*Soil hydraulic properties are described using the van Genuchten-Mualem model [35]an important role regarding the fate of SOM [17]. This expanded degradation scheme requires modification of the massbalance equations for Network 1 and the definition of an additional mass balance equation for SOM-L, resulting in:
RSOM−P1 ¼ i F þ rB −rSOM−P1

ð5Þ

RSOM−P2 ¼ vSOM−P1→P2 rSOM−P1 −rSOM−P2

ð6Þ

RB ¼ vSOM−P1→B rSOM−P1 þ vSOM−P2→B rSOM−P2
þ vSOM−L→B rSOM−L −rB

ð7Þ

RSOM−L ¼ vSOM−P1→SOM−L rSOM−P1
þ vSOM−P2→SOM−L rSOM−P2 −rSOM−L −rS

ð8Þ

RTIC ¼ vSOM−P1→TIC rSOM−P1 þ vSOM−P2→TIC rSOM−P2
þ vSOM−L→TIC rSOM−L

ð9Þ

where rS (=−RaDOM ) is the rate describing desorption kinetics (see Section 2.3). The stoichiometry of the degradation
reactions for network 2 is given in Table 1.

2.2 Organic matter degradation rates
Three different equations describing SOM decomposition are implemented to account for the influence of
biomass: (i) linear dependence only on substrate (nonexplicit biomass), (ii) linear dependence on substrate
and biomass (linear interaction), and (iii) linear dependence on substrate and exponential dependence on biomass (exponential interaction). In the first case, the biomass is a non-limiting factor and, thus, the rate of SOM
degradation is not dependent on the simulated biomass
concentration [9]:rX ¼ k α C X

where kα is the first-order degradation rate [s−1] and Cx isthe concentration of the substrate [mol L−1 H2O]. However,this simple approximation is often unrealistic, because thedegradation of recalcitrant SOM is more strongly affected by
the quantity of enzymes associated with biomass than by theconcentration of the substrate [23]. Thus, a linear dependence
on both quantities of microbial biomass and substrate, describing a co-limitation by decomposers and substrate, provides analternative formulation [32]:ð10ÞTable 8 Boundary and initial conditions—principal problem, totalcomponent concentrations
Components

Recharge

InitialUnit

Shallow soil
H+
Ca2+

3.70 × 10−3
1.00 × 10−10

3.70 × 10−3
1.00 × 10−10

mol L−1 H2O
mol L−1 H2O

1.90 × 10−3
3.50 × 10−4
1.75 × 10−4
7.00 × 10−5
7.00 × 10−5
1.33 × 10−5
2.33 × 10−3

1.90 × 10−3
3.50 × 10−4
1.75 × 10−4
7.00 × 10−5
7.00 × 10−5
1.33 × 10−5
2.33 × 10−3

mol L−1 H2Okg C kg−1 bulk
kg C kg−1 bulk
kg C kg−1 bulk
kg C kg−1 bulk
kg C kg−1 bulk
mol L−1 H2O

–
–
–
–
–
–
–
–
–

3.70 × 10−3
1.00 × 10−10
1.90 × 10−3
4.30 × 10−10
4.30 × 10−10
4.30 × 10−10
4.30 × 10−10
4.30 × 10−10
1.00 × 10−10

mol L−1 H2O
mol L−1 H2O
mol L−1 H2Okg C kg−1 bulk
kg C kg−1 bulk
kg C kg−1 bulk
kg C kg−1 bulk
kg C kg−1 bulk
mol L−1 H2O

CO32−
FOM
SOM-P1
SOM-P2
Biomass
SOM-LaDOMa
Deep soil
H+
Ca2+
CO32−
FOM
SOM-P1
SOM-P2
Biomass
SOM-LaDOMa
aOnly for network 2

Comput GeoscirX ¼ k β C X C B

where kγ is the first-order degradation term [s−1] and c [LH2O mol−1] is a shape parameter. The decay of microbial
biomass (rB) is always modeled as first-order kinetics.

ð11Þ

where kβ is the second-order degradation rate [L H2O
mol−1 s−1] and CB is the concentration of the microbial
biomass [mol L−1 H2O]. Although this formulation considers the quantity of microbial biomass, this linear dependence does not adequately describe the non-linear effects ofbiomass concentrations on SOM degradation. Alternatively,
we can consider an exponential formulation that follows theformulation of Wutzler and Reichstein [12]:2.3 Sorption of dissolved organic matter
Kinetic adsorption-desorption reactions involving SOM-L
and DOM are described using a first-order kinetic rate equation:

rS ¼ −k S C SOM−L −C eq
SOM−L

ð12Þ


rX ¼ k γ 1−e−cC B C X

0.05
0

500

1000
Time[day]

0.15
0.1
0.05

0.3

1000
Time [days]

1500

0.2
0.15
0.1
0.05
0

500

1000
Time [day]1500

Fig. 2 Results for supporting problem C1 based on SOM degradation for
networks 1 (a–c) and 2 (d–f) for linear substrate dependence only, linear
biomass-substrate interaction, and exponential biomass-substrate

0.25
0.2
0.15

0.1

0.1

0.05

0.05

0

0
500

1000
Time[day]

1500
0.3

(e) Linear interaction

0.2

0.2

0.1

0.1

0

0
0

0.3

(c) Exponential interaction

0.25

0

0.15

0

SOM / Biomass [g dm-3 bulk]

0.2

500

0.2

0.3

(b) Linear interaction

0

DOM
SOM-P1
SOM-P2
Biomass
SOM-L

-3

1500

0.25

0

0.25

DOM [mol L-1 water]

0.1

Network 2

DOM [mol L-1 water]

0.15

0.3

(d) Non-explicit

500

1000
Time[day]

1500
0.3

(f) Exponential interaction

0.2

0.2

0.1

0.1

0

DOM [mol L-1 water]

0.2

SOM / Biomass [g dm bulk]

Network 1
SOM-P1
SOM-P2
Biomass

0.3
SOM / Biomass [g dm-3 bulk]

0.3(a) Non-explicit

0.25

0

SOM / Biomass [g dm-3 bulk]

where kS is the first-order rate constant [s−1] and C eq
SOM−L isthe equilibrium concentration [mol L−1 H2O]. The equilibrium

SOM / Biomass [g dm-3 bulk]

SOM / Biomass [g dm-3 bulk]

0.3

ð13Þ

0
0

500

1000

1500

Time [day]interaction. Note that symbols represent HP1 results, while solid lines
represent MIN3P-THCm resultsComput Geosciconcentration is calculated based on the Langmuir isotherm[33]:
K L S max C DOM
C eq
SOM−L ¼
1 þ K L C DOM

ð14Þ

where KL is the coefficient of binding affinity [L H2O mol−1]
and Smax is the maximum adsorption capacity [mol L−1 H2O].
The formulation of sorption of DOM was modified for the
MIN3P-THCm simulations to allow its inclusion based on theexisting geochemical reaction network available in the code [29].
To this end, Eqs. 13 and 14 were combined and written in themathematically equivalent form of a Michaelis-Menten rate expression:rS ¼ −k 1 C SOM−L þ k 2

C DOM
K þ C DOM

ð15Þ

where k1 = kS and k2 = kS × Smax. K is the half saturation
constant in the Michaelis-Menten term and equals 1/KL.

2.4 Mass conservation, solute, and gas transport
In addition to degradation processes, transport of geochemical
species is considered in the benchmark. The mass conservation equation for component j present in the aqueous andgaseous phases and described in terms of total component
concentrations takes the form:
i ∂h
i
∂h
S a φT aj þ
S g φT gj
∂t
∂t
h
i h
i h
i
a
þ ∇⋅ qa T j −∇⋅ S a φDa ∇T aj −∇⋅ S g φDg ∇T gj −Q j¼0where T aj [mol L−1 H2O] defines the total aqueous component
concentrations; T gj [mol L−1 gas] defines the total component
concentrations in the gas phase; φ [m3 void m−3 bulk] is theporosity; Sa [m3 H2O m−3 void] and Sg [m3 gas m−3 void] aresaturations of the aqueous and gas phases, respectively; qa [m
0.45

0.45

MIN3P-THCm
HP1

0.4

Water content [-]Water content [-]

0.4
0.35
0.3
0.25

0.35
0.3
0.25
0.2

0.2

d = 3 cm
0.15

ð16Þ

0

100

200

0.15

300

d = 6 cm
0

100

0.45

0.4

0.4

Water content [-]Water content [-]

0.45

0.35
0.3
0.25
0.2
0.15

0.3
0.25
0.2

d = 10 cm
0

100

200

0.15

300

d = 25 cm
0

100

200

300

Time [days]

0.45

0.45

0.4

0.4

Water content [-]Water content [-]

300

0.35

Time [days]

0.35
0.3
0.25
0.2
0.15

200

Time [days]

Time[days]

0.35
0.3
0.25
0.2

d = 50 cm
0

100

200

Time [days]

300

0.15

d = 100 cm
0

100

200

300

Time [days]

Fig. 3 Results for supporting problem C2a showing transient evolution of water contents at selected depths (3, 6, 10, 25, 50, and 100 cm from top tobottom)Comput Geoscis−1] is the Darcy flux vector; and Da [m2 s−1] is the hydrodynamic dispersion tensor applicable to all dissolved species in thesystem. The tortuosity implicitly contained in the pore gas diffusion coefficient Dg is calculated based on the classical Millington
and Quirk [34] model. Finally, Qj are source-sink terms [moldm−3 bulk s−1] due to kinetically controlled reactions and fluxes
across the boundaries.

geochemical equilibrium and kinetic reactions, as well as about
numerical methods used, are given in [31]. Input and databasefiles, as well as selected output files for the two codes HP1 and
MIN3P-THCm, are provided as text files as part of the SI.

3 Benchmark definitions
and parameterization

2.5 Reactive transport codes

0.006

0.004

0.002

0.000

d = 3 cm
0

100

200

0.008

-1

MIN3P-THCm
HP1

The benchmark consists of a principal problem and two
supporting problems (Table 2). The first supporting problem(C1) (Section 3.1) focuses on the simulation of kinetic reaction
networks for SOM degradation. The second supporting problem(C2) (Section 3.2) considers transient water flow and solute transport, as well as gas generation and gas transport. Finally, theprincipal problem (P) (Section 3.3) combines the processes ofthe supporting problems C1 and C2 in a reactive transport modeling framework, explicitly describing SOM degradation undervariably saturated flow conditions.

Conservative solute [mol L H2O]

0.008

-1

Conservative solute [mol L H2O]

The benchmark problems were simulated using two reactive
transport codes: HP1 [25, 26] and MIN3P-THCm [29]. The mostimportant difference between HP1 and MIN3P-THCm is thesolution technique of the coupled reactive-transport equations.
While HP1 solves the coupled reactive-transport equations using
sequential non-iterative approach (SNIA), MIN3P-THCm solves
them simultaneously using the global implicit method (GIM).
Details about the governing equations for variably saturated water flow, advective-dispersive solute transport, gas diffusion, and

300

0.006

0.004

0.002

d = 6 cm

0.000
0

100

-1

0.006
0.005
0.004
0.003
0.002
0.001
0.000

d = 10 cm
0

100

200

300

0.005
0.004
0.003
0.002
0.001

d = 25 cm

0.000
0

100

0.002

0.001

d = 50 cm
200

Time [days]

300

Conservative solute [mol L-1 H2O]

-1

Conservative solute [mol L H2O]

0.003

100

200

300

Time [days]

0.004

0

300

0.006

Time [days]

0.000

200

Time [days]
Conservative solute [mol L H2O]

-1

Conservative solute [mol L H2O]Time [days]

0.004

0.003

0.002

0.001

0.000

d = 100 cm
0

100

200

300

Time [days]

Fig. 4 Results for supporting problem C2a showing transient evolution of conservative solute concentrations at selected depths (3, 6, 10, 25, 50, and100 cm from top to bottom)

Comput Geosci

The first supporting problem (C1) verifies the implementationof two different SOM degradation networks, each evaluated forthe three decomposer equations defined above (Eqs. 10–12) ina batch simulation, i.e., without transport. The parameters forthe two degradation networks are loosely based on literature
data [9, 17, 22] and are provided in Tables 3 and 4. Initialconditions for the supporting problem C1 are provided inTable 5. The equilibrium reaction network for inorganic C speciation is provided in the SI (Section S.1). All geochemical
reactions and associated thermodynamic data for aqueous components (i.e., H+, Ca2+, CO2−
3 , and DOM), secondary aqueous
species (i.e., OH−, HCO−3 , H2CO3, and CaHCOþ
3 ), and equilibrium gas partitioning (i.e., CO2) with water are provided in
the SI (Tables S.1–S.3). The simulations were run for 1825 days
(i.e., 5 years) with a maximum time step of 0.1 days.

While transient variably saturated water flow is a highly nonlinear problem and thus difficult to solve numerically, thenumerical solution of the advection-dispersion equation foradvection-dominated transport is subject to numerical dispersion. To ensure that flow and transport processes are solved
accurately under dynamic flow conditions, model verificationis conducted first without the added complexity of SOM turnover models. The second supporting problem (C2) was designed to focus on the simulation of transient water flow, conservative solute transport, gas generation, and diffusive gastransport. The supporting problem C2 is subdivided into twoparts: (i) transient water flow and conservative solute transport(C2a) and (ii) transient water flow, solute transport, gas generation, and gas diffusion (C2b). A summary of the initial andboundary conditions for the supporting problems C2a and

0.04

(a)

MIN3P-THCm
HP1

0.03
0.02
0.01
0
0

100

200

300

Actual water outflux [m day-1]

3.2 Transient water flow, solute transport, and gas
generation (supporting problem C2)

Actual water influx [m day-1]

3.1 Kinetic networks for SOM degradation(supporting problem C1)

(b)

0.04
0.03
0.02
0.01
0
0

100

0.05

(c)

0.04
0.03
0.02
0.01
0
0

100

200

300

(d)

0.006
0.004
0.002
0
0

100

0.8
0.6
0.4
0.2
0
300

Fig. 5 Results for supporting problem C2a showing transient evolutionof a actual influx of water across the surface and b actual outflux of water
across both the surface and boundary, c actual influx of a conservative
solute across the surface and d actual outflux of a conservative solute

Cum. solute outflux [mol m-2]

Cum. solute influx [mol m-2]

(e)

200
Time [days]

200

300

Time [days]

1

100

300

0.008

Time [days]

0200
Time [days]

Actual solute outflux [mol m-2 day-1]

Actual solute influx [mol m-2 day-1]

Time [days]

0.3

(f)

0.2

0.1

0
0

100

200
Time [days]

300

across the bottom of the domain, and e cumulative influx of aconservative solute across the surface and f cumulative outflux of aconservative solute across the bottom of the domain

Comput Geosci

C2b is provided in Table 6. All geochemical parameters andassociated thermodynamic data for the supporting problems
C2a and C2b are provided in the SI (Tables S.1–S.3). Gas
generation due to soil respiration for the supporting problem
C2b is included via a zero-order source term.
3.2.1 Transient water flow and conservative solute transport(supporting problem C2a)3.2.2 Transient water flow, solute transport, gas generation,and gas diffusion (supporting problem C2b)
The supporting problem C2b builds on C2a and adds gas
diffusion and the generation of CO2 resulting from soil
respiration. The flow and transport models are similar to
C2a with some modifications to the boundary conditions
and initial solute concentrations. A mixed-type boundary
condition is similar to the third-type boundary condition

CO2 (aq) [mol L-1 H2O]

MIN3P-THCm
HP1

0.0025-1CO2 (aq) [mol L H2O]

The supporting problem C2a focuses on the verificationof transient water flow and conservative solute transport
in a partially saturated, one-dimensional soil profile with adepth of 100 cm. Table 7 summarizes all parameters forthe flow and transport simulations. The top boundary forflow is defined as a prescribed flux based on daily potential water fluxes (calculated as the difference between
rainfall and potential evaporation). A free drainage boundary condition is applied at the bottom. The initial condition is set to a pressure head of − 150 cm. For solutetransport, a third-type (Cauchy) boundary condition isapplied at the top boundary for recharge conditions(Table 6). During rainfall events, the top boundary is defined as a prescribed flux for solute transport, while underevaporative conditions, the top boundary changes and isclosed for solute transport, leading to solute being
retained and causing concentration build-up. At the bottom, a free exit boundary condition is set. Initially, soluteconcentrations are negligible throughout the profile. The
simulation was run for 365 days with a maximum time
step of 0.025 days.

0.002

0.00150.00250.002

0.0015d = 3 cm0

d = 6 cm
100

200

300

0

100

Time [days]

CO2 (aq) [mol L-1 H2O]

0.0025-1CO2 (aq) [mol L H2O]

300

0.003

0.003

0.002

0.00150.00250.002

0.0015d = 25 cm

d = 10 cm
0

100

200

0

300

100

Time [days]

200

300

Time [days]

0.003

0.003

CO2 (aq) [mol L-1 H2O]

CO2 (aq) [mol L-1 H2O]

200

Time [days]

0.00250.002

0.00150.00250.002

0.0015d = 100 cm

d = 50 cm
0

100

200

Time [days]

300

0

100

200

300

Time [days]

Fig. 6 Results for supporting problem C2b showing transient evolution of aqueous CO2 concentrations at selected depths (3, 6, 10, 25, 50, and 100 cm
from top to bottom)

Comput Geosci

3.3 Principal problems

0.03
0.02
0.01
0
100

200
Time [days]

0.06
0.04
0.02
0
0

100

200
Time [days]

300

(e)

0.08
0.06
0.04
0.02
0
0

100

200
Time [days]

0.03
0.02
0.01
0
0

(c)

0.08

(b)

0.04

300

300

Fig. 7 Results for supporting problem C2b showing transient evolutionof a actual influx of water across the surface and b actual outflux of water
across both the surface and bottom of the domain, c actual influx of a

Actual solute outflux [mol m-2 day-1]

Actual solute influx [mol m-2 day-1]

0

Actual gas influx [mol m-2 day-1](a)MIN3P-THCm
HP1

0.04

Actual gas outflux [mol m-2 day-1]

Actual water influx [m day-1]

The objective of the principal problem is to verify the integration of biogeochemistry, flow, and transport for the case of anon-explicit biomass rate equation only (Eq. 10), because thisform is the most common way of implementing SOM degradation models [12]. The first principal problem (P1) includes
the processes of the supporting problems C1a (non-explicitcase) and C2b, using the same parameters, except that thefirst-order kinetic gas production term of C2b is replaced by
the SOM degradation network. The second principal problem(P2) is similar to P1, but also includes kinetic adsorption ofDOM, following reaction network 2. The relevant adsorption
parameters (i.e., the rate constant and the half-saturation constant for the Michaelis-Menten term) were slightly modified toincrease the sensitivity of the results to changes in water contents. The principal problems include advective-dispersive
solute transport, diffusive gas transport, kinetic reaction networks for SOM degradation, kinetic sorption/desorption ofDOM (in the case of P2), aqueous complexation, and gas
exchange reactions. The principal problems provide an opportunity to evaluate the ability of generic reactive transport
codes to simulate SOM degradation under transient variably
Actual water outflux [m day-1]

in C2a for solute transport. However, influx and outfluxof CO2 are also permitted through the gas phase (Table 6).
The initial conditions throughout the domain are identicalto the boundary condition (Table 6). In addition, soil respiration is considered in the upper 30 cm of the soil profile. The CO2 generation term is described by a zero-order
reaction with a rate constant of 5.79 × 10−10 mol dm−3
bulk s−1. The simulation was run for 365 days with amaximum time step of 0.0005 day. Note that C2b provides the foundation for the principal problem, since both
the decomposition of substrates and microbial biomass
produce CO2, resulting in the release of CO2 to the atmosphere, depending on soil physical properties and moisture dynamics.

100

200
Time [days]

300

(d)

0.08
0.06
0.04
0.02
0
0

100

200
Time [days]

300

(f)

0.08
0.06
0.04
0.02
0
0

100

200
Time [days]

300

solute across the surface and d actual outflux of a solute across the bottomof the domain, and e actual influx of gas across the bottom and f actual
outflux of gas across the surface of the domain

1.0
0.0

-2

100

200
Time [days]

1.0
0.0
100

200

300

Time [days]

4.0

2.0

0.0
200

2.0

0

(c)

100

(b)

3.0

300

6.0

0

Cum. solute outflux [mol m-2]

2.0

0Cum. gas influx [mol m ](a)

MIN3P-THCm
HP1

3.0

Cum. gas outflux [mol m-2]

Cum. solute influx [mol m-2]

Comput Geosci

(d)

6.0

4.0

2.0

0.0

300

0

100

Time [days]

200

300

Time [days]

Fig. 8 Results for supporting problem C2b showing transient evolution of a cumulative influx of a solute across the surface and b cumulative outflux of asolute across the bottom of the domain and c cumulative influx of gas across the bottom and d cumulative outflux of gas across the surface of the domain-3

SOM / Biomass [g dm bulk]

SOM-P1
Biomass
SOM-P2

0.3
0.2
0.1

d = 3 cm
0

100

200
Time[day]

0.3
0.2
0.1

d = 10 cm
0

100

200
Time[day]300

Fig. 9 Results for principal problem P1 for SOM degradation network 1(no biomass dependence) showing transient evolution of SOM andbiomass concentrations at selected depths (3, 6, 10, and 25 cm from top

0.3
0.2
0.1
0

300

0.4

0

SOM / Biomass [g dm-3 bulk]

0.4

0.4

0

principal problems are provided in Table 8 and SI
Tables S.1–S.3, respectively. The simulation time and time
stepping are the same as in the supporting problems C2.

SOM / Biomass [g dm-3 bulk]

SOM / Biomass [g dm-3 bulk]

saturated flow conditions. The associated kinetic degradation
and adsorption parameters are provided in Tables 3 and 4.
Initial conditions and thermodynamic parameters for thed = 6 cm
0

100

200
Time [days]

300

100

200
Time [days]

300

0.4
0.3
0.2
0.1
0

d = 25 cm
0

to bottom). Note that symbols represent HP1 results, while solid lines
represent MIN3P-THCm results

200Time[day]d= 10 cm

0.3

0.3

0.2

0.2

0.1

0.1

0

0

100

200
Time[days]

0

300

0.005

MIN3P-THCm
HP1

0.004
0.003
0.002

0.1

0

200
Time [days]

0.003
0.002

0.1

0.1

0.002

d = 50 cm
0

100

200
Time [days]

0

100

300

200
Time [day]

0

300

0.005
0.004
0.003
0.002
0.001

d = 6 cm
0

100

200
Time [days]

300

100

200
Time [days]

300

200
Time [days]

300

0.005
0.004
0.003
0.002
0.001 d = 25 cm
0

-1

0.003

0.001

0

6, 10, and 25 cm from top to bottom). Note that symbols represent HP1
results, while solid lines represent MIN3P-THCm results

CO2 (aq) [mol L H2O]

0.004

0.3

0.2

300

0.005

0

300

0.2

d = 10 cm
200
Time [days]

200
Time[day]d = 25 cm

0.3

-1

0.004

100

100

300

0.005

0.001

0

DOM [mol L-1 water]

0.1

CO2 (aq) [mol L H2O]

-1

CO2 (aq) [mol L H2O]

100

0

-1

0.2

0.001 d = 3 cm
0

CO2 (aq) [mol L H2O]

0.2

CO2 (aq) [mol L-1 H2O]

CO2 (aq) [mol L-1 H2O]

Fig. 10 Results for principal problem P2 based on SOM degradation
network 2 (no biomass dependence) showing transient evolution ofDOM, SOM, biomass, and SOM-L concentrations at selected depths (3,SOM / Biomass [g dm-3 bulk]

0

300

0.3

-1

100

d =6 cm
0.3

DOM [mol L water]0

SOM / Biomass [g dm bulk]

0.1

0.1

0

0.2

-1

0.2

0.3
DOM [mol L water]

DOM
SOM-P1
SOM-P2
Biomass
SOM-L-3

d = 3 cm

0.3

DOM [mol L-1 water]

SOM / Biomass [g dm-3 bulk]

SOM / Biomass [g dm-3 bulk]

Comput Geosci

0.005
0.004
0.003
0.002
0.001

d = 100 cm
0100

Fig. 11 Results for principal problem P1 showing transient evolution of aqueous CO2 concentrations at selected depths (3, 6, 10, 25, 50, and 100 cm
from top to bottom)

Comput Geosci

4 Results and discussion
4.1 Kinetic networks for SOM degradation(supporting problem C1)

0.005

CO2 (aq) [mol L-1 H2O]

MIN3P-THCm
HP1

0.004
0.003
0.002
0.001

d = 3 cm
100

0.004
0.003
0.002
0.001

200
Time [days]

0.002
0.001 d = 6 cm

0.003
0.002
0.001 d = 50 cm
100

200
Time [days]

300

200
Time [days]

300

100

200
Time [days]

300

100

200
Time [days]

300

0.004
0.003
0.002
0.001 d = 25 cm
0

-1

0.004

100

0.005

300

0.005

0

0.003

0

d = 10 cm
100

0.004

300

0.005

0

CO2 (aq) [mol L-1 H2O]

200
Time [days]

CO2 (aq) [mol L-1 H2O]

CO2 (aq) [mol L-1 H2O]

0

0.005

CO2 (aq) [mol L H2O]

CO2 (aq) [mol L-1 H2O]

The evolution of SOM-P1, SOM-P2, and microbial biomassin network 1 is similar for various formulations of degradation
as a function of its dependency on biomass (Fig. 2a–c). Two
distinct periods can be observed during the evolution of microbial biomass. The concentration of microbial biomass increases rapidly during the first 100 days, followed by a gradual decline and an approach towards quasi-steady-state conditions after 500 days. The rapid accumulation of biomass occurs because a fraction of the decomposed substrates (i.e.,
SOM-P1 and SOM-P2) is rapidly assimilated by the microbial
biomass at high substrate concentrations, although the decayof microbial biomass concurrently leads to the production of
SOM-P1. After 500 days, the rate of substrate assimilation
gradually decreases due to a decline of substrate concentrations. Also, SOM-P1 and SOM-P2 reach constant concentrations after 500 days, when quasi-steady-state conditions arereached. Wutzler and Reichstein [12] also found that SOM
pools eventually reach steady-state levels for the degradation
rate formulations (Eqs.10–12) used in this study. Similar results are obtained for network 2 (Fig. 2d–f). As a result of thefirst-order implementation of the sorption process (Eqs.13–15), the concentrations of SOM-L and DOM slowly approach
equilibrium conditions as shown in Fig. 2d–f. In addition, Fig.
2 illustrates the effect of different degradation rate formulations as a function of the dependency on biomass on approaching the quasi-steady-state conditions. For example,
the concentrations of SOM-P1 are slightly higher for the
non-explicit formulation in comparison to the linear and exponential equations. These differences are related to the dependency of the FOM assimilation rate on SOM degradation
rates. More detailed analysis of this topic is beyond the scopeof the current contribution, and the reader is referred to [12].
Figure 2 shows a near perfect agreement between the results obtained by HP1 and MIN3P-THCm, with almost identical transient evolution of all organic pools for networks 1 and2. The results for the inorganic pool also showed excellent
agreement (results not shown).

0.005
0.004
0.003
0.002
0.001 d = 100 cm
0

Fig. 12 Results for principal problem P2 showing transient evolution of aqueous CO2 concentrations at selected depths (3, 6, 10, 25, 50, and 100 cm
from top to bottom)

Comput Geosci

0.03
0.02
0.01
0
0

100

200
Time [days]-2

Actual gas influx [mol m-2 day-1]

0.06
0.04
0.02
0
0

100

200
Time [days]

300

(e)

0.30
0.20
0.10
0.00
0

100

200
Time [days]

-1

0.03
0.02
0.01
0
0

(c)

0.08

(b)

0.04

300

-1

Actual solute influx [mol m day ](a)

MIN3P-THCm
HP1

300

Fig. 13 Results for principal problem P1 showing transient evolution of a
actual influx of water across the surface and b actual outflux of water
across both the surface and bottom of the domain, c actual influx of a

Actual solute outflux [mol m-2 day-1]

0.04

Actual gas outflux [mol m-2 day-1]-1Actual water influx [m day ]Simulation results for supporting problem C2a show that transient evolution of water contents (Fig. 3) and conservative
solute concentrations (Fig. 4) are affected by rainfall and evaporation (Fig. 5a, b). Significantly different temporal responses
in water contents and solute concentrations can be observed in
the shallow (0–30 cm) and deep (30–100 cm) soil. Whilehigh-frequency temporal variability can be observed in shallow depths, low-frequency variability can be observed indeeper horizons. This behavior is the result of highfrequency temporal variability of rainfall and evaporation that
directly propagates to water content, leading to substantial
fluctuations of solute concentrations in the shallow soil, but
a more muted response at greater depths. Both models agree
well in describing both pronounced fluctuations in water contents and conservative solute concentrations near the surface
as well as in more dampened responses deeper in the soil
profile (Figs. 3 and 4). Moreover, there is a near perfectagreement between the two codes in simulated actual conservative solute influxes and outfluxes (Fig. 5c, d), and cumulative conservative solute influxes and outfluxes (Fig. 5e, f)
across the surface and bottom of the soil profile.
The supporting problem C2b considers soil respiration, exsolution of CO2, diffusive gas transport of CO2, and CO2
emissions to the atmosphere. In the simulations, diffusion isconsidered to be the dominant transport process in the gas
phase, and respiration is restricted to the shallow soil (0–30cm). The simulation results for water contents are identical to
the results for supporting problem C2a (results not shown).
The concentrations of aqueous CO2 slightly increase withdepth (Fig. 6), which can be attributed to the aqueous CO2transported downwards with infiltrating water and longer diffusion paths towards the ground surface. The largest mass
fluxes for CO2 occur in the gas phase and lead to emissions
to the atmosphere (Figs. 7f and 8d). An overall excellent
agreement is reached between the results of HP1 and
MIN3P-THCm for aqueous CO2 concentrations in the soil
profile (Fig. 6), actual water fluxes (Fig. 7a–b), and actual/
Actual water outflux [m day ]4.2 Transient water flow, solute transport, and gas
generation (supporting problem C2)

100

200
Time [days]

300

(d)

0.08
0.06
0.04
0.02
0
0

100

200
Time [days]

300

(f)

0.30
0.20
0.10
0.00
0

100

200
Time [days]

300

solute across the surface and d actual outflux of a solute across the bottomof the domain, and e actual influx of gas across the bottom and f actual
outflux of gas across the surface of the domain

Comput Geosci

1.0

0.0

Cum. gas influx [mol m-2]

0

100

200
Time [days]

300

(c)

60
40
20
0
0

100

200
Time [days]

300

3.0

(b)

2.0

1.0

0.0
0

Cum. gas outflux [mol m-2]

2.0

The evolution of different SOM pools shows similar
trends as was seen in the supporting problems C1.
However, the evolution of DOM is smoothened due towater content fluctuations. As depicted in Figs. 9 and
10 for networks 1 and 2, respectively, both models arecapable of describing SOM degradation during variably
saturated flow conditions and results agree well. CO2
produced by SOM degradation is partially released tothe atmosphere via diffusion in the gas phase (Figs.
13f and S.4f), while CO2 is also transported to the bottom of the soil profile by infiltrating water (Figs. 11
and 12). Note that CO2 concentrations close to the soil
surface are relatively low due to the surface boundary
condition, which mimics rapid CO2 release to the atmosphere. Agreement between the results of HP1 and
MIN3P-THCm is excellent for water content (Figs. S.2
and S.3) in the soil profile, actual water fluxes (Figs.
13a, b and S.4a, b), actual aqueous/gaseous CO2 fluxes(Figs. 13c–f and S.3c–f), and cumulative aqueous/gaseous CO2 fluxes (Figs. 14a–d and S.5a–d) acrossthe surface and bottom of the soil profile. However,slight differences are observed for actual gaseous CO2
outfluxes (Figs. 13f and S.4f). Analogous to the
supporting problem C2b, the reason for these differences can be attributed to small operator splitting errors
in HP1. However, differences are again small and do
not affect the interpretation of modeling results.

-2

(a)

MIN3P-THCm
HP1

4.3 Principal problems

Cum. solute outflux [mol m ]

3.0

-2

Cum. solute influx [mol m ]cumulative CO2 fluxes in the aqueous and gaseous phases
across the surface and bottom of the soil profile (Figs. 7c–f
and 8, respectively). Slight discrepancies exist in actual gaseous CO2 outfluxes and cumulative aqueous CO2 influxes(Figs. 7f and 8a, respectively) across the soil-atmosphereboundary. However, the effect on the mass balance is negligible, since differences in solute influxes are cancelled out by
the differences in gas effluxes. To assess these differences inmore detail, a sensitivity analysis on the effect of time
stepping was conducted with the HP1 and MIN3P-THCm
codes (see SI, Fig. S.1). The maximum time step was varied
by one order of magnitude. The sensitivity analysis shows thatthe results for MIN3P-THCm do not change significantly as afunction of time stepping, while the HP1 simulations show
more pronounced differences as a function of time stepping.
In addition, the sensitivity analysis shows that the results of
HP1 approach those of MIN3P-THCm for smaller time steps.
Therefore, the main reason for the small observed differences
can be attributed to the different coupling techniques (GIM in
MIN3P-THCm and SNIA in HP1), which introduces small
yet acceptable operator splitting errors in the HP1 solution
for larger time steps. The impact of these differences on theinterpretation of results is negligible. In addition, the current
benchmark study does not account for the ingress of atmospheric O2. However, this is not a limitation of RTMs. Both
HP1 and MIN3P-THCm are able to simulate atmospheric O2
ingress and soil CO2 egress, as was demonstrated in a previous
benchmarking study on the generation and attenuation of acidrock drainage [36].

100

200
Time [days]

300

(d)

60
40
20
0
0

100

200
Time [days]

300

Fig. 14 Results for principal problem P1 showing transient evolution of a cumulative influx of a solute across the surface and b cumulative outflux of asolute across the bottom of the domain and c cumulative influx of gas across the bottom and d cumulative outflux of gas across the surface of the domain

Comput Geosci

5 Conclusions

4.

A benchmark problem was established to evaluate the implementation of SOM degradation models under variably saturated flow conditions in generic reactive transport models. Tworeactive transport codes were involved in the benchmark study(namely HP1 and MIN3P-THCm). Both codes were capableof capturing the degradation of SOM and associated masstransport processes. Overall, excellent agreement was obtained between the two codes. Small differences in simulated
results can be attributed to different solution methods (i.e.,
SNIA and GIM). However, differences between the simulation results are practically negligible and do not affect theinterpretation of modeling results. Consequently, good agreement between the results provides confidence in the use ofmulticomponent reactive transport codes for implementing
C-cycling models in dynamic environments. In comparison
to classical SOM degradation models, multicomponent reactive transport codes have several distinct advantages such as
the providing of a flexible framework for implementing different SOM degradation model structures and the tightcoupling with flow and mass transport processes. In addition,
reactive transport codes also allow accounting for (i) soil statevariables in degradation models (e.g., Jacques et al., 2018[26]), (ii) geochemical reactions for pH and redox buffering,(iii) soil interactions with a single root and root architecture(e.g., Gérard et al., 2017 [37], Nowak et al., 2006 [38]), (iv)the effect of SOM on contaminant mobility (e.g., Leterme
et al., 2014 [39]), (v) bioturbation (e.g., Jacques et al., 2018[26]), and (vi) soil structure and its effects on SOM degradation and the transport of O2 (e.g., Kuka et al., 2007, [40];Keiluweit et al., 2016 [41]). The results presented here indicate that computer codes capable of solving this benchmark
provide a promising platform for addressing similar problems
involving a range of climate conditions, soil types, soil moisture dynamics, and biogeochemical reactions.

5.

Funding information Funding for this research was provided by
Agriculture and Agri-Food Canada through the project Valuing diversity
in agro-ecosystems: The interplay of natural habitat, integrated BMPs,
and field cropping systems on GHG emissions and carbon stocks,
supporting M. Jia with a Research Assistantship and through an
NSERC (Natural Sciences and Engineering Research Council of
Canada) discovery grant held by K. Ulrich Mayer.

References
1.
2.

3.

Paul, E. A.: Soil Microbiology, Ecology and Biochemistry, 3rd Ed.
Academic Press, M.A. (2007).
Jenkinson, D.S., Adams, D.E., Wild, A.: Model estimates of CO2
emissions from soil in response to global warming. Nature. 351,
304–306 (1991)
Dungait, J.A.J., Hopkins, D.W., Gregory, A.S., Whitmore, A.P.:
Soil organic matter turnover is governed by accessibility not recalcitrance. Glob Chang Biol. 18, 1781–1796 (2012)

6.

7.

8.
9.

10.
11.

12.

13.

14.

15.

16.

17.

18.

19.

20.

21.

22.

23.

Campbell, E.E. and Paustian, K.: Current developments in soil organic matter modeling and the expansion of model applications: areview. Environ Res Lett 10. (2015).
Vereecken, H., et al.: Modeling soil processes: review, key challenges, and new perspectives. Vadose Zone J 15(5), (2016).
Adair, E., Parton, W., Del Grosso, S., Silver, W., Harmon, M., Hall,
S., Burke, I., Hart, S.: Simple three-pool model accurately describes
patterns of long-term litter degradation in diverse climates. Glob
Change Biol. 14, 2636–2660 (2008)
Manzoni, S., Porporato, A.: Soil carbon and nitrogen mineralization: theory and models across scales. Soil Biol Biochem. 41,
1355–1379 (2009)
Jenkinson, D.: The turnover of organic carbon and nitrogen in soil.
Philos Trans R Soc Lond B. 329, 361–368 (1990)
Porporato, A., D’Odorico, P., Laio, F., Rodriguez-Iturbe, I.:
Hydrologic controls on soil carbon and nitrogen cycles. I.
Modeling scheme. Adv Water Resour. 26, 45–58 (2003)
Bosatta, E., Ågren, G.I.: Theoretical-analysis of degradation of heterogeneous substrates. Soil Biol Biochem. 17, 601–610 (1985)
Sierra, C.A., Müller, M., Trumbore, S.E.: Models of soil organic
matter degradation: the SoilR package, version 1.0. Geosci Model
Dev. 5, 1045–1060 (2012)
Wutzler, T., Reichstein, M.: Colimitation of decomposition by substrate and decomposers—a comparison of model formulations.
Biogeosciences. 5, 749–759 (2008)
Grant, R.F., Juma, N.G., McGill, W.B.: Simulation of carbon andnitrogen transformations in soil: microbial biomass and metabolic
products. Soil Biol Biochem. 25, 1331–1338 (1993)
Herbst, M., Hellebrand, H.J., Bauer, J., Huisman, J.A., Simunek, J.,
Weihermuller, L., Graf, A., Vanderborght, J., Vereecken, H.:
Multiyear heterotrophic soil respiration: evaluation of a coupled
CO2 transport and carbon turnover model. Ecol Model. 214, 271–
283 (2008)
Del Grosso, S.J., Parton, W.J., Mosier, A.R., Ojima, D.S., Kulmala,A.E., Phongpan, S.: General model for N2O and N2 gas emissions
from soils due to denitrification. Glob Biogeochem Cycles. 14,
1045–1060 (2000)
Parton, W.J., Holland, E.A., Del Grosso, S.J., Hartman, M.D.,
Martin, R.E., Mosier, A.R., Ojima, D.S., Schimel, D.S.:
Generalized model for NOx and N2O emissions from soils. J
Geophys Res. 106, 17403–17417 (2001)
Braakhekke, M.C., Beer, C., Hoosbeek, M.R., Reichstein, M.,
Kruijt, B., Schrumpf, M., Kabat, P.: SOMPROF: A vertically explicit soil organic matter model. Ecol Model. 222(10), 1712–1730(2011)
Fang, C., Moncrieff, J.B.: A model for soil CO2 production andtransport 1: model development. Agr Forest Meteorol 95, 225–
236 (1999)
Jassal, R.S., Black, T.A., Drewitt, G.B., Novak, M.D., GaumontGuay, D., Nesic, Z.: A model of the production and transport of
CO2 in soil: predicting soil CO2 concentrations and CO2 efflux
from a forest floor. Agric For Meteorol. 124, 219–236 (2004)Šimůnek, J., Suarez, D.L.: Modeling of carbon dioxide transportand production in soil, 1, Model development. Water Resour Res.
29(2), 487–497 (1993)
Pansu, M., Sarmiento, L., Rujano, M.A., Ablan, M., Acevedo, D.,
Bottner, P.: Modeling organic transformations by microorganismsof soils in six contrasting ecosystems: validation of the MOMOS
model. Glob Biogeochem Cycles 24 (2010).
Batlle-Aguilar, J., Brovelli, A., Porporato, A., Barry, D.A.:
Modelling soil carbon and nitrogen cycles during land use change.
A review. Agron Sustain Dev. 31(2), 251–274 (2011)
Fontaine, S., Barot, S.: Size and functional diversity of microbe
populations control plant persistence and long-term soil carbon accumulation. Ecol Lett. 8, 1075–1087 (2005)

Comput Geosci
24.

Thaysen, E.M., Jacques, D., Jessen, S., Andersen, C.E., Laloy, E.,
Ambus, P., Postma, D., Jakobsen, I.: Inorganic carbon fluxes acrossthe vadose zone of planted and unplanted soil mesocosms.
Biogeosciences. 11, 7179–7192 (2014)
25. Jacques, D., Šimůnek, J., Mallants, D., van Genuchten, M.T.:
Modeling coupled hydrologic and chemical processes: long-term
uranium transport following phosphorus fertilization. Vadose Zone
J. 7, 698–711 (2008)
26. Jacques, D., Šimůnek, J., Mallants, D., van Genuchten, M.T.: The
HPx software for multicomponent reactive transport during
variably-saturated flow: recent developments and applications. J
Hydrol Hydromech 66, 211–226 (2018)27. Steefel, C.I., Yabusaki, S.B., Mayer, K.U.: Reactive transportbenchmarks for subsurface environmental simulation. Comput
Geosci. 19, 439–443 (2015)
28. Carrayrou, J., Hoffmann, J., Knabner, P., Krautle, S., de
Dieuleveult, C., Erhel, J., van der Lee, J., Lagneau, V., Mayer,
K.U., Macquarrie, K.T.B.: Comparison of numerical methods for
simulating strongly nonlinear and heterogeneous reactive transport
problems-the MoMaS benchmark case. Comput Geosci. 14(3),483–502 (2010)
29. Mayer, K.U., Frind, E.O., Blowes, D.W.: Multicomponent reactive
transport modelling in variably saturated porous media using a generalized formulation for kinetically controlled reactions. Water
Resour Res. 37, 1174 (2002)
30. Bea, S.A., Wilson, S.A., Mayer, K.U., Dipple, G.M., Power, I.M.,
Gamazo, P.: reactive transport modeling of natural carbon sequestration in ultramafic mine tailings. Vadose Zone J. 11, 1–17 (2012)
31. Steefel, C.I., Appelo, C.A.J., Arora, B., Jacques, D., Kalbacher, T.,
Kolditz, O., Lagneau, V., Lichtner, P.C., Mayer, K.U., Meeussen,
J.C.L., Molins, S., Moulton, D., Shao, H., Simunek, J., Spycher, N.,
Yabusaki, S.B., Yeh, G.T.: Reactive transport codes for subsurface
environmental simulation. Comput Geosci. 19, 445–478 (2015)32.

Manzoni, S., Porporato, A.: Theoretical analysis of nonlinearities
and feedbacks in soil carbon and nitrogen cycles. Soil Biol
Biochem. 39, 1542–1556 (2007)
33. Moore, W. J.: Physical Chemistry, 4th Ed. Prentice-Hall, N. J
(1972).
34. Millington, R.J.: Gas diffusion in porous media. Science. 130, 100–
102 (1959)
35. van Genuchten, M.T.: Closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J. 44,
892–898 (1980)
36. Mayer, K.U., Alt-Epping, P., Jacques, D., Arora, B., Steefel, C.I.:
Benchmark problems for reactive transport modeling of the generation and attenuation of acid rock drainage. Comput Geosci. 19,
599–611 (2015)
37. Gérard, F., Blitz-Frayret, C., Hinsinger, P., Pagès, L.: Modelling theinteractions between root system architecture, root functions andreactive transport processes in soil. Plant Soil. 413, 161–180 (2017)38. Nowack, B., Mayer, K.U., Oswald, S.E., van Beinum, W., Appelo,
C.A.J., Jacques, D., Seuntjens, P., Gérard, F., Jaillard, B., Schnepf,
A., Roose, T.: Verification and intercomparison of reactive transport
codes to describe root-uptake. Plant Soil. 14, 305–321 (2006)
39. Leterme, B., Blanc, P., Jacques, D.: A reactive transport model for
mercury fate in soil—application to different anthropogenic pollution sources. Environ Sci Pollut Res. 21(21), 12279–12293 (2014)
40. Kuka, K., Franko, U., Ruhlmann, J.: Modelling the impact of porespace distribution on carbon turnover. Ecol Model. 208, 295–306(2007)
41. Keiluweit, M., Nico, P.S., Kleber, M., Fendorf, S.: Are oxygen
limitations under recognized regulators of organic carbon turnover
in upland soils? Biogeoch. 127, 157–171 (2016)
Publisher’s note Springer Nature remains neutral with regard tojurisdictional claims in published maps and institutional affiliations.