RESEARCH ARTICLEThe Semi-Analytical Solution for Non-Equilibrium Solute
Transport in Dual-Permeability Porous MediaKey Points:

Abhimanyu Sharma1
Jirka Šimůnek3

10.1029/2020WR029370• S emi-analytical solution for
non-equilibrium solute transport
in a dual-permeability porous
system with an immobile region ispresented• The global sensitivity analysis iscarried out for input parameters fora non-reactive solute• An insight into intra- and interaggregate mass transfers is presented
Correspondence to:
D. Swami,
dr.deepak.iitm@gmail.com
Citation:
Sharma, A., Swami, D., Joshi, N.,
Chandel, A., & Šimůnek, J. (2021).
The semi-analytical solution for
non-equilibrium solute transport
in dual-permeability porous
media. Water Resources Research,
57, e2020WR029370. https://doi.
org/10.1029/2020WR029370
Received 4 DEC 2020
Accepted 21 APR 2021

, Deepak Swami1, Nitin Joshi2, Aman Chandel1

, and

1

School of Engineering, Indian Institute of Technology Mandi, Mandi, India, 2Department of Civil Engineering, Indian
Institute of Technology Jammu, Jammu, India, 3Department of Environmental Sciences, University of California
Riverside, Riverside, CA, USA

Abstract This study presents a conceptual Dual-Permeability Non-Equilibrium (DPNE) model that
accounts for both physical and chemical non-equilibria to describe the reactive solute transport through a
porous medium. A semi-analytical solution of the DPNE model is derived in the Laplace domain, whichis then numerically inverted to obtain concentrations in different domains at different times and depths.
The derived semi-analytical solution is validated using experimental data and existing analytical (forsimple problems) and numerical solutions. The Global Sensitivity Analysis (GSA) is performed to identify
model parameters with an impact on nonreactive tracer breakthrough curves (BTC) for a selected data
set. A sensitivity-based calibration of the DPNE model’s parameters is carried out to simulate a BTC withmultiple inflection points. The DPNE model can better describe the experimental datasets than the model
based on the dual-porosity (MIM) concept with first-order mass transfer. Since the DPNE model considers
detailed sorption and mass transfer dynamics, it is suitable for process-based investigations.
1. Introduction
Transport of contaminants through porous media often experiences non-equilibrium conditions and remains a challenging problem for groundwater hydrologists. The porous media heterogeneity is associatedwith irregular stratification, dead-end pores, fractures, and lenses with different hydraulic conductivities(Bai & Roegiers, 1994, 1997; Bai et al., 1993; Kartha & Srivastava, 2008). It has been widely recognizedthat heterogeneities in porous media’s physical and chemical properties significantly influence transport
processes and result in an anomalous transport behavior (Brusseau et al., 1989; Brusseau & Rao, 1989;
Coats & Smith, 1964). The solute transport through such porous media is distinguished by asymmetry, early
breakthrough, and extensive tailing in breakthrough curves (BTC) (Brusseau, 1992; Brusseau et al., 1989).
Classical advection-dispersion equation-based models fail to capture the solute plume behavior under thenon-Fickian and nonequilibrium conditions. Brusseau and Srivastava (1997) reported several factors affecting the transport of organic solutes, including physical and chemical heterogeneity and nonlinear and
rate-limited sorption. Several conceptual models have been developed to analyze the fate of solutes in heterogeneous porous media by considering nonequilibrium, multi-rate mass transfer, and/or non-Fickian transport behavior (Brusseau et al., 1989; Haggerty & Gorelick, 1995; Pot et al., 2005; Selim et al., 1999; Sharma
et al., 2015, 2016; Srivastava & Brusseau, 1996; van Genuchten & Wierenga, 1976).
The nonequilibrium transport can be categorized as physical nonequilibrium (PNE) or chemical nonequilibrium (CNE). PNE is associated with transport-related non-equilibrium and affects both sorbing and nonsorbing solutes (Brusseau et al., 1989; van Genuchten et al., 1977; van Genuchten & Wierenga, 1976). PNE
is often modeled using dual- (multi-) porosity and/or dual- (multi-) permeability approaches (Šimůnek
et al., 2003; Šimůnek & van Genuchten, 2008).

© 2021. American Geophysical Union.
All Rights Reserved.

SHARMA ET AL.

The dual-porosity models divide the porous medium into two interacting continua (Coats & Smith, 1964;
Gaudet et al., 1977; Gerke & van Genuchten, 1993, 1996; Goltz & Roberts, 1986; Gwo et al., 1995; Selim
et al., 1999; Selim & Ma, 1998). While the first region is associated with macropores or fractures, the second region is associated with micropores or rock matrix. The dual-porosity models assume that advective flow is limited to the fracture domain and that dissolved solutes can move in and out of microporesor rock matrix due to molecular diffusion or first-order mass transfer (van Genuchten & Wierenga, 1976;
Swami et al., 2014, 2016, 2018). The most popular model based on the dual-porosity approach is the
1 of 23

Water Resources Research

10.1029/2020WR029370

mobile-immobile model (MIM) (Gaudet et al., 1977; Skopp et al., 1981; van Genuchten & Wierenga, 1976;
Warren & Root, 1963).
The dual-permeability models consider flow and transport in both macropore and micropore regions (Gerke& van Genuchten, 1993). These models assume that water flows relatively faster in the macropore regionand slower in the micropore region. The dual-permeability models can be further extended to account for
the immobile region in the matrix domain (e.g., Šimůnek & van Genuchten, 2008). The solute exchange
between the fracture and matrix regions and between the matrix and its immobile region occurs due tomolecular diffusion. Most studies considering the dual-permeability concept solve the governing equations
using numerical techniques (Gerke & van Genuchten, 1996; Zimmerman et al., 1993).
A few conceptual models incorporate CNE into a dual-permeability model (Kartha & Srivastava, 2008;Šimůnek & van Genuchten, 2008). CNE is often incorporated using two-site sorption conceptualization,
which divides the sorption sites into two fractions with equilibrium and kinetic sorption (van Genuchten& Wagenet, 1989). Kartha and Srivastava (2008) developed a dual-permeability model that considers porous domains with slow and fast-moving liquid and an immobile zone. Their model’s conceptualization
was analogous to the multi-process non-equilibrium (MPNE) model of Brusseau et al. (1989). The two-siteconceptualization represented the sorption process. They assumed that both equilibrium and kinetic sites
exist in the immobile and slow regions and that no sorption occurs in the fast region. A dual-permeability
model with an immobile zone inside of the matrix domain and two-site sorption in both domains was also
incorporated into the HYDRUS-1D software (Gerke & van Genuchten, 1993; Šimůnek et al., 2003, 2008)
and successfully applied to multiple experimental datasets (e.g., Dousset et al., 2007; Köhne et al., 2006; Potet al., 2005). Equilibrium and kinetic sorption can take place in both macropore and micropore regions in
HYDRUS-1D.
Several studies developed analytical or semi-analytical solutions for the dual-porosity model incorporating
PNE and CNE (Carnahan & Remer, 1984; Gao et al, 2009, 2010; Goltz & Roberts, 1986; Joshi et al., 2012; Leijet al., 1993; Neville et al., 2000; Toride et al., 1993; van Genuchten & Wagenet, 1976, 1989). However, only
a few studies developed analytical/semi-analytical solutions for the dual-permeability model (Leij & Bradford, 2013; Leij et al., 2012; Liang et al., 2019; Skopp et al., 1981). To the best of our knowledge, an analytical model considering a dual-permeability system with an immobile region and chemical non-equilibrium(referred to here as the DPNE model) has not yet been reported. The present study aims to fill this research
gap and develop a semi-analytical solution for the DPNE model. The developed semi-analytical solution isverified against the semi-analytical solution for the MIM and MPNE models (i.e., simpler problems) andthe numerical solution in HYDRUS-1D (i.e., a more complex problem). To identify the most sensitive DPNE
model parameters, the global sensitivity analysis was performed, and its results were then applied to analyze the existing experimental data set from the literature.

2. Conceptualization of the DPNE Model and Its Semi-Analytical Solution
The concept of the DPNE model with an immobile liquid zone in porous media and its implications for thesolute transport is mathematically described in this section. The water movement in the porous media iscategorized into fast, slow, and immobile liquid fractions. The fast and slow-moving liquids exist in separate
hydraulically connected pore regions having vfs and vsl velocities.
The conceptual approach of the DPNE model is similar to that of the multi-process non-equilibrium model
developed by Brusseau et al. (1989). The conceptual model considers advection, dispersion, first-order masstransfer, and sorption processes. The schematic of the conceptual model is presented in Figure 1. The soluteis assumed to coexist in the liquid and sorbed phases in the fast, slow, and immobile pore regions (zones).
The existence of three liquid domains results in three different liquid phase solute concentrations: Immobile (cim), slow (csl), and fast (cfs) domain concentrations (Figure 1).
These three solute concentrations (in different pore regions) have corresponding instantaneous andrate-limited sorbed concentrations since sorption is assumed to occur in all three zones. Additionally, thesolute mass transfer between the liquid and sorbed phases in the same porous region and between liquid
phases in different pore regions is considered. The mass transfer between sorbed phases in two different
SHARMA ET AL.

2 of 23

Water Resources Research

10.1029/2020WR029370

Figure 1. Conceptualization of the DPNE model. DPNE, dual-permeability non-equilibrium.

pore regions is not considered. The mass transfer between fast and slow-moving liquids is governed by the
first-order rate process, ɷsf (t−1). Similarly, as shown in Figure 1, the mass transfer between slow and immobile liquids is governed by the first-order rate process, ɷim (t−1). Instantaneous sorption in the three poreregions is defined by the three distribution coefficients Kim, Ksl, and Kfs. The mass transfer for rate-limited
sorption in the three pore regions is defined by kim, ksl, and kfs. The model assumes that the mass transfer
rates ɷsf and ɷim are constant and independent of time and space. The upper-case K represents distribution
coefficients for linear sorption, the lower case k represents rate-limited sorption, the symbol ω represents
mass transfer coefficients, and subscripts fs, sl, and im represent fast, slow, and immobile zones.

2.1. Governing Equation of DPNE
Three coupled one-dimensional equations represent the transport processes in porous media with different
flow rates in the fast, slow, and immobile regions. The processes of advection and dispersion in the longitudinal direction in both slow and fast regions and mass transfer are considered. The first-order mass transfer
is the coupling process between two different pore regions. Also, both instantaneous equilibrium sorption
and first-order rate-limited sorption occur in all pore regions, corresponding to the classical two-site sorption model of van Genuchten and Wierenga (1976).
The governing solute transport equation for the fast region can be given as:
S fs2
C fs
C fs
 2C fs
(1) q fs  fs D fs
 sf C fs  Csl  f fs b
 fs  f fs Ffs  b K fs
2
t
x
t
x









The solute transport equation for the slow region can be given as:
Ssl 2CCsl
 2Csl sf Csl  C fs  im  Csl  Cim   fsl b sl Dsl
(2)sl  fsl Fsl b Ksl  tsl  qsl
2
x
t
x





The mass balance equation for the immobile zone in the slow region can be given as:
S
C
(3)
im Csl  Cim im  Fim fim b Kim  tim  fim b imt 2 

The rate-limited sorption kinetics in fast, slow, and immobile regions can be represented as:
S fs2(4)k fs  1  Ffs K fsC fs  S fs2 


t



SHARMA ET AL.



3 of 23

Water Resources Research

10.1029/2020WR029370

Ssl2
(5)
ksl 1  Fsl  K slCsl  Ssl2 

t
Sim2
(6)
kim 1  Fim  KimCim  Sim2 trespectively. The kinetic sorption model described by Equations 4–6 allows for different sorption and desorption rates.
2.2. Initial and Boundary Conditions
A semi-infinite domain is considered for obtaining the analytical solution. Initially, it is assumed that themedia is free of contamination in all regions and phases.
C
C
C
S
S
S
0
fs  x ,0 
sl  x ,0 
im  x ,0 
fs  x ,0 
sl  x ,0 
im  x ,0 
(7)The inflow boundary condition for the fast and slow regions are:C fs
 0,t (8)v f C fs
  Df
 v fsC0 1   t  t0  0,t 
x





Csl
 0,t (9)vslCsl  Dsl vslC0 1   t  t0  0,t xwhere C0 is the source concentration, δ = 0 represents a constant concentration (Dirichlet) boundary condition, and δ = 1 represents a constant flux (Cauchy) boundary condition.
For the outflow boundary, a semi-infinite domain is considered. The outflow boundary condition for the
fast region is:
C fs  , t  0
(10)and for the slow region is:
Csl  , t  0(11)2.3. Solution in the Laplace Domain
The analytical solution is derived in the Laplace domain. A detailed derivation is given in Appendix A. The
numerical inversion of the Laplace domain solution is done using de Hoog’s algorithm (de Hoog et al., 1982).
de Hoog’s algorithm is used extensively for numerous solute transport problems (Joshi et al., 2012; Nevilleet al., 2000; Sudicky & Frind, 1982; Tang et al., 1981). This algorithm uses Fourier series approximation
based on the quotient difference algorithm for the approximation of the inverse Laplace transform and isgiven as:C2M
fs  x , p 1
exp  pt  Re   C fs ik  
C fs  x ,t 
(12)
 2
T
k 1  x , p 

2 



where T defines the period of approximating Fourier series and T = 0.80t (Neville et al., 2000). C fs represents the fast region concentration in the Laplace domain.
ln  Er (13)s  
2T

SHARMA ET AL.

4 of 23

Water Resources Research
Table 1
Values of Different DPNE Model Parameters Used to Compare the
Proposed Model With Other (MIM, MPNE, and HYDRUS-1d) Models

No.

Parameter

Units

MIM

MPNE

HYDRUS(without
immobile
water)

1

L

cm

30

30

30

30

2

ρ

g/cm3

1.30

1.36

1.25

1.25

3

D

cm2/day

30

–

4

α

cm

–

0.34

5.74

0.03

5

θfs

–

0.26

0.439

0.082

0.002

6

θsl

–

0.14

0.034

0.316

0.063

7

θim

–

0

0

0

0.333

8

qfs

cm/day

10

5.11

6.19

0.771

9

qsl

cm/day

0

0

1.20

6.62

10

ωsf

day−1

0.15

0.075

1.063

0.072

11

ωim

day−1

0

0

0

0.470

12

ffs

–

0.4

0.4

0

0

13

fsl

–

0.6

0.6

0

0

14

f im

–

0

0

0

0

15

Ffs

–

1

0.5

0

0

16

Fsl

–

1

0.5

0

0

17

Fim

–

1

0

0

0

18

kfs

day−1

0

0.663

0

0

19

ksl

day−1

0

0.663

0

0

20

kim

day−1

0

0

0

0

21

Kfs

cm3/g

0.5

0.429

0

0

22

Ksl

cm3/g

0.5

0.416

0

0

Kim

3

0

0

0

0

23

cm /g

HYDRUS
(with
immobile
water)–

10.1029/2020WR029370

de Hoog et al. (1982) suggested the following values of constants for accelerating the convergence of the Fourier series: ϵ = 0, Er = 0.001, andM = 7.

3. Verification of the Developed Solution
This section demonstrates the DPNE model’s capabilities in reproducing BTCs obtained by a range of models based on the dual-porosity anddual-permeability concepts. Two analytical solutions from the literature(Neville et al., 2000; van Genuchten & Wierenga, 1976) are used to verify
the DPNE model for the dual-porosity system. For this verification, theDPNE model is converted/simplified into the MIM and MPNE models
by selecting specific values of various coefficients. The numerical HYDRUS-1D model is used to fully verify the DPNE model for the dual-permeability system with or without an additional immobile region.
The analytical solution for the MIM model considered here is limited to
instantaneous equilibrium sorption. Therefore, the governing equationsof the DPNE model are simplified by taking Ffs = Fsl = Fim = 1. WhenFfs = Fsl = Fim = 0, sorption in the fast, slow, and immobile porous regions
is reduced to first-order kinetics. Only physical non-equilibrium without
sorption is considered for the dual-permeability system and the comparison with HYDRUS-1D. Therefore, the governing equations of the DPNE
model are simplified by taking Kfs = Ksl = Kim = kfs = ksl = kim = 0. Table 1
shows a list of input parameters for all three modeling setups.
3.1. Comparison With the MIM Solution Developed by van
Genuchten and Wierenga (1976)
van Genuchten and Wierenga (1976) developed an analytical solution
for a physical non-equilibrium model considering a bi-continuum liquid
phase concentration, that is, a dual-porosity model. The concentrationdifference between two hydraulically connected regions governs thefirst-order mass transfer between them. The DPNE model is converted
to the MIM model by considering ωim = 0 in the governing equations of
the DPNE model (shown in Appendix C). In the dual-porosity approach,the mobile and immobile pore regions are hydraulically connected. The
fast region processes (in the DPNE model) are analogous to the mobile
region processes (in the MIM model), while the slow region processes are analogous to the immobile region processes. This is achieved byconsidering qsl = 0 and the reaction parameters Ffs = Fsl = Fim = 1 and
kfs = ksl = kim = 0 as the rate-limited sorption is not considered. The listof additional parameters is listed in the MIM column of Table 1. Figure 2
shows that the results of the downscaled DPNE model and the MIM
model are identical.
3.2. Comparison With the Semi-Analytical Solution for the MultiProcess Non-Equilibrium Model (MPNE)

Figure 2. Verification of the downscaled DPNE model with the MIM(Mobile-Immobile) model. Experimental data are from van Genuchten
and Wierenga (1976). DPNE, dual-permeability non-equilibrium; MIM,
mobile-immobile model.

SHARMA ET AL.

The MPNE model is an extension of the bi-continuum model by including site-specific sorption. The MPNE model thus considers both physicaland chemical non-equilibria. A semi-analytical solution for the MPNE
model incorporating instantaneous and rate-limited sorption was presented by Neville et al. (2000). The semi-analytical solution of the DPNE
model can be reduced to the MPNE model by switching off the immobile
5 of 23

Water Resources Research

10.1029/2020WR029370

region by substituting ωim = 0 and qsl = 0, similarly as in the previous
section. To convert the DPNE model to the MPNE model, additional
instantaneous and rate-limited sorption parameters were considered inboth mobile and immobile domains. A list of input parameters is shownin Table 1. Figure 3 shows an excellent agreement between the DPNE and
MPNE models.

3.3. Comparison With the Numerical Solution for the DualPermeability System With Immobile Water in the Matrix
The DPNE model is further validated in this section using the numerical solution for solute transport in a dual-permeability porous medium
incorporated in the HYDRUS-1D model (Šimůnek et al., 2003, 2008). A
soil column experiment involving the transport of the CaCl2 solution conFigure 3. Verification of the downscaled DPNE model with the MPNE
taining Br (Pot et al., 2005) is used for the validation. The undisturbed soilmodel. Experimental data are from Neville et al. (2000). DPNE, dualcolumn 30 cm long was collected from a soil layer with hydromorphic silt
permeability non-equilibrium; MPNE, multi-process non-equilibrium.
loam (with 24% clay, 40% silt, and 36% sand). The experiment was conducted under unsaturated, steady-state flow conditions. The effluent wascollected at 0.5 h intervals. More details can be found in Pot et al. (2005),
who analyzed the collected BTC using the dual-permeability model with or without an immobile region in
the soil matrix implemented in HYDRUS-1D. The list of parameters used to validate the DPNE model thatconsiders the dual-permeability system with and without an immobile region in the soil matrix is shown in
Table 1 (the last column). Figure 4 shows an excellent agreement between the results of the DPNE modeland those simulated by the HYDRUS-1D model.

4. Global Sensitivity Analysis
Being a complex, process-based model, the DPNE model requires a large number of input parameters. In
its full complexity (i.e., multiple porosities, multiple processes), many of these model parameters may not
be identifiable using typically collected information (i.e., outflow breakthrough curves). Under such conditions, various assumptions about different parameters need to be considered. For example, some parameters
can be eliminated to simplify the model by neglecting one or two pore domains or considering only instantaneous or kinetic sorption (as shown in the examples above). For less sensitive parameters, typical values
may be obtained from the literature. Similar parameters can be linked,
such as assuming the same sorption and/or degradation parameters in
different phases or regions. Multiple parameter optimization steps can beused to estimate various parameters. For example, Pot et al. (2005) firstused the non-reactive tracer BTCs to estimate flow parameters (fluxes andmass transfer between different regions). The reactive tracer BTCs werethen used to estimate reaction parameters (for equilibrium and kinetic
sorption). Such strategies and assumptions will depend on the system’s
complexity and the information content of collected data.

Figure 4. Verification of the DPNE model with the Hydrus-1D model (thedual-permeability model with (DP_IM) and without (DP) immobile waterin the matrix). Experimental data are from Pot et al. (2005). DPNE, dualpermeability non-equilibrium.

SHARMA ET AL.

In the case of a process-based model, the simulated output may not be
similarly sensitive to all input parameters simultaneously. Before estimating or calibrating the model for a particular system (e.g., MIM, MPNE,
DPNE) against the field/laboratory data, it is important to identify themost sensitive parameters, that is, those that have the largest impact on
the results and are thus pivotal for the model deployment, especially incase of a complex modeling approach (Ceriotti et al., 2019). This can beaccomplished by performing the Global Sensitivity Analysis (GSA). Such
practice outlines a basic framework for the model-based interpretation to
design the experimental work (Fajraoui et al., 2011).

6 of 23

Water Resources Research
Table 2
Parameter Values Used in the Sensitivity Analysis
Parameters
No.

Parameter

Units

Lower limit

Mean

Upper limit

1

qfs(cm/min)

0.19

0.28

0.37

2

qsl

(cm/min)

0.26

0.38

0.49

3

ωsf

(min−1)

1.7E−3

2.4E−3

3.2E−34

ωim(min−1)

3.6E−55.2E−5

6.8E−5

5

α

(cm)

0.41

0.58

0.92

10.1029/2020WR029370

Different sensitivity analysis (SA) approaches have been proposed for various research areas related to groundwater and environmental sciences(e.g., Ciriello et al., 2013; Dell'Oca et al., 2020; Hamby, 1994; Norton, 2008;
Pianosi et al., 2016), hydrogeological models (Wainwright et al., 2014),
complex kinetic systems (Turanyi, 1990), building energy analysis(Tian, 2013), and radioactive waste disposal models (Helton, 1993). For
a complex model, an increase in the number of input parameters is imperative. It is important to opt for an appropriate sensitivity analysis (SA)approach, specific to the type of the model input-output configuration(Nguyen & Reiter, 2015).

The Fourier Amplitude Sensitivity Test (FAST) is a variance-based method to quantify the effect of different parameters on the output of the
non-linear models (Cukier et al., 1973, 1975). The FAST model decomposes the model output variance into fractions within a Fourier amplitude integration framework and attributes the model output variance to particular parameters (assumed to be uncorrelated) or their subsets.
The details of the FAST formulation and calculation details are given in Appendix B.
In this section, the global sensitivity analysis for the DPNE model is carried out to determine the most sensitive model parameters for a particular selected data set so that these parameters can be estimated effectively.
The sensitivity analysis is based on Miyamoto et al.’s (2003) data set involving conservative solute transport
experiments on the Andisol soil column. Andisols, developed from volcanic ash, constitute 17% of Japan’sland and are commonly used for agriculture. Due to unique physical and chemical properties, Andisol’s
flow and transport properties are of prime interest. Andisol has a distinct interaggregate and intra-aggregate
pore space arrangement due to its composition of non-crystalline minerals. The Andisol experiment was
conducted on a 26 cm long soil column with an internal diameter of 4.4 cm. An aqueous solution of CaCl2
of 0.06M (considered to be a non-reactive solute) was eluded at a steady-state saturated flow rate of 0.56 cm/
min. Four electrical conductivity probe sensors were inserted horizontally at 7.72 cm, 9.96 cm, 14.05 cm,
and 16.87 cm along the column to monitor the electrical conductivity (Toride et al., 2003) and experimental
breakthrough curves. All depths were considered for investigating the sensitivity of parameters. Please refer
to Section 4 of Leij et al. (2012) for further details about the experiment.
Table 2 lists the DPNE model parameters, their means, and lower and upper limits considered in the global
sensitivity analysis. The mean values of these parameters were obtained by preliminary fitting of BTCs at
all four depths simultaneously. The lower and upper limits of these parameters were selected to ensurethat wide intervals in the parameter space were explored and to minimize the chance of sampling negativevalues.
The quantitative results illustrated in this section show the general relative importance of the input parameters and that model output sensitivity varies spatially and temporally for different input transport parameters. The results may vary based on the flow characteristics and transport properties.
Figure 5 compares the breakthrough curve profiles reported by Miyamoto et al. (2003) and related total sensitivity indices of the DPNE model parameters. Figures 5 A-D show relative sensitivities for five unknown
parameters qs, qf, ɷsf, ωim, and α. Figures 5 E-H then show relative sensitivities of three parameters ɷsf, ωim,and α, when the seepage velocities in the slow and fast flow regions are assumed to be known.
Figures 5A–5D show the FAST sensitivity indices for five DPNE parameters (qs, qf, ɷsf, ωim, and α) at distances of 7.72, 9.96, 14.05, and 16.87 cm. In all cases, the sensitivity indices for qs and qf are significantly larger
than those for the other three parameters. It is also interesting to point out that the sensitivity indices for qs
and qf show an inverse correlation, identifiable by local peaks. Figure 5 also indicates that the dispersivity α
has negligible sensitivity during the initial and late stages of the BTC.
Figures 5E–5H show the FAST sensitivity indices for three less sensitive parameters (ɷsf, ωim, and α), assuming that water fluxes are known. In Figures 5E and 5F, the sensitivity is governed by α, revealing the importance of dispersivity at earlier times. At farther distances (Figures 5G and 5H), ɷsf dominates the sensitivity
index. The sensitivity index for ɷsf attains sharp peaks at the BTC’s inflection points due to the impact of

SHARMA ET AL.

7 of 23

Water Resources Research

10.1029/2020WR029370

Figure 5. Breakthrough curves (Miyamoto et al., 2003) (dashed line) and Fast Sensitivity Indices (solid lines) associated with parameters optimized using
BTCs from all depths (AE – 7.72 cm, BF – 9.96 cm, CG – 14.05 cm, and DH – 16.87 cm). (A–D) Relative sensitivities of qs, qf, ɷsf, ωsim, and α. (E–H) Relative
sensitivities of ɷsf, ωsim, and α when the seepage velocities in the slow and fast regions are assumed to be known.

diffusive processes between primary and secondary porosities. Figures 5E–5H also illustrate that the model
output is sensitive to the parameter ɷim during later stages of the BTC.
Figures 5E–5H reveal important information about ɷsf, ωim, and α. Figures 5E–5H clearly show that at all
distances, overall sensitivity is mainly affected by α. Both parameters ɷsf and α show independent behavior,and the sensitivity of BTC at later stages is governed by both ωim and α. A careful observation of Figures 5E–
5H reveals temporal locations of inflection points and sharp peaks and drops in the FAST index of ɷsf and α.
At a short distance (7.72 cm), the sensitivity index of ɷsf has a single point of inflection. It has two at middle
distances (9.96 and 14.05 cm) and three at a large distance (16.9 cm). The sensitivity associated with ωim is
reduced at larger distances. The findings of Figure 5 imply that all five parameters, that is, qs, qf, ɷsf, ωim, andα are relevant. When the BTC has multiple peaks/inflections, these parameters contain essential information about the model response variance for model calibration.

5. Simulation of Experimental Data
The experimental data set for Andisol (Miyamoto et al., 2003), considered above in the global sensitivity
analysis, is analyzed in this section to evaluate the applicability of the developed model. The calibrationof the DPNE model parameters for BTCs collected at four observed locations is carried out based on GSA’s
results discussed above. A set of sensitive parameters for a nonreactive tracer was identified in the previous
section using the sensitivity analysis. Breakthrough curves at all observation points were assigned equal
weights.

SHARMA ET AL.

8 of 23

Water Resources Research

10.1029/2020WR029370

Table 3
Parameters Optimized for the DPNE and MIM Models Using PEST From Experimental Data

DPNE
qfs = 0.27 cm/mina
qsl = 0.29 cm/minaMIM
q = 0.56 cm/mina

Parameter

Mean and 95% confidence intervalα (cm)

0.61[1.08
0.20

ωsf (min−1)

E 3
2.51E−3[ 4.4
7.38 E  4

ωim (min−1)

E 5
3.1E-5[5.57.82 E  6

θfs

0.1[0.18
0.05

θsl

0.36[0.48
0.15

θim

0.11 [0.16
0.09

θm

0.21[0.32
0.12

θim

0.36[0.51
0.22

α (cm)

0.81[1.34
0.33

ω (min−1)

8.89 E 1
0.45[1.5E
2aFlux is fixed to minimize the number of optimized parameters.

The observed data were analyzed using the DPNE model, and its results were compared with those of thedual-porosity (MIM) model. The MIM and DPNE transport parameters were estimated using PEST (Doherty & Christensen, 2011) and are shown in Table 3. Since Darcy’s discharge (q) and total porosity wereexperimentally measured (Leij et al., 2012), they were fixed at 0.56 cm/min and 0.57, respectively. The dispersivity (α), mass transfer coefficient (ω), mobile (θm) and immobile (θim) water contents were optimized
for the MIM model. For the DPNE model, six parameters were optimized (θsl, θfs, θim, ɷsf, ɷim, and α). The
estimation was carried out with two constraints, that is, the total flux q = qsl + qfs and the total porosityθ = θsl + θfs + θim, which were defined in the objective function in PEST. The estimated parameters areshown in Table 3.
The DPNE model parameters acquired by simultaneously optimizing BTCs at four depths are shown in
Table 3. Figure 6 demonstrates that the DPNE model performs better than the MIM model, specificallyfor BTCs with multiple inflections at greater depths. It appears that the contaminant moves through three
hydraulically connected domains more or less independently, which is evident from small values of masstransfer coefficients. It can also be observed that the volume fractions of the fast and slow regions are equal.
Most of the solute transport is governed by the slow region, which is apparent from the estimated θsl. Masstransfer between the fast and slow regions is faster than between the slow and immobile regions, as shownby the corresponding mass transfer coefficients' values. It is important to highlight the MIM model’s limitations to capture the complex behavior of BTCs affected by multiregional transport phenomena. The optimized values of dispersivity and mass transfers coefficient are higher for MIM compared to DPNE.
The two models’ performance and suitability are compared in Table 4 to identify whether a simpler MIM
model can better describe the experimental data than a more complex DPNE model while using fewer
computational resources and fewer parameters. This is assessed by calculating the root mean square error (RMSE), determination coefficient (R2), Nash-Sutcliffe efficiency (NSE), Akaike information criterion(AIC), and Modified Akaike information criterion (AICc) using the following set of equations:
2
1 n
RMSE    Ciobs  Cisim (14)n i 1SHARMA ET AL.

9 of 23

Water Resources Research

10.1029/2020WR029370

Figure 6. Experimental data (Miyamoto et al., 2003) and their fits at distances of 7.72, 9.96, 14.05, and 16.87 cm using
the DPNE and MIM models with parameters optimized using data from all depths. DPNE, dual-permeability nonequilibrium; MIM, mobile-immobile model.





2


 
 
(15)
R 
2
2
 n C iobs  Ciobs  in 1 Cisim  Cisim  
i 1


n
Ciobs  Ciobs in 1 Cisim  Cisim
 i 1
2in 1  Ciobs  Cisim 
NSE  1 (16)2in 1 Ciobs  Ciobs
2





 


AIC N . Ln  2  2
(17)2  1
AICc N . Ln  2  2 (18) n    1 

Table 4The Goodness of Fit for Different (MIM and DPNE) Models to Simulate Experimental Data

Distance(cm)

SHARMA ET AL.

R2
DPNE

RMSE
MIM

DPNE

MIM

NSE
DPNE

MIM

AIC
DPNE

AICc
MIM

DPNE

MIM

Pearsoncoefficient
DPNE

MIM

7.72

0.993

0.998

0.026

0.016

0.993

0.997

−317.6−360.51−315.6−358.5

0.997

0.999

9.96

0.997

0.995

0.02

0.03

0.996

0.992

−338.8−309.5−336.8−307.5

0.999

0.981

14.05

0.994

0.96

0.038

0.051

0.988

0.969

−282.5−256.2−280.5−254.2

0.997

0.97

16.87

0.996

0.94

0.023

0.091

0.996

0.94

−326.7−208.3−324.8−206.3

0.998

0.95

10 of 23

Water Resources Research

10.1029/2020WR029370

where Ciobs and Cisim are ith observed and simulated concentrations, respectively, Ciobs and Cisim are the means of observed and simulated concentrations, respectively, n is the number of observations taken at any
particular observation point, N is the sample size,  is the number ofparameters in the model plus 1, and σ2 is the residual sum of squares/N.
The DPNE model, which also considers transfer processes in the immobile zone, can better describe observed BTCs than the MIM model when
the dual-porosity concept is used. To evaluate whether it is beneficial touse a more complex model, statistical measures such as AIC and AICc are
shown in Table 4. It is concluded that the DPNE model outperforms the
MIM model at all distances in simulating BTCs with multiple inflection
points. AIC and AICc also show the maximum likelihood criteria are infavor of DPNE.

Figure 7. Comparison of the zeroth temporal moments for experimental
data (Miyamoto et al., 2003) (Exp) and the DPNE and MIM model fits.
DPNE, dual-permeability non-equilibrium; MIM, mobile-immobile model.

6. Temporal Moments

This section compares temporal moments for the MIM and DPNE models derived from observed BTCs. Temporal moments can describe the
distribution of solute mass within a porous medium and have provenuseful in analyzing the solute mass response to various transport processes. The zeroth temporal moment gives the mass of solute recovered at the sampling location. The first(T1) and second (T11) temporal moments normalized by a zeroth moment give the mean arrival time andvariance of BTCs, respectively. The numerical expressions for temporal moments are as follows (Kendall &Stuart, 1977):T(19)M n  x    t nC x ,t  dt
0M
T1  x   1
(20)
Mo
2

M2  M1 (21)T11 xMo  Mo 

where Mn(x) represents the nth absolute temporal moment, T1(x) is thenormalized first temporal moment (by zeroth temporal moment), andT11(x) is the normalized second temporal moment or variance of the BTC.
Figures 7–9 show the moment analysis for the models compared. The solute mass in the advective domain at various distances is shown using thezeroth temporal moment in Figure 7. Note that the solute mass decreases
with distance as some solute is transferred from the advective domain to
the immobile domain. The mass accumulated at all distances is predicted precisely by the DPNE model. The estimation of mass at 14.05 and9.96 cm is more precise by MIM compared to DPNE. At other distances,
MIM overestimated the mass. The zeroth temporal moments show thatthe DPNE model, due to its consideration of detailed transport processes,
performs best for the BTC with complex behavior, which is observed at
16.87 cm. At other distances, the estimation of mass by the MIM and
DPNE models is comparable.

Figure 8. Comparison of the first temporal moments for experimental
data (Miyamoto et al., 2003) (Exp) and the DPNE and MIM model fits.
DPNE, dual-permeability non-equilibrium; MIM, mobile-immobile model.

SHARMA ET AL.

The residence time of the solute mass within the system and its mean
arrival time at different depths capture various transport processes occurring in the porous medium. The mean arrival time of reactive solutes
at any given observation location will be higher than for tracers. For
11 of 23

Water Resources Research

10.1029/2020WR029370

nonreactive solutes, such as the one considered here, the mass transfer
processes between different pore regions will determine the mean arrival
time. Shorter mean arrival times are expected for soils with multiple porosities and limited mass exchange. Such a phenomenon can be observed
using the first temporal moment. The MIM model slightly overestimatesmean arrival time at longer distances, whereas the DPNE model estimates the exact mean arrival time. This is due to the inclusion of physicaltransport processes and secondary and immobile pore regions. Hydraulically connected fast-slow and slow-immobile regions actively participate
in the mass transfer processes, which causes a delayed mean arrival time.

Figure 9. Comparison of the second temporal moments for experimental
data (Miyamoto et al., 2003) (Exp) and the DPNE and MIM model fits.
DPNE, dual-permeability non-equilibrium; MIM, mobile-immobile model.

Preferential flow through a complex porous medium causes larger spatial spreading of the solute plume. This spreading will be higher if PNEprocesses such as solute mass transfer between domains are involved.
For nonreactive solutes, these processes are mostly diffusion-based andinclude intra- and inter-domain mass transfers and intra-particle diffusion. For reactive solutes, in addition to diffusion-based processes, there
are chemical transformation reactions that will also cause larger spatial
spreading.

Better estimation of variance is done using DPNE compared to MIM, as
shown in Figure 9. Both DPNE and MIM models estimate higher variances at shorter distances and lower variances at larger distances. A comparison of the MIM and DPNE
parameters in Table 3 shows the overestimation of the MIM parameters, which resulted in higher varianceestimation. The lower variance estimated by the DPNE model at 16.87 cm is due to the higher liquid fractionof the slow region than the fast region (Figure 6).
From the analysis of temporal moments, it is evident that the DPNE model outperforms the MIM model.
The inclusion of additional physically based parameters is required to fully describe experimental BTCs,
which cannot be accomplished by the MIM model, which neglects relevant transport processes.

7. Conclusions
A semi-analytical solution for solute transport with non-equilibrium physical and chemical processes understeady-state water flow conditions in a dual-permeability porous medium is presented in this manuscript.
The semi-analytical solution was derived in the Laplace domain and numerically inverted in the time domain. The semi-analytical solution was validated using exiting analytical (MIM and MPNE) solutions for
simpler conditions and numerical solutions (HYDRUS-1D) for more complex conditions. Excellent agreements between the new semi-analytical solution and existing analytical and numerical solutions were obtained in all evaluated cases.
The global sensitivity analysis was used to identify model parameters, which the model output is the mostsensitive to, using an experimental data set of Miyamoto et al. (2003). The results of the sensitivity analysis
were then used in model calibration. The fluxes in slow and fast regions, that is, qs and qf, were identified as
the most important parameters during the entire experiment. On the other hand, the ɷsf and α parameters
were found to affect the BTC mainly during the initial breakthrough and at inflection points, respectively. It
was concluded that all five DPNE parameters played an important role when analyzing BTCs with multiple
inflection points.
The DPNE and MIM model parameters were estimated for an experimental data set of Miyamoto et al. (2003)
involving BTCs with multiple inflection points. The DPNE model outperformed the MIM model. The temporal moment analysis illustrated that the DPNE model could better capture the mass, mean arrival time,and spreading of the BTC than the DPM model.
Since the DPNE model incorporates two distinct flow regions and considers both physical and chemical
nonequilibria, it can be readily applied to analyze experimental data collected at the field scale. Additionally, the semi-analytical solution presented here can be used to verify complex numerical models considering
SHARMA ET AL.

12 of 23

Water Resources Research

10.1029/2020WR029370

multiple overlapping porous flow regions with simultaneous physical and chemical nonequilibrium
processes.

Appendix A: Analytical Solution of the DPNE Model in the Laplace Domain
Taking the Laplace transform of Equation 6 and solving it leads to:
pSim2 kim 1  Fim  KimCim  Sim2 (A1)where p is the Laplace operator.
Further solving Equation A1 we get:kim 1  Fim  KimCim
Sim2 (A2) p  kim 

Taking the Laplace transform of Equation 3 and solving it results in:

Rim pCim  fim b pSim2 
im Csl  Cim
(A3)

R
im  Fim fim b Kim
(A4)im
Cim  A1Csl
(A5)

where



im
A1 
(A6)
pRim  J1  imfim b pkim 1  Fim  K im
J1 
(A7)
 p  kim 

Taking the Laplace transform of Equation 5 and solving it leads to: ksl 1  Fsl  K sl  Csl
Ssl2  (A8)p  ksl

Taking the Laplace transform of Equation 2 and solving it: 2Csl
Csl
qsl
 qsl
 A2Csl  sf C fs 
0
(A9)
2
x
x

where:
A2  Rsl p  sf  im  im A1  J 2
(A10)

R sl  fsl b Fsl K sl
(A11)
sl
fsl b pksl 1  Fsl  K sl
(A12)J2 pSolving Equation A9, we get:sf
d 2Csl dCsl
A2


Csl 
C fs

(A13)
dx
qsl
qsl
dx 2
SHARMA ET AL.

13 of 23

Water Resources Research

10.1029/2020WR029370

Similarly, the solution for the fast region can be obtained as:2

d C fs dC fs
sf
A3



C fs Csl
(A14)
dx
q fs
qsl
dx 2

where
A3  R fs p  sf  J3
(A15)
R
 fs  f fs b Ffs K fs
(A16)fs



f fs b pk fs 1  Ffs K fs
J3 
(A17)p  k fsNow, Equations A13 and A14 are solved simultaneously using Eigenvalue and Eigenvector as shown:
L C  A.C
(A18)

where operator L is:d2d(A19)

L 

dx dx
 A2
sf 



q
qsl 
A   sl(A20)
 sf
A3 


q
q
fs
 fs


As the transformed concentration vector C is Csl

C fs  , Eigenvalues (λ1 and λ2) and Eigenvectors ( 1and  2)corresponding to matrix A are obtained using the Mathematica software. They are given as: A  A5  A6 
1   4
(A21)

2 q fs .qsl


 A  A5  A6 
2   4
(A22)

2 q fs .qsl
  A  A  A 
5
6  11     7



v1 
2 sf qsl
(A23)  
 21  1

  A  A  A 
5
6  12     7



v2 
2sf qsl
(A24)  
 22  
1



where
SHARMA ET AL.

14 of 23

Water Resources Research

10.1029/2020WR029370


A4 A2 q fs  A3qsl
(A25)


A5 A22 q 2fs  2 A2 A3q fs qsl
(A26)

A6 4 sf2 q fs qsl  A32 qsl2
(A27)

A7 A2 q fs  A3qsl
(A28)

There exists a modal matrix P, using which matrix A is diagonalizable because matrix A will have distinct
roots. These roots can be calculated as:
(A29)A  I 
0
Modal matrix P is defined from the Eigen vectors as: 11  12 P(A30)
 21  22 

Also, the modal matrix P has a property that allows the diagonalization of matrix A: 1 0 
P 1 AP  (A31)
0 2 
Now, a setting is done:C  PY(A32)

Substituting C in Equation A18:
L  PY   A  PY 
(A33)As P is linear:
PL Y   APY
(A34)

Pre-multiplying Equation A34 with P−1 leads to:
L Y   P 1 APY
(A35)

Using Equation A19 and A31, Equation A35 becomes:
 d 2Y1dY   2  1   0  Y 
1dx   1 dx
(A36) 
 d 2Y2dY2   0 2  Y2 



dx 
 dx 2

Solving the above equation for Y1, we get:d 2YdY 21  1 
1Y1
(A37)dxdx D 2  D  1 Y1 
0
(A38)So,
SHARMA ET AL.

15 of 23

Water Resources Research

10.1029/2020WR029370

1  1  41
(A39)
D1,2 
2
1  1  42
(A40)
D3,4 
2
x

x

DD
(A41)
Y1 m1e 1  m2e 2

and
x

x

D
D
(A42)
Y2 m3e 3  m4 e 4For solving Equation A41, we need the boundary condition in terms of Y. Since C  PY : 11  12  Y1 C 
(A43) 
 21  22  Y2 

At the boundary condition corresponding to C, from Equations 7–11, it follows:C 0,t   C0
(A44)C
0
(A45) , t  xTherefore, at the boundary condition corresponding to C :CC 0, p   0
(A46)p
C
0
(A47) , p  xWhen the boundary conditions are transformed corresponding to Y, since C = PY, we get:(A48)Y  P 1CNow
adj PP 1 
(A49)P

Using the values of P from Equation A30, we get:  22  12 P 1    11 22   12 21 
(A50)
  21  11 
Substituting P−1 in Equation A48:
Y1   22
 12  Csl      11 22   12 21 
(A51)
Y2    21  11  C fs 
gives a value of:
SHARMA ET AL.

16 of 23

Water Resources Research

10.1029/2020WR029370

 22Csl   12C fs
Y1 (A52)
 11 22   12 21
and 11C fs   21Csl
Y2 
(A53)
 11 22   12 21
Dx

The outlet conditions suggest that m1 = m3 = 0 and thus Y1  m2e 2 .
At x = 0:
  22   12
 C0
m
Y
(A54)
.
21
  11 22   12 21  p

Therefore: D2x
C   22   12
Y1  0 
(A55).ep   11 22   12 21 

and
 D4x
C   11   21
Y2  0 
(A56).ep   11 22   12 21 

Using C=PY:
Csl   11  12  Y1  (A57)  22  Y2 C fs   21Csl  11Y1   12Y2
(A58)
C fs  22Y2   21Y1
(A59)

Appendix B: Details About the FAST Sensitivity Index
The expected value of the output E{y} for a function f (x1, x2,…, xm) containing x1, x2,…, xm parameters can beevaluated by:E  y  ...  f  x1, x2 ,..., xm  P  x1, x2 ,..., xm  dx1, dx2 ,..., dxm
(B1)where P is the m-dimensional probability density for x. The FAST method expands the model into a Fourier
series where Fourier coefficients are used to calculate the mean and the variance of the model.
f  X     ...  Ck1 .k2 ...km e  1 1 2 2 m m (B2)k1   k2  

j 2 k x  k x ... k x

km  

with 1 1 2 2 m m  dX
Ck1k2 ...km   f  X  . e
(B3) j 2 k x  k x ... k x

Im

SHARMA ET AL.

17 of 23

Water Resources Research

10.1029/2020WR029370

Using Equation 74, component fi1 …is (xi1 …, xis) of the analysis of variance (ANOVA), decomposition canbe presented as Fourier series having ki1…kis as the only non-null indices (Sobol’, 1990, 1993). Applying thismethod, the variance can be expressed as the sum of modulus of the Fourier coefficients.









k1  

kis  var fi1 ...is  ...  Cki ...ki
(B4)
s
1

To find the Fourier coefficients, the multi-dimensional integral in Equation B3 needs to be solved. Ergodic
theorem is used to convert this equation into a one-dimensional integral by expressing every parameter as afunction of a new independent variable ‘s' (Saltelli et al., 1999).
1xi  s(B5) 2  arcsin sin i s Here, ωi is the frequency assigned to xi, and s is the search variable. As Equation B5 is periodic in nature, svarying between −∞ and ∞ traverses the entire curve (McRae et al., 1982). Hence, the expected output value
is calculated as
1 
E 
y
 f  s  . ds
(B6)
2 

Similarly, the variance of the first-order function, which only depends on the input parameter, can be represented by Fourier coefficients.
E  y | xi    Cki
(B7)iwhere coefficients can be calculated as:
1 C f  s  .e  j 2 ki i s ds
(B8)
ki
2 

Results from Equations B6 and B7 are used to calculate the first-order sensitivity index, denoted as Si, whichmeasures the effect of input parameter xi on the output variance, ignoring the interaction with other parameters (Sobol', 1993). It is estimated to be:

var E  y | xi 
(B9)
Si 
var  yThe value of Si varies from 0 to 1. An increase in this value represents an increase in the corresponding
parameter’s importance in the model outcome.

Appendix C: Comparison of the DPNE Model With Simpler Models
Comparison with the multi-process non-equilibrium model (MPNE) presented by Nevilleet al. (2000)
We start with the governing equations in the DPNE model. The governing solute transport equations for thefast, slow, and immobile regions are given as follows, respectively:
C fs
C fs
 2C fs
S fs 2
(C1) fs  f fs Ffs  b K fs
 q fs  fs D fs
 sf C fs  Csl  f fs b
2
t
x
t
x



SHARMA ET AL.







18 of 23

Water Resources Research

10.1029/2020WR029370

CCsl
 2CslS sl Dsl sf Csl  C fs  im  Csl  Cim   fsl  b sl 2(C2)sl  fsl Fsl b Ksl  tsl  qsl
2
x
t
x





Sim2
Cim
(C3)
im  fim Fim b Kim
 fim bim  Csl  Cim 
t
t





The rate-limited sorption kinetics in the fast, slow, and immobile regions is represented as follows,
respectively:S fs2
(C4)k fs  1  Ffs K fsC fs  S fs2 


t





Ssl2
(C5)
ksl 1  Fsl  K slCsl  Ssl2 

t
Sim2
(C6)
kim 1  Fim  KimCim  Sim2 tNow, we convert the DPNE model into the MPNE model. First, the immobile region in the slow domain isturned off by substituting ωim = 0. As a result, Equation C3 is eliminated, and Equation C2 is modified asfollows:
C
Csl
 2CslS sl Dsl sf Csl  C fs  fsl b sl 2
(C7)sl  fsl Fsl b Ksl  tsl  qsl
x
t
x 2





Second, the slow region Equation C7 is converted into an equation describing the immobile region by substituting qsl = 0 (which also leads to Dsl = 0):
C
S
(C8)
sl  fsl Fsl b K sl  tsl sf C fs  Csl  fsl b slt 2

S fs 2Sand sl 2 are calculated from Equatttions C4 and C5, respectively. The conceptual model of the MPNE model is shown in Figure C1.

Equations C1 and C8 represent the MPNE model, in which

Figure C1. Conceptualization of the MPNE model.

SHARMA ET AL.

19 of 23

Water Resources Research

10.1029/2020WR029370

Comparison With the MIM Model Developed by van Genuchten and Wierenga (1976)
The DPNE model can be further simplified into the MIM model. The MIM model developed by van
Genuchten and Wierenga (1976) considers only instantaneous equilibrium sorption, that is, kfs = ksl = 0 in
Equations C4 and C5 and Ffs = Fsl = 1 in Equations C1 and C8, resulting in:
C fs
C fs
 2C fs
(C9) fs  f fs b K fs
 q fs   fs D fs
 sf C fs  Csl
t
x
x 2









C
sf C fs  Csl
(C10)
sl  fsl b K sl  tsl Additionally, ffs=(1-ffs). Equations C9 and C10 represent the MIM model with instantaneous sorption. The
simplified conceptual model is shown in Figure C2.

Figure C2. Conceptualization of the MIM model.

List of Symbols
Subscript “fs”		
Parameters of the fast region
Subscript “sl”		
Parameters of the slow region
Subscript “im”
Parameters of the immobile regionD		
Dispersion coefficient (L2T−1)ωsf		
Mass transfer coefficient between slow and fast regions (T−1)ωim		
Mass transfer coefficient between slow and immobile regions (T−1)θ		PorosityQ		
Darcy’s velocity (LT−1)
V		
Seepage velocity (LT−1)
F		Fraction of sorption sites available
F		
Fraction of sites for which sorption is instantaneousC		
Volume-averaged dissolved concentration in the liquid phase (ML−3)
Concentration in the Laplace domain
C 		
ρb		
Bulk density of the porous medium (ML−3)α		Dispersivity (L)Kd		
Sorption coefficient (L3M−1)
K		
First-order rate-limited sorption coefficient (T−1)
S		
Sorbed phase concentration (ML−3)Pe		
Peclet numberDm		
Damkohler’s numbert0		
Duration of the pulse-type input (T)tmax		
Maximum experimental duration (T)SHARMA ET AL.

20 of 23

Water Resources Research
T		
L		
P		
S 		
C 		
N		
σ2		
		
λ		
 		
 		

10.1029/2020WR029370

Simulation time (T)Length of any location from the inflow (L)Laplace operator
Sorption in the Laplace domain
Concentration in the Laplace domain
Sample size or observations
Residual sum of squares/N
Number of parameters in the model + 1
Eigen values
Eigen vector
Laplace inversion algorithm coefficient

Data Availability Statement
Datasets for this research are included in this paper (and its supplementary information files): Miyamoto
et al. (2003) and referred at s. no. 39.

Acknowledgment
The corresponding author acknowledges the contribution of all co-authors to
help in preparing the manuscript forsuch a prestigious journal.

SHARMA ET AL.

References
Bai, M., Elsworth, D., & Roegiers, J.-C. (1993). Multiporosity/multipermeability approach to the simulation of naturally fractured reservoirs. Water Resources Research, 29(6), 1621–1633. https://doi.org/10.1029/92WR02746
Bai, M., & Roegiers, J.-C. (1994). On the correlation of nonlinear flow and linear transport with application to dual-porosity modeling.
Journal of Petroleum Science and Engineering, 11(2), 63–72. https://doi.org/10.1016/0920-4105(94)90029-9
Bai, M., & Roegiers, J.-C. (1997). Triple-porosity analysis of solute transport. Journal of Contaminant Hydrology, 28(3), 247–266. https://doi.org/10.1016/S0169-7722(96)00086-1
Brusseau, M. L. (1992). Transport of rate-limited sorbing solutes in heterogeneous porous media: Application of a one-dimensional multifactor nonideality model to field data. Water Resources Research, 28(9), 2485–2497. https://doi.org/10.1029/92WR00907
Brusseau, M. L., Jessup, R. E., & Rao, P. S. C. (1989). Modeling the transport of solutes influenced by multiprocess nonequilibrium. Water
Resources Research, 25(9), 1971–1988. https://doi.org/10.1029/WR025i009p01971
Brusseau, M. L., & Rao, P. S. C. (1989). The influence of sorbate-organic matter interactions on sorption nonequilibrium. Chemosphere,18(9–10), 1691–1706. https://doi.org/10.1016/0045-6535(89)90453-0
Brusseau, M. L., & Srivastava, R. (1997). Nonideal transport of reactive solutes in heterogeneous porous media 2. Quantitative analysis of the Borden natural-gradient field experiment. Journal of Contaminant Hydrology, 28(1–2), 115–155. https://doi.org/10.1016/
S0169-7722(97)00036-3
Carnahan, C. L., & Remer, J. S. (1984). Nonequilibrium and equilibrium sorption with a linear sorption isotherm during masstransport through an infinite porous medium: Some analytical solutions. Journal of Hydrology, 73(3–4), 227–258. https://doi.
org/10.1016/0022-1694(84)90002-7
Ceriotti, G. A, Russian, D., & Porta, G. (2019). Russian, a double-continuum transport model for segregated porous media Derivationand sensitivity analysis-driven calibration. Advances in Water Resources, 128, 206–217. https://doi.org/10.1016/j.advwatres.2019.04.003
Ciriello, V., Guadagnini, A., Di Federico, V., Edery, Y., & Berkowitz, B. (2013). Comparative analysis of formulations for conservativetransport in porous media through sensitivity-based parameter calibration. Water Resources Research, 49(9), 5206–5220. https://doi.
org/10.1002/wrcr.20395
Coats, K. H., & Smith, B. D. (1964). Dead-end pore volume and dispersion in porous media. Society of Petroleum Engineers Journal, 4(1),73–84. https://doi.org/10.2118/647-PA
Cukier, R. I., Fortuin, C. M., Shuler, K. E., Petschek, A. G., & Schaibly, J. H. (1973). Study of the sensitivity of coupled reaction systems to
uncertainties in rate coefficients. I Theory. Journal of Chemical Physics, 59(8), 3873–3878. https://doi.org/10.1063/1.1680571
Cukier, R. I., Schaibly, J. H., & Shuler, K. E. (1975). Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients.
III. Analysis of the approximations. Journal of Chemical Physics, 63(3), 1140–1149. https://doi.org/10.1063/1.431440
de Hoog, F. R., Knight, J. H., & Stokes, A. N. (1982). An improved method for numerical inversion of Laplace transforms. SIAM Journal on
Scientific and Statistical Computing, 3(3), 357–366. https://doi.org/10.1137/0903022
Dell’Oca, A., Riva, M., & Guadagnini, A. (2020). Global sensitivity analysis for multiple interpretive models with uncertain parameters.
Water Resources Research, 56(2), 0–3. https://doi.org/10.1029/2019WR025754
Doherty, J., & Christensen, S. (2011). Use of paired simple and complex models to reduce predictive bias and quantify uncertainty. Water
Resources Research, 47(12), 1–21. https://doi.org/10.1029/2011WR010763
Dousset, S., Thevenot, M., Pot, V., Šimunek, J., & Andreux, F. (2007). Evaluating equilibrium and non-equilibrium transport of bromide and isoproturon in disturbed and undisturbed soil columns. Journal of Contaminant Hydrology, 94(3–4), 261–276. https://doi.
org/10.1016/j.jconhyd.2007.07.002
Fajraoui, N., Ramasomanana, F., Younes, A., Mara, T. A., Ackerer, P., & Guadagnini, A. (2011). Use of global sensitivity analysis andpolynomial chaos expansion for interpretation of nonreactive transport experiments in laboratory-scale porous media. Water Resources
Research, 47(2), 1–14. https://doi.org/10.1029/2010WR009639
Gao, G., Feng, S., Zhan, H., Huang, G., & Mao, X. (2009). Evaluation of anomalous solute transport in a large heterogeneous soil column
with mobile-immobile model. Journal of Hydrologic Engineering, 14(9), 966–974. https://doi.org/10.1061/(asce)he.1943-5584.0000071
Gao, G., Zhan, H., Feng, S., Fu, B., Ma, Y., & Huang, G. (2010). A new mobile-immobile model for reactive solute transport with scale-dependent dispersion. Water Resources Research, 46(8), 1–16. https://doi.org/10.1029/2009WR008707

21 of 23

Water Resources Research

10.1029/2020WR029370

Gaudet, J. P., Jégat, H., Vachaud, G., & Wierenga, P. J. (1977). Solute transfer, with exchange between mobile and stagnant water, throughunsaturated sand1. Soil Science Society of America Journal, 41(4), 665. https://doi.org/10.2136/sssaj1977.03615995004100040009x
Gerke, H. H., & van Genuchten, M. T. (1993). A dual-porosity model for simulating the preferential movement of water and solutes instructured porous media. Water Resources Research, 29(2), 305–319. https://doi.org/10.1029/92WR02339
Gerke, H. H., & van Genuchten, M. T. (1996). Macroscopic representation of structural geometry for simulating water and solute movement in dual-porosity media. Advances in Water Resources, 19(6), 343–357. https://doi.org/10.1016/0309-1708(96)00012-7
Goltz, M. N., & Roberts, P. V. (1986). Three-dimensional solutions for solute transport in an infinite medium with mobile and immobile
zones. Water Resources Research, 22(7), 1139–1148. https://doi.org/10.1029/WR022i007p01139
Gwo, J. P., Jardine, P. M., Wilson, G. V., & Yeh, G. T. (1995). A multiple-pore-region concept to modeling mass transfer in subsurface media.
Journal of Hydrology, 164(1–4), 217–237. https://doi.org/10.1016/0022-1694(94)02555-P
Haggerty, R., & Gorelick, S. M. (1995). Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale
heterogeneity. Water Resources Research, 31(10), 2383–2400. https://doi.org/10.1029/95WR10583
Hamby, D. M. (1994). A review of techniques for parameter sensitivity analysis of environmental models. Environmental Monitoring andAssessment, 32(2), 135–154. https://doi.org/10.1007/BF00547132
Helton, J. C. (1993). Uncertainty and sensitivity analysis techniques for use in performance assessment for radioactive waste disposal.
Reliability Engineering & System Safety, 42(2–3), 327–367. https://doi.org/10.1016/0951-8320(93)90097-I
Joshi, N., Ojha, C. S. P., & Sharma, P. K. (2012). A nonequilibrium model for reactive contaminant transport through fractured porous
media: Model development and semianalytical solution. Water Resources Research. 48. 1–19. https://doi.org/10.1029/2011WR011621
Kartha, S. A., & Srivastava, R. (2008). Effect of slow and fast moving liquid zones on solute transport in porous media. Transport in PorousMedia, 75(2), 227–247. https://doi.org/10.1007/s11242-008-9219-7
Kendall, M., & Stuart, A. (1977). The Advanced Theory of Statistics (Vol. 1). Distribution theoryGriffin
Köhne, J. M., Köhne, S., & Šimůnek, J. (2006). Multi-process herbicide transport in structured soil columns: Experiments and model analysis. Journal of Contaminant Hydrology, 85(1–2), 1–32. https://doi.org/10.1016/j.jconhyd.2006.01.001
Leij, F. J., & Bradford, S. A. (2013). Colloid transport in dual-permeability media. Journal of Contaminant Hydrology, 150, 65–76. https://doi.org/10.1016/j.jconhyd.2013.03.010
Leij, F. J., Toride, N., Field, M. S., & Sciortino, A. (2012). Solute transport in dual-permeability porous media. Water Resources Research,48(4), 1–13. https://doi.org/10.1029/2011WR011502
Leij, F. J., Toride, N., & van Genuchten, M. T. (1993). Analytical solutions for non-equilibrium solute transport in three-dimensional porous
media. Journal of Hydrology, 151(2–4), 193–228. https://doi.org/10.1016/0022-1694(93)90236-3
Liang, X., Zhang, Y. K., Liu, J., Ma, E., & Zheng, C. (2019). Solute transport with linear reactions in porous media with layered structure: A
semi-analytical model. Water Resources Research, 55(6), 5102–5118. https://doi.org/10.1029/2019WR024778
McRae, G. J., Tilden, J. W., & Seinfeld, J. H. (1982). Global sensitivity analysis-a computational implementation of the Fourier Amplitude
Sensitivity Test (FAST). Computers & Chemical Engineering, 6(1), 15–25. https://doi.org/10.1016/0098-1354(82)80003-3
Miyamoto, T., Annaka, T., & Chikushi, J. (2003). Soil aggregate structure effects on dielectric permittivity of an andisol measured by time
domain reflectometry. Vadose Zone Journal, 2(1), 90–97. https://doi.org/10.2136/vzj2003.900010.2113/2.1.90
Neville, C. J., Ibaraki, M., & Sudicky, E. A. (2000). Solute transport with multiprocess nonequilibrium: A semi-analytical solution approach.
Journal of Contaminant Hydrology, 44(2), 141–159. https://doi.org/10.1016/S0169-7722(00)00094-2
Nguyen, A.-T., & Reiter, S. (2015). A performance comparison of sensitivity analysis methods for building energy models. Building Simulation, 8(6), 651–664. https://doi.org/10.1007/s12273-015-0245-4
Norton, J. P. (2008). Algebraic sensitivity analysis of environmental models. Environmental Modelling & Software, 23(8), 963–972. https://doi.org/10.1016/j.envsoft.2007.11.007
Pianosi, F., Beven, K., Freer, J., Hall, J. W., Rougier, J., Stephenson, D. B., & Wagener, T. (2016). Sensitivity analysis of environmental
models: A systematic review with practical workflow. Environmental Modelling & Software, 79, 214–232. https://doi.org/10.1016/j.
envsoft.2016.02.008
Pot, V., Šimůnek, J., Benoit, P., Coquet, Y., Yra, A., & Martínez-Cordón, M.-J. (2005). Impact of rainfall intensity on the transport of two
herbicides in undisturbed grassed filter strip soil cores. Journal of Contaminant Hydrology, 81(1–4), 63–88. https://doi.org/10.1016/j.
jconhyd.2005.06.013
Saltelli, A., Tarantola, S., & Chan, K. P.-S. (1999). A quantitative model-independent method for global sensitivity analysis of model output.
Technometrics, 41(1), 39–56. https://doi.org/10.1080/00401706.1999.10485594
Selim, H. M., & Ma, L. (1998). Physical nonequilibrium in soils-Modeling and application. Ann Arbor Press.
Selim, H. M., Ma, L., & Zhu, H. (1999). Predicting solute transport in soils second-order two-site models. Soil Science Society of America
Journal, 63(4), 768–777. https://doi.org/10.2136/sssaj1999.634768x
Sharma, P. K., Ojha, C. S. P., Swami, D., Joshi, N., & Shukla, S. K. (2015). Semi-analytical solutions of multiprocessing non-equilibriumtransport equations with linear and exponential distance-dependent dispersivity. Water Resources Management, 29(14), 5255–5273.
https://doi.org/10.1007/s11269-015-1116-6
Sharma, P. K., Shukla, S. K., Choudhary, R., & Swami, D. (2016). Modeling for solute transport in mobile-immobile soil column experiment. ISH Journal of Hydraulic Engineering, 22(2), 204–211. https://doi.org/10.1080/09715010.2016.1155181
Šimůnek, J., Jarvis, N. J., Van Genuchten, M. T., & Gärdenäs, A. (2003). Review and comparison of models for describing non-equilibrium
and preferential flow and transport in the vadose zone. Journal of Hydrology, 272, 14–35. https://doi.org/10.1016/S0022-1694(02)00252-4
Šimůnek, J., & van Genuchten, M. T. (2008). Modeling nonequilibrium flow and transport processes using HYDRUS. Vadose Zone Journal,
7(2), 782–797. https://doi.org/10.2136/vzj2007.0074
Šimůnek, J., van Genuchten, M. T., & Šejna, M. (2008). Development and applications of the HYDRUS and STANMOD software packages
and related codes. Vadose Zone Journal, 7(2), 587–600. https://doi.org/10.2136/vzj2007.0077
Skopp, J., Gardner, W. R., & Tyler, E. J. (1981). Solute movement in structured soils: Two-region model with small interaction. Soil Science
Society of America Journal, 45(5), 837–842. https://doi.org/10.2136/sssaj1981.03615995004500050002x
Soboĺ, I. M. (1990). Quasi-Monte Carlo methods. Progress in Nuclear Energy, 24(1–3), 55–61. https://doi.org/10.1016/0149-1970(90)90022-W
Sobol, I. M. (1993). Sensitivity analysis for non-linear mathematical models. Mathematical Modelling in Civil Engineering, 1, 407–414.
https://ci.nii.ac.jp/naid/10027137978/en/
Srivastava, R., & Brusseau, M. L. (1996). Nonideal transport of reactive solutes in heterogeneous porous media: 1. Numerical model development and moments analysis. Journal of Contaminant Hydrology, 24(2), 117–143. https://doi.org/10.1016/S0169-7722(96)00039-3
Sudicky, E. A., & Frind, E. O. (1982). Contaminant transport in fractured porous media: Analytical solutions for a system of parallel fractures. Water Resources Research, 18(6), 1634–1642. https://doi.org/10.1029/WR018i006p01634

SHARMA ET AL.

22 of 23

Water Resources Research

10.1029/2020WR029370

Swami, D., Sharma, A., Sharma, P. K., & Shukla, D. P. (2016). Predicting suitability of different scale-dependent dispersivities for reactive
solute transport through stratified porous media. Journal of Rock Mechanics and Geotechnical Engineering, 8(6), 921–927. https://doi.
org/10.1016/j.jrmge.2016.07.005
Swami, D., Sharma, P. K., & Ojha, C. S. P. (2014). Simulation of experimental breakthrough curves using multiprocess non-equilibrium
model for reactive solute transport in stratified porous media. Sadhana, 39(6), 1425–1446. https://doi.org/10.1007/s12046-014-0287-9
Swami, D., Sharma, P. K., Ojha, C. S. P., Guleria, A., & Sharma, A. (2018). Asymptotic behavior of mass transfer for solute transport throughstratified porous medium. Transport in Porous Media, 124(3), 699–721. https://doi.org/10.1007/s11242-018-1090-6
Tang, D. H., Frind, E. O., & Sudicky, E. A. (1981). Contaminant transport in fractured porous media: Analytical solution for a single fracture. Water Resources Research, 17(3), 555–564. https://doi.org/10.1029/WR017i003p00555
Tian, W. (2013). A review of sensitivity analysis methods in building energy analysis. Renewable and Sustainable Energy Reviews, 20,
411–419. https://doi.org/10.1016/j.rser.2012.12.014
Toride, N., Inoue, M., & Leij, F. J. (2003). Hydrodynamic dispersion in an unsaturated dune sand. Soil Science Society of America Journal,
67(3), 703. https://doi.org/10.2136/sssaj2003.0703
Toride, N., Leij, F. J., & van Genuchten, M. T. (1993). A comprehensive set of analytical solutions for nonequilibrium solute transport with
first-order decay and zero-order production. Water Resources Research, 29(7), 2167–2182. https://doi.org/10.1029/93WR00496
Turanyi, T. (1990). Sensitivity analysis of complex kinetic systems. Tools and applications. Journal of Mathematical Chemistry, 5(3), 203–
248. https://doi.org/10.1007/BF01166355
van Genuchten, M. T., & Wagenet, R. J. (1989). Two-site/two-region models for pesticide transport and degradation: Theoretical development and analytical solutions. Soil Science Society of America Journal, 53(5), 1303–1310. https://doi.org/10.2136/
sssaj1989.03615995005300050001x
van Genuchten, M. T., & Wierenga, P. J. (1976). Mass transfer studies in sorbing porous Media I. Analytical solutions. Soil Science Societyof America Journal, 40(4), 473–480. https://doi.org/10.2136/sssaj1976.03615995004000040011x
van Genuchten, M. T., Wierenga, P. J., & O'Connor, G. A. (1977). Mass transfer studies in sorbing porous media: III. Experimental evaluation with 2,4,5-T. Soil Science Society of America Journal, 41(2), 278–285. https://doi.org/10.2136/sssaj1977.03615995004100020023x
Wainwright, H. M., Finsterle, S., Jung, Y., Zhou, Q., & Birkholzer, J. T. (2014). Making sense of global sensitivity analyses. Computers &Geosciences, 65, 84–94. https://doi.org/10.1016/j.cageo.2013.06.006
Warren, J. E., & Root, P. J. (1963). The behavior of naturally fractured reservoirs. Society of Petroleum Engineers Journal, 3(3), 245–255.
https://doi.org/10.2118/426-pa
Zimmerman, R. W., Chen, G., Hadgu, T., & Bodvarsson, G. S. (1993). A numerical dual-porosity model with semianalytical treatment of
fracture/matrix flow. Water Resources Research, 29(7), 2127–2137. https://doi.org/10.1029/93WR00749

SHARMA ET AL.

23 of 23