The Linear Boltzmann Equation in
Column Experiments of Porous Media

Kenji Amagai, Motoko Yamakawa,
Manabu Machida & Yuko Hatano

Transport in Porous MediaISSN 0169-3913
Transp Porous Med
DOI 10.1007/s11242-020-01393-1

1 23

Your article is protected by copyright andall rights are held exclusively by Springer
Nature B.V.. This e-offprint is for personaluse only and shall not be self-archived
in electronic repositories. If you wish toself-archive your article, please use theaccepted manuscript version for posting on
your own website. You may further deposit
the accepted manuscript version in any
repository, provided it is only made publicly
available 12 months after official publicationor later and provided acknowledgement isgiven to the original source of publication
and a link is inserted to the published article
on Springer's website. The link must be
accompanied by the following text: "The final
publication is available at link.springer.com”.

1 23

Author's personal copy
Transport in Porous Media
https://doi.org/10.1007/s11242-020-01393-1

The Linear Boltzmann Equation in Column Experimentsof Porous Media
Kenji Amagai1 · Motoko Yamakawa1 · Manabu Machida2 · Yuko Hatano1
Received: 5 March 2019 / Accepted: 19 February 2020
© Springer Nature B.V. 2020

Abstract
The use of the linear Boltzmann equation is proposed for transport in porous media in acolumn. By column experiments, we show that the breakthrough curve is reproduced by
the linear Boltzmann equation. The advection–diffusion equation is derived from the linear
Boltzmann equation in the asymptotic limit of large propagation distance and long time.
Keywords Anomalous diffusion · Linear Boltzmann transport · Column experiment ·Porous media

1 Introduction
Mass transport in porous media is one of the central topics in hydrology. For example, masstransport has been intensively studied in the contexts of carbon dioxide capture and storage(Benson and Orr 2008) and enhanced oil recovery (Thomas 2008). Moreover, migration of
the soil pollution of radioactive cesium in Fukushima, Japan, which is considered as masstransport in porous media, is an urgent issue (Matsuda et al. 2015). The phenomenon of masstransport in porous media has been analyzed by the advection–diffusion equation (ADE) since
1950s. However, different examples in which ADE is not applicable have been found. Adams
and Gelhar showed that the mass transport predicted by ADE was significantly different frommeasured values in a field scale experiment (Adams and Gelhar 1992). In the laboratory scale,
Berkowitz, Scher, and Silliman found that the transport in a laboratory flow cell does not obey
ADE, but it shows anomalous diffusion (or non-Fickian transport) and deviates from Fick’slaw (Berkowitz et al. 2000). In such flow of anomalous diffusion, unlike the prediction by
ADE, the peak of concentration barely moves and its distribution has a long tail. The deviationfrom ADE was investigated by a random particle simulation (Kennedy and Lennox 2001).
* Yuko Hatano
hatano@risk.tsukuba.ac.jp
Manabu Machidamachida@hama‑med.ac.jp
1

Graduate School of Systems and Information Engineering, University of Tsukuba,
Tsukuba 305‑7361, Japan

2

Institute for Medical Photonics Research, Hamamatsu University School of Medicine,
Hamamatsu 431‑3192, Japan

13

Vol.:(0123456789)

Author's personal copy
K. Amagai et al.

To reproduce anomalous diffusion, continuous time random walk (CTRW) has been developed since the first proposal by Montroll and Weiss (Montroll and Weiss 1965). Anomalous
diffusion in porous media in different situations can be reproduced by CTRW when waitingtime and jump distance are suitably adjusted. Although CTRW is widely used in the studyof porous media (Berkowitz and Scher 1998; Hatano and Hatano 1998; Levy and Berkowitz
2003; Nissan et al. 2017), the relation between the movement of tracer particles and the waiting time and jump distance of CTRW is not apparent. By taking the correlation in waiting time
into account, correlated continuous time random walk was proposed (Montero and Masoliver2007). In addition to CTRW, the following models were proposed. In tempered anomalous
diffusion, strong anomalous diffusivity is suppressed by the introduction of an exponential
cutoff for waiting times (Meerschaert et al. 2008). The stochastic hydrology was developed(Gelhar 1986; Rubin 2003). The mobile/immobile (MIM) model and fractional MIM model
divide mass transport into the moving part and immovable part (Van Genuchten and Wierenga
1976; Schumer et al. 2003). It is known that diffusion equations with fractional derivatives are
obtained in the asymptotic limit of CTRW (Metzler and Klafter 2000). In fractional advection–dispersion equation, time and spatial derivatives are replaced by fractional derivatives(Benson et al. 2000). The use of fractional derivatives has attracted attention as a theoretical
tool for anomalous diffusion. Fractional derivative models are compared (Wei et al. 2016; Sunet al. 2017), parameters in fractional equations were estimated (Chakraborty et al. 2009; Kelly
et al. 2017), the physical meaning of fractional derivatives was discussed (Liang et al. 2019),
and orders of fractional derivatives were taken to be variables (Sun et al. 2009).
To study transport in porous media, in addition to the concentration in mass transport,
actual flow was visualized (Moroni et al. 2009; Sen et al. 2012) and the flow was reproducedby large-scale simulation (De Anna et al. 2013; Liu et al. 2016). These early studies show that
tracer particles and fluid particles change velocity as a function of time when along the flow
path they make a detour around obstacles in the porous medium.
Although CTRW and fractional equations have achieved certain success, the physicalmeaning of the random walk process and physical background of fractional-order derivatives
are not easy to understand. In this paper, we propose to use the linear Boltzmann equation(LBE), in which the physical scattering process is more apparent. Figure 1 illustrates the transport in porous media. Tracer particles in a porous medium flow along fissures and fractures.

Fig. 1  The transport of tracer
particles in a porous medium.
This transport is regarded asnoninteracting particles which
undergo scatterings by scatterers

13Author's personal copy
The Linear Boltzmann Equation in Column Experiments of Porous…

We regard the process of a tracer particle changing its direction as a scattering. Thus for theensemble of tracer particles, the transport in the porous medium can be viewed as the transport of noninteracting particles which are scattered by scatterers in different directions with acertain probability distribution.
The rest of the paper is organized as follows. In Sect. 2, the linear Boltzmann equation is
introduced. We explore the advection–diffusion approximation in Sect. 3. In Sect. 4, we solve
the linear Boltzmann equation. Section 5 is devoted to column experiments. Conclusions are
given in Sect. 6.

2 Linear Boltzmann Transport
The linear Boltzmann equation, radiative transport equation, or transport equation governs thetransport of noninteracting particles with scattering and absorption. The equation describes
neutron transport in a reactor (Case and Zweifel 1967) and light transport in random media
such as clouds and the interstellar medium (Apresyan and Kravtsov 1996; Chandrasekhar
1960; Ishimaru 1978; Sobolev 1976; Thomas and Stamnes 1999). Let us regard the transportof tracer particles along flow paths as the transport of particles in a homogeneous random
medium. See the schematic figure in Fig. 1. There is a random arrangement of flow paths inbeads and sand in column experiments. At a branching point where several paths split into
different directions, tracer particles are scattered to these directions with certain probabilities.
The average distance that a tracer particle travels along a flow path between branching points
is an effective mean free path v0 ∕𝜎s [see (1) for v0 , 𝜎s]. Indeed, the use of the linear Boltzmann
equation for the flow of tracer particles in porous media was proposed by Williams (1992a, b,
1993a, b).
A remark is necessary for our Boltzmann transport. In recent years, the lattice Boltzmann
method (Succi 2001) has become one of commonly used numerical methods to investigate theflow in porous media. Instead of fluid in porous media, in this work we focus on the Boltzmann transport of tracer particles in the fluid. Hence, the proposed equation (1) is a linear
equation.
Having in mind the transport in column experiments which is described in Sect. 5, we consider the one-dimensional linear Boltzmann equation. Let 𝜇 ∈ [−1, 1] be the cosine of thepolar angle and v0 > 0 be the inherent particle speed. Although the coefficient of the spatial
derivative is v0 𝜇 in Williams’ study (Williams 1992a, b) , in this paper, we take into accountadvection, which will be denoted by u ≥ 0. Hence, the velocity v is written as v = u + v0 𝜇. Let
𝜓(x, 𝜇, t) be the angularly dependent number density of particles at position x ≥ 0 in direction
𝜇 at time t ≥ 0. We can write the linear Boltzmann equation as follows:1𝜎𝜕𝜕𝜓(x, 𝜇, t) d𝜇.
𝜓(x, 𝜇, t) + (u + v0 𝜇) 𝜓 + (𝜎a + 𝜎s )𝜓 = s
𝜕t
𝜕x
2 ∫−1

(1)

where isotropic scattering is assumed in the integral on the right-hand side. Here, 𝜎a > 0
and 𝜎s > 0 are the absorption and scattering coefficients, respectively. Let n0 be the initial
particle number density. Tracer particles are injected to the column at x = 0 in the normal
direction. The initial and boundary conditions are given by𝜓(x, 𝜇, 0) = 0,

(2)

13

Author's personal copy
K. Amagai et al.

𝜓(0, 𝜇, t) = n0 𝛿(𝜇 − 1),

𝜇 > −𝜂,(3)

where 𝛿(𝜇 − 1) is the Dirac delta function and

u
.
v0

𝜂=

(4)

We note that unlike the usual half-range boundary condition of 𝜇 > 0 (Case and Zweifel1967), the boundary value in (3) is specified for 𝜇 ∈ (−𝜂, 1] due to the presence of theadvection. We regard the column as the half space (0 < x < ∞) and have 𝜓 → 0 as x → ∞.
The particle number density n(x, t) is calculated as1n(x, t) =∫−1

𝜓(x, 𝜇, t) d𝜇.

(5)Then, the concentration is given by

C(x, t) = 𝛼n(x, t),(6)where 𝛼 is the mass per particle.

3 Advection–Diffusion Approximation
Let us derive the advection–diffusion equation from the linear Boltzmann equation. According
to the usual procedure of the diffusion approximation (Duderstadt and Martin 1979; Ishimaru
1978), we assume that the angular dependence of 𝜓(x, 𝜇, t) is weak and written as𝜓(x, 𝜇, t) =3
1
C(x, t) +J(x, t)𝜇,
2𝛼
2𝛼v0

(7)

where
1

C(x, t) = 𝛼∫−1𝜓(x, 𝜇, t) d𝜇,(8)𝜇𝜓(x, 𝜇, t) d𝜇.

(9)1J(x, t) = 𝛼v0∫−1

In (7), we implicitly assumed a large propagation distance and long time. If absorption
is strong, 𝜓 becomes zero before the form (7) is achieved. Hence, we assume 𝜎a is small.
Indeed in Sect. 5, we will see that 𝜎a is negligibly small. The approximation in (7) is called
the P1 approximation (Case and Zweifel 1967). Below, we will proceed further and obtain
the diffusion equation. We substitute the form (7) for 𝜓 in (1). By integrating the resulting
equation over 𝜇, we obtain

𝜕C
𝜕C 𝜕J
+u
+
+ 𝜎a C = 0.
𝜕t𝜕x𝜕x
By integrating the equation over 𝜇 after multiplying 𝜇, we obtain

13

(10)

Author's personal copy
The Linear Boltzmann Equation in Column Experiments of Porous…

1 𝜕J
u 𝜕J v0 𝜕C 𝜎a + 𝜎s+++
J = 0.
v0 𝜕t v0 𝜕x
3 𝜕x
v0

(11)

We assume that the first term and second term of the above equation are small and can be
dropped. Then, we have

𝜕C
,
𝜕x

(12)

v0 𝓁 ∗
,
3

(13)

v
v0
≃ 0.
𝜎a + 𝜎s𝜎s(14)

J = −D
where

D=

𝓁∗ =

Hence, we obtain the following advection–diffusion equation by substituting (12) into (10):

𝜕2 C
𝜕C
𝜕C
−D 2 +u+ 𝜎a C = 0,
𝜕t
𝜕x
𝜕x

(15)

For the photon diffusion equation, 𝓁 ∗ is known to be independent of the absorption coefficient (Furutsu and Yamada 1994). The diffusion approximation holds when x is sufficiently
larger than the transport mean free path 𝓁 ∗.
The diffusion equation (15) reduces to the usual advection–diffusion equation when 𝜎a = 0.
Since we found 𝜎a ≈ 0 in our column experiments, the solution by Ogata and Banks (Ogata
and Banks 1961) is used in this paper when the concentration calculated from the linear Boltzmann equation is compared to the advection–diffusion equation. In Appendix 2, the solution
to (15) is presented.
Finally, we consider the boundary condition. Let us substitute (7) for 𝜓(x, 𝜇, t) in theboundary condition (3). By integrating both sides of the resulting equation over 𝜇 from −𝜂 to1, we obtain

C(0, t) −2𝛼n03D𝜂(1 − 𝜂) D 𝜕C.
(0, t) =
2
u 𝜕x
1+𝜂

(16)

Thus, we obtain the Robin boundary condition. In this paper, we neglect the flux term andwrite

C(0, t) = C0 ,(17)

where C0 > 0 denotes the boundary source. Since the diffusion approximation does nothold on the boundary, the actual proportionality constant between 𝛼n0 and C0 is a fitting
parameter. We introduce

𝛽=

𝛼n0
.
C0

(18)

13

Author's personal copy
K. Amagai et al.

4 Concentration from the Linear Boltzmann Equation
Let us consider the Laplace transform:

∫0

𝜓(x,
̂ 𝜇, p) =∞e−pt 𝜓(x, 𝜇, t) dt.

(19)

We introduce

𝜇t =𝜎a + 𝜎s + p
,
v0

𝜇s =

𝜎s
.
v0

(20)

Then, our transport equation is written as(𝜂 + 𝜇)1𝜇𝜕
𝜓(x,̂ 𝜇, p) d𝜇,𝜓(x,̂ 𝜇, p) + 𝜇t 𝜓(x,
̂ 𝜇, p) = s
𝜕x
2 ∫−1

(21)

with the boundary condition:𝜓(0,̂ 𝜇, p) =n0
𝛿(𝜇 − 1),p

𝜇 > −𝜂.

(22)

We will compute 𝜓̂ with the analytical discrete ordinates (ADO) method (Siewert andWright 1999; Barichello et al. 2000; Barichello and Siewert 2001). Unlike the standard
ADO, the coefficient 𝜇t is a complex number since the variable p is complex, and moreover, the boundary condition is given for 𝜇 > −𝜂 instead of 𝜇 > 0. We begin by decomposing 𝜓̂ into the ballistic term 𝜓̂ b and scattering term 𝜓̂ s as𝜓(x,̂ 𝜇, p) = 𝜓̂b (x, 𝜇, p) + 𝜓̂s (x, 𝜇, p).
Here, 𝜓̂ b obeys(𝜂 + 𝜇)𝜕𝜓̂ + 𝜇t 𝜓̂b = 0,
𝜕x bx > 0, −1 ≤ 𝜇 ≤ 1,

(23)

(24)

with the boundary condition𝜓̂b (0, 𝜇, p) =n0
𝛿(𝜇 − 1),p

𝜇 > −𝜂,(25)

and 𝜓̂ s obeys(𝜂 + 𝜇)1𝜇𝜕𝜓̂s + 𝜇t 𝜓̂s = s𝜓̂ d𝜇 + q,𝜕x2 �−1 sx > 0, −1 ≤ 𝜇 ≤ 1,

(26)

with the boundary condition𝜓̂s (0, 𝜇, p) = 0,

𝜇 > −𝜂,

(27)

where

q(x, 𝜇, p) =We obtain

13

n 𝜇
𝜇s 1
𝜓̂ (x, 𝜇, p) d𝜇 = 0 s e−x𝜇t ∕(𝜂+𝜇) .
2 ∫−1 b2p

(28)

Author's personal copy
The Linear Boltzmann Equation in Column Experiments of Porous…

𝜓̂b (x, 𝜇, p) =n0 −x𝜇 ∕(𝜂+𝜇)e t𝛿(𝜇 − 1).
p

(29)

Let us drop p and write 𝜓̂ s (x, 𝜇) = 𝜓̂ s (x, 𝜇, p) and q(x, 𝜇) = q(x, 𝜇, p). For the computationof 𝜓̂ s, we discretize the integral by the Gauss–Legendre quadrature and obtain𝜕𝜓̂ (x, 𝜇i ) + 𝜇t 𝜓̂s (x, 𝜇i )
𝜕x s
N
]
𝜇 ∑ [
wj 𝜓̂s (x, 𝜇j ) + 𝜓̂s (x, −𝜇j ) + q(x, 𝜇i ),
= s
2 j=1

(𝜂 + 𝜇i )(30)

where 𝜇i , wi (i = 1, 2, … , 2N ) are abscissas and weights, respectively. These 𝜇i , wiwith 0 < 𝜇1 < ⋯ < 𝜇N < 1 and 𝜇N+i = −𝜇i (i = 1, … , N ) can be calculated by the
Golub–Welsch algorithm (Golub and Welsch 1969). In the numerical calculation in Sect. 5,
we set(31)

N = 30.

Furthermore, we introduce N𝜂 as the largest integer such that −𝜂 < 𝜇N𝜂 . The scattering part𝜓̂ s is obtained as𝜓̂s (x, 𝜇i ) =

2N
∑

∫
j=1 0

∞

(32)

G(x, 𝜇i ;x� , 𝜇j )q(x� , 𝜇j ) dx� ,where the Green’s function defined for each p satisfies(𝜂 + 𝜇i )2N
𝜇 ∑𝜕G(x, 𝜇i ;x0 , 𝜇i0 ) + 𝜇t G(x, 𝜇i ;x0 , 𝜇i0 ) = sw G(x, 𝜇j ;x0 , 𝜇i0 )𝜕x2 j=1 j(33)

+ 𝛿(x − x0 )𝛿ii0 ,where 𝛿ii0 is the Kronecker delta. Here, the boundary condition is given by

G(0, 𝜇i ;x0 , 𝜇i0 ) = 0,𝜇i > 𝜇N𝜂 .

(34)

Let us consider the following homogeneous equation to calculate the Green’s function withADO (Siewert and Wright 1999; Barichello et al. 2000; Barichello and Siewert 2001).
Note that the equation will be solved for each p:((𝜂 + 𝜇i )

2N
)
𝜇 ∑
𝜕
+ 𝜇t 𝜓(x,̂ 𝜇i ) = sw 𝜓(x,
̂ 𝜇j ).
𝜕x
2 j=1 j(35)

We note that 𝜓̂ depends on p through 𝜇t . With separation of variables, we can write 𝜓̂ as𝜓(x,̂ 𝜇i ) = 𝜙(𝜈, 𝜇i )e−x∕𝜈 ,(36)

where 𝜈 is the separation constant. The function 𝜙(𝜈, 𝜇i ) satisfies the normalization
condition:
2N
∑
i=1

wi 𝜙(𝜈, 𝜇i ) =

N
∑

(
)
wi 𝜙(𝜈, 𝜇i ) + 𝜙(𝜈, −𝜇i ) = 1.

(37)

i=1

13

Author's personal copy
K. Amagai et al.

For the integral equation in which the sum on the right-hand side of (35) is replaced by the1𝜇integral 2s ∫−1 𝜓̂ d𝜇, this 𝜙 is called the singular eigenfunction (Case 1960). However, by
the discretization, we obtain

𝜇s 𝜈
1
,
2 𝜇t 𝜈 − 𝜇i − 𝜂𝜙(𝜈, 𝜇i ) =(38)

assuming 𝜈 ≠ (𝜇i + 𝜂)∕𝜇t . It is known that 𝜇 ≠ 𝜇i ∕𝜇t if 𝜂 = 0 and p is real (Siewert andWright 1999). We can show the following orthogonality relation:2N∑

wi (𝜇i + 𝜂)𝜙(𝜈, 𝜇i )𝜙(𝜈 � , 𝜇i ) = N(𝜈)𝛿𝜈𝜈 � ,

(39)

i=1

where

N(𝜈) =2N∑

wi (𝜇i + 𝜂)𝜙(𝜈, 𝜇i )2 .

(40)

i=1

We can find 2N eigenvalues 𝜈 = 𝜈n (n = 1, 2, … , 2N ). Moreover, there are N𝜂 eigenvalues
with positive real parts. See Appendix 1 for the eigenvalues.
Taking the fact that G vanishes as x → ∞, we can write

G(x, 𝜇i ;x0 , 𝜇i0 ) = G0 (x, 𝜇i ;x0 , 𝜇i0 )∑+B(𝜈;x0 , 𝜇i0 )𝜙(𝜈, 𝜇i )e−x∕𝜈 ,(41)ℜ𝜈n >0where coefficients B(𝜈;x0 , 𝜇i0 ) will be determined later. The free-space Green’s function G0
is calculated as

G0 (x, 𝜇i ; x0 , 𝜇i0 ) = ±∑wi0

N(𝜈n )±ℜ𝜈 >0𝜙(𝜈n , 𝜇i0 )𝜙(𝜈n , 𝜇i )e−(x−x0 )∕𝜈n ,

(42)

n

where upper signs are chosen for x > x0 and lower signs are used for x < x0. Since G vanishes at x = 0 for 𝜇i > 𝜇N𝜂 , we have from (41):∑B(𝜈; x0 , 𝜇i0 )𝜙(𝜈n , 𝜇i ) =ℜ𝜈n >0∑wi0

N(𝜈n )ℜ𝜈 <0

𝜙(𝜈n , 𝜇i0 )𝜙(𝜈n , 𝜇i )ex0 ∕𝜈n .

(43)

n

Let us multiply (43) by exp(−x0 𝜇t ∕(𝜂 + 1)), then integrate both sides with x0 and take thesum with respect to i. We obtain∑E(𝜈n )𝜙(𝜈n , 𝜇i ) =ℜ𝜈n >0∑−𝜈n (𝜂 + 1)
𝜙(𝜈n , 𝜇i ),N(𝜈n )(𝜂 + 1 − 𝜈n 𝜇t )
ℜ𝜈 <0

(44)

n

where

E(𝜈n ) =2N∑∫
i =1 0∞B(𝜈n ;x0 , 𝜇i0 )e−𝜇t x0 ∕(𝜂+1) dx0 .

0We can numerically obtain E(𝜈n ) from the linear system (44). Finally, we obtain

13

(45)

Author's personal copy
The Linear Boltzmann Equation in Column Experiments of Porous…

][)( −𝜇 x∕(𝜂+1)
n0 𝜇s ∑
𝜈n
𝜂+1
−x∕𝜈n
t
𝜙(𝜈n , 𝜇i )𝜓̂s (x, 𝜇i ) =e−e
2p ℜ𝜈 >0 N(𝜈n )𝜂 + 1 − 𝜈n 𝜇t
n[]
n0 𝜇s ∑
𝜈n
𝜂+1
−𝜇t x∕(𝜂+1)+𝜙(𝜈n , 𝜇i )e
2p ℜ𝜈 <0 N(𝜈n )𝜂 + 1 − 𝜈n 𝜇t
n
n 𝜇 ∑E(𝜈n )𝜙(𝜈n , 𝜇i )e−x∕𝜈n .
+ 0 s
2p ℜ𝜈 >0(46)n

In this way, 𝜓̂ s is computed using ADO.
The Laplace transform of the particle number density n̂ (x, p) is obtained as
1n̂ (x, p) =

∫−1

𝜓(x,
̂ 𝜇, p) d𝜇11𝜓̂ (x, 𝜇, p) d𝜇
∫−1 s
∫−1
n
= 0 e−x𝜇t ∕(𝜂+1)p[]( −x𝜇 ∕(𝜂+1)
)
n 𝜇 ∑
𝜈n
𝜂+1
e t
+ 0 s
− e−x∕𝜈n
2p ℜ𝜈 >0 N(𝜈n ) 𝜂 + 1 − 𝜈n 𝜇t
n[]
n0 𝜇s ∑
𝜈n
𝜂+1
e−x𝜇t ∕(𝜂+1)+
2p ℜ𝜈 <0 N(𝜈n ) 𝜂 + 1 − 𝜈n 𝜇t
n
n0 𝜇s ∑+E(𝜈n )e−x∕𝜈n .
2p ℜ𝜈 >0=𝜓̂b (x, 𝜇, p) d𝜇 +

(47)

n

We note that the particle number density n(x, t) is given by𝛾+i∞n(x, t) =1
ept n̂ (x, p) dp.
2𝜋i ∫𝛾−i∞(48)

The Bromwich integral in the inverse Laplace transform is numerically evaluated by thetrapezoidal rule. In general, deformation of the contour can be considered for the inversion(Weideman and Trefethen 2007). We found 𝛾 = 0.04 is suitable. Hence, the concentration
is obtained as{
𝛼n0 𝜇s 𝛾+i∞ ept 2 −x𝜇 ∕(𝜂+1)e t
C(x, t) =4𝜋i ∫𝛾−i∞ p 𝜇s[]∑
( −x𝜇 ∕(𝜂+1)
)
𝜈n
𝜂+1
−x∕𝜈n
t
+
e
−e
N(𝜈n ) 𝜂 + 1 − 𝜈n 𝜇tℜ𝜈n >0
]
[
(49)
∑
𝜈n
𝜂+1
e−x𝜇t ∕(𝜂+1)+N(𝜈n ) 𝜂 + 1 − 𝜈n 𝜇tℜ𝜈n <0}∑+E(𝜈n )e−x∕𝜈n dp.
ℜ𝜈n >013

Author's personal copy
K. Amagai et al.

5 Column Experiments
Column experiments are often used for the study of transport in random media (Cortis
et al. 2004). We here present results of a series of tracer breakthrough experiments conducted in a one-dimensional flow field. Fluid with tracer particles is injected by the peristaltic pump (MP-1000, Eyela). The flow rate of injected tracer solution is controlled
by the peristaltic pump. Ultra-pure water is used. First, water is injected for 24 hours ormore until steady flow field is achieved. Next, tracer solution is injected to displace thefresh water. The discharged solution is collected by the fraction collector (CHF161RA,
Advantec) at regular intervals. Column experiments are performed at room temperatureof 25–26 ◦ C.
Tracer experiments are conducted until the discharge concentration becomes equal tothe influent concentration. For the sake of the accuracy of measurements, we correct themeasurement time error from the tube volume and inlet volume of the column by subtracting the time lag due to the switchover from the breakthrough time.

5.1 Column Experiments with Non‑adsorbed Solute
First, we use non-adsorbed solute. We use two different sizes of glass beads made withsoda glass: (Run 1) smaller beads with diameter ranging from 0.177 to 0.250 mm and(Run 2) larger beads with diameter from 0.500 to 0.710 mm . The column length is
20.0 cm with the internal diameter of the column 2.9 cm . In order to prevent entrainmentof air, we carry out beads packing under saturated conditions, in which the water surface is kept at a constant height from the top of the beads layer. We stirred beads whilefilling in order to remove small bubbles attached to the beads.
The inlet and the outlet ends of the column are separated from the porous medium
of a glass filter with the pore size 10 μm, originally for the purpose of the liquid chromatography. The glass filter serves to enhance tracer mixing and to dampen
small flow pulses on injection. The tracer concentration is measured by the EC meter(LAQUAtwin-EC-33B, Horiba). The precision of the meter is ±2%. The temperaturecompensation circuit is installed.
We prepare the NaCl solution of 0.1% with the salt of purity 99.5% (Wako). Injectionof the solution is done by the peristaltic pump. The measured EC values are given by
the unit of Sv/cm; we made the calibration curve and obtain the concentration.

5.1.1 Run 1: Small Beads
The column is filled with small beads. The measured porosity is 𝜂 = 0.386 . The bedheight is 18.0 cm with section area 6.61 cm2 . The flow rate is 0.958 cm3 ∕min.
Figure 2 shows the chloride breakthrough curve from Run 1. The vertical axis shows
the relative concentration C∕C0 , and the horizontal axis is time t. The measured values from the column experiment (black dots) are compared to C∕C0 calculated by
the linear Boltzmann equation (LBE) (red solid line) and by the ADE model (Ogata
and Banks 1961) (blue dashed line). The parameters in both equations are determined by nonlinear least-squares fitting (i.e., the minimization of the sum of squaresof the residuals between the experimental and theoretical curves); they are obtained
as 𝜎a = 1.0 × 10−8 min−1, 𝜎s = 5.1645 min−1, v0 = 5.3073 cm∕min , u = 1.6445 cm∕min ,13

Author's personal copy
The Linear Boltzmann Equation in Column Experiments of Porous…
1
0.8

C/C0

Fig. 2  The relative concentration
C∕C0 for Run 1 (black dots) is
plotted with calculated values
from LBE (red solid line) andADE (blue dashed line)

0.6
0.4
Run 1
ADE
LBE

0.2
0

0

10

20

30

40

50

40

50

40

50

t (min)Fig. 3  Same as Fig. 2, but
1 − C∕C0 is plotted in the semilog graph

0

10

-1

10

-2

1-C/C0

10

10-3

Run 1
ADE
LBE

0

10

20

30

t (min)

1
0.8

C/C0

Fig. 4  The relative concentration
C∕C0 for Run 2 (black dots) is
plotted with calculated values
from LBE (red solid line) andADE (blue dashed line)

0.6
0.4
Run 2
ADE
LBE

0.2
0

0

10

20

30

t (min)𝛽 = 0.09130 for LBE, and u = 1.2886 cm∕min , D = 1.8379 cm2 ∕min for ADE. In Fig. 2,
experimental results agree well with both LBE and ADE. Figure 3 shows semilog plotsof 1 − C∕C0 as a function of t.

5.1.2 Run 2: Large Beads
The column is filled with large beads. The measured porosity is 𝜂 = 0.367. (This value ofthe porosity is smaller than that of Run 1. This may be due to small bubbles attached on

13Author's personal copy
K. Amagai et al.

the beads surface during Run 1.) The bed height is 18.5 cm with section area 6.61 cm2. The
flow rate is 0.875 cm3 ∕min.
Figure 4 for Run 2 is the same as Fig. 2, but large glass beads are used. By nonlinear least-squares fitting, the parameters are obtained as 𝜎a = 1.0 × 10−8 min−1,
𝜎s = 4.6778 min−1, v0 = 8.1753 cm∕min, u = 2.4999 cm∕min, 𝛽 = 0.08579 for LBE, andu = 1.6189 cm∕min, D = 3.8519 cm2 ∕min for ADE. In Fig. 4, experimental results aredescribed well by LBE and ADE. The longtime behavior is shown in Fig. 5, in which
1 − C∕C0 is plotted.

5.2 Column Experiments with Adsorbed Solute
Here, we use adsorptive solute (zinc solution). The filling material is the standard sand(Tohoku silica sand No. 4, Kitanihon Sangyo). The median diameter of the sand is 750 μm.
As a preparation, we eliminate the organic matter that may have been contained in the sandby soaking it in HNO3 solution. The zinc solution is 2 ppm. We set a filter on the top of thesand bed made with glass wool. We also put the same filter at the bottom of the column.
A short column of length 12.0 cm with internal diameter 3.1 cm was used. We perform
a blank test beforehand to make sure that the zinc is not absorbed on the surface of the column wall. The concentration is measured with the atomic absorption photometer (Z-2300,
Hitachi High-Technologies), the compressor (SC820, Koki Holdings), and the Neo Cool
Circulator (CF700, Yamato).

5.2.1 Run 3
The measured porosity is 0.295, and the flow rate is 9.34 cm3 ∕h. The bed height is 10.7 cm,
and the section area is 7.54 cm2.
Figure 6 shows the relative concentration C∕C0 of the zinc breakthrough curve for Run
3 (black dots) with calculated values from the linear Boltzmann equation (LBE) (red solidline) and the Ogata–Banks ADE model (blue dashed line). By nonlinear least-squares fitting, we obtain 𝜎a = 1.0 × 10−7 h−1, 𝜎s = 2.8134 h−1, v0 = 5.0663 cm∕h, u = 1.9876 cm∕h,𝛽 = 0.1739 for LBE, and u = 1.1406 cm∕h, D = 4.3864 cm2 ∕h. The dashed curve for ADE
partially deviates from the measured values. The semilog plot in Fig. 7 shows a clear discrepancy between the experimental result and ADE.

Fig. 5  Same as Fig. 4, but
1 − C∕C0 is plotted in the semilog graph
1-C/C0

10

0

10

-1

10

-2

10

-3

Run 2
ADE
LBE

0

10

20

30

t (min)

13

40

50

Author's personal copy
The Linear Boltzmann Equation in Column Experiments of Porous…
1
0.8

C/C0

Fig. 6  The relative concentration
C∕C0 for Run 3 (black dots) is
plotted with calculated values by
LBE (red solid line) and by ADE(blue dashed line)

0.6
0.4
Run 3
ADE
LBE

0.2
0

Fig. 7  Same as Fig. 6, but
1 − C∕C0 is plotted in the semilog graph
1-C/C0

10

5

10

t (h)

15

20

25

20

25

20

25

0

Run 3
ADE
LBE

10-1

10

-2

0

5

10

t (h)

15

1
0.8

C/C0

Fig. 8  The relative concentration
C∕C0 for Run 4 (black dots) is
plotted with calculated values by
LBE (red solid line) and by ADE(blue dashed line)

0

0.6
0.4
Run 4
ADE
LBE

0.2
0

0

5

10

t (h)

15

5.2.2 Run 4
Figure 8 for Run 4 is the same as Fig. 6, but we repeat another run for the standard sandand adsorbed solute. The measured porosity is 0.289, and the flow rate is 11.25 cm3 ∕h.
The bed height is 10.7 cm, and the section area is 7.54 cm2. By nonlinear least-squares fitting, we obtain 𝜎a = 1.0 × 10−7 h−1, 𝜎s = 2.6773 h−1, v0 = 4.1166 cm∕h, u = 1.7999 cm∕h,𝛽 = 0.2091 for LBE, and u = 1.1384 cm∕h, D = 3.9814 cm2 ∕h for ADE. The discrepancy
of ADE is even more apparent. The semilog plot in Fig. 9 shows a discrepancy between theexperimental result and ADE.

13

Author's personal copy
K. Amagai et al.
Fig. 9  Same as Fig. 8, but
1 − C∕C0 is plotted in the semilog graph
1-C/C0

10

0

Run 4
ADE
LBE

10-1

-2

10

0

5

10

t (h)

15

20

25

6 Conclusions
Not only for field observation reported in studies such as Adams and Gelhar 1992 and Benson et al. 2000, 2001, anomalous diffusion also appears in laboratory-scale experiments(Cortis et al. 2004). In Sect. 5, ADE does not precisely reproduce the breakthrough curves
even for beads, which form a medium material with non-adsorbed solution. Williams proposed to use the linear Boltzmann equation for transport in porous media (Williams 1992a,b, 1993a, b). It was shown that the proposed model of the linear Boltzmann equation reproduces the whole breakthrough curves of experimental results including the non-monotonic
decay of 1 − C∕C0 in the tail.
The diffusion approximation holds when the propagation distance of tracer particles issufficiently larger than the transport mean free path 𝓁 ∗, which is given in (14). In the caseof light transport, the photon diffusion equation starts to work when the propagation distance becomes ten times larger than the transport mean free path (Yoo et al. 1990). Using
the relation (13), we can calculate the diffusion coefficient, which is denoted by D′, fromparameters in the linear Boltzmann equation. We let D denote the fitted diffusion coefficient in the advection–diffusion equation. Let us compare Run 1 and Run 2. For Run 1
in Sect. 5.1.1, 𝓁 ∗ = 1.03 and D� = 1.818, whereas D = 1.8379. The relative difference
is |D − D� |∕D = 0.0108. For Run 2 in Sect. 5.1.2, 𝓁 ∗ = 1.75 and D� = 4.762, whereas
D = 3.8519. The relative difference is |D − D� |∕D = 0.236. In either case, the propagation
distance of 18 cm or 18.5 cm is much larger than 𝓁 ∗ and the transport is expected to bein the diffusion regime. Indeed, Figs. 2 and 4 show that both LBE and ADE reproduce
the experimentally obtained breakthrough curves. Next, let us also look at Run 3 and Run
4. For Run 3 in Sect. 5.2.1, 𝓁 ∗ = 1.80 and D� = 3.041, whereas D = 4.3864. The relative
difference is |D − D� |∕D = 0.307. For Run 4 in Sect. 5.2.2, 𝓁 ∗ = 1.54 and D� = 2.110,
whereas D = 3.9814. The relative difference is |D − D� |∕D = 0.470. The propagation distance of 10.7 cm is not as long as the distance in Run 1 and Run 2 in terms of the transportmean free path. The discrepancy of the curve for ADE in Figs. 6 and 8 might also attributeto anisotropic scattering for the sand. The transport mean free path becomes larger when
the scattering in LBE becomes anisotropic. Moreover, the discrepancy may come from thefact that the surface of the sand is rough and the particle speed cannot be modeled by aconstant v0.
Although we assumed isotropic scattering in (1), it is possible to take anisotropic scattering
1
into account and write the integral term in a more general form as 𝜎s ∫−1 p(𝜇, 𝜇� )𝜓(x, 𝜇 � , t) d𝜇�,
where p(𝜇, 𝜇� ) is the probability that tracer particles moving in the direction 𝜇′ change their

13

Author's personal copy
The Linear Boltzmann Equation in Column Experiments of Porous…

direction to the direction 𝜇 by scattering. We can similarly develop the analytical discrete ordi∑
nates method for anisotropic scattering by writing p(𝜇, 𝜇� ) = 12 l 𝛽l Pl (𝜇)Pl (𝜇 � ) with someconstant coefficients 𝛽l and Legendre polynomials Pl (l = 0, 1, …) (Barichello 2011; Siewert
2000).
In this paper, we have shown that the mass transport in column experiments, whichshow anomalous diffusion, obeys the linear Boltzmann equation. Furthermore, the linear
Boltzmann equation in the transport regime appears at the mesoscopic scale and the advection–diffusion equation in the diffusion regime is derived from the linear Boltzmann equation in the asymptotic limit when the propagation distance is sufficiently larger than 𝓁 ∗ at
the macroscopic scale. The physical origin of anomalous diffusion in a heterogeneous random medium with more complex structure of fissures is still an open problem.
Acknowledgements The seed of this research was planted on the occasion of the Study Group Workshop (Department of Mathematical Sciences, The University of Tokyo, December 2010), which is greatly
appreciated. The research was restarted with the support of the Focusing Collaborative Research by theInterdisciplinary Project on Environmental Transfer of Radionuclides (University of Tsukuba and Hirosaki
University). The research was partially supported by the JSPS A3 Foresight Program: Modeling and Computation of Applied Inverse Problems. MM acknowledges support from Grant-in-Aid for Scientific Research(17H02081, 17K05572, 18K03438) of the Japan Society for the Promotion of Science. YH acknowledges
support from Grant-in-Aid for Scientific Research (15H05740, 17H01478) of the Japan Society for the Promotion of Science.

Appendix 1: Eigenvalues
Let us write the homogeneous equation as
(
)
)
𝜇 (
1
𝜇t I − Ξ+ Φ+ = s W Φ+ + Φ− ,
𝜈
2
)
(
)
𝜇s (
1
𝜇t I − Ξ− Φ− = W Φ+ + Φ− ,
𝜈
2

(50)

where I is the identity, Ξ± and W are matrices, and Φ± are vectors defined as⎞⎛±𝜇1 + 𝜂⎟⎜± 𝜇2 + 𝜂
Ξ± = ⎜⎟,⋱⎟⎜± 𝜇N + 𝜂 ⎠⎝We obtain

SM
where

S=

( −1
Ξ+

( )
( )
Φ+
1 Φ+
=
,
Φ−
𝜈 Φ−

)
Ξ−1
−

,

⎛ 𝜙(𝜈, ±𝜇1 ) ⎞⎜ 𝜙(𝜈, ±𝜇2 ) ⎟Φ± = ⎜⎟,⋮⎟⎜⎝𝜙(𝜈, ±𝜇N )⎠M = 𝜇t (s)

( )
(
)
𝜇 W W
I
− s
.
I
2 W W

{W}ij = wj . (51)

(52)

(53)

We label the eigenvalues as 𝜈n (n = 1, 2, … , 2N ).
By deforming the Bromwich contour, we obtain a contour (−∞, 𝛾), 𝛾 > 0, on which s isreal. In this case, we have13

Author's personal copy
K. Amagai et al.

(54)

SMw0 = 𝜆0 w0 ,where 𝜆0 is real and w0 a 2N-dimensional real vector. We note that M = ℜM in this caseand moreover we can write

𝜇t (s) =𝜎a + 𝜎s + ℜs
.
v0

(55)

Since M is a symmetric real matrix, by the Cholesky decomposition we can write
M = U T U , where U is a real triangular matrix with positive diagonal entries. Hence,

USU T w
̃0 = 𝜆̃w0 ,w
̃0 = Uw0 .

(56)

Since U is nonsingular, S and USU T have the same inertia, i.e., the same number of positive, negative, and zero eigenvalues according to Sylvester’s law of inertia. Clearly, S has
N𝜂 positive eigenvalues. Therefore, there are N𝜂 eigenvalues 𝜆0.
Next, we assume that s has a small imaginary part. Then, we can write
(
( ))
ℑs I
S ℜM + iw = 𝜆w.
(57)
I
v0
We treat the imaginary part as perturbation and express the matrix-vector equation as
(
( ))
) (
)(
)
(
ℑs I
S ℜM + i𝜖
w0 + 𝜖w1 + ⋯ = 𝜆0 + 𝜆1 + ⋯ w0 + 𝜖w1 + ⋯ . (58)Iv0
By collecting terms of order O(𝜖 0 ), we have S(ℜM)w0 = 𝜆0 w0. From terms of O(𝜖 1 ), we
have

S(ℜM)w1 + iℑs
Sw = 𝜆0 w1 + 𝜆1 w0 .
v0 0

(59)

Let us multiply wT0 on both sides of the above equation. We obtainwT0 S(ℜM)w1 + iℑs T
w Sw = 𝜆0 wT0 w1 + 𝜆1 wT0 w0 .
v0 0 0(60)

Using wT0 (S(ℜM))T = wT0 SℜM = 𝜆0 wT0 , we obtain

𝜆1 = i(ℑs)wT0 Sw0v0 wT0 w0.

(61)

That is, 𝜆1 is pure imaginary and the number of positive ℜ𝜆 does not change. This fact
implies that there are N𝜂 eigenvalues 𝜈n (n = 1, … , 2N ) such that ℜ𝜈n > 0.

Appendix 2: Advection–Diffusion Equation with Absorption
Let us consider the following advection–diffusion equation with the absorption term:

13

Author's personal copy
The Linear Boltzmann Equation in Column Experiments of Porous…

(

)
𝜕2
𝜕
𝜕
− D 2 + u + 𝜎a C(x, t) = 0,
𝜕t
𝜕x
𝜕x

(62)

C(x, 0) = 0,(63)

C(0, t) = C0 .

(64)

Let us introduce Γ(x, t) as(

C(x, t) = Γ(x, t) exp

)
u2 t
ux
−
− 𝜎a t .
2D 4D

(65)

We have𝜕2𝜕Γ(x, t) − D 2 Γ(x, t) = 0,
𝜕t
𝜕x

(66)

Γ(x, 0) = 0,

(67)

( 2
)
u t
+ 𝜎a t .
Γ(0, t) = C0 exp
4D

(68)

Let us consider the Laplace transform:̂ p) =
Γ(x,∫0∞e−pt Γ(x, t) dt,

u2+ 𝜎a .
4D

(69)

C0
.
u2
p − 4D − 𝜎a

(70)

ℜp >

Then, we obtainp̂
𝜕2 ̂
Γ(x, p) = Γ(x,
p),
D
𝜕x2

̂ p) =
Γ(0,We obtain√pC0̂ p) =
Γ(x,

−

x

.

(71)

f (t − 𝜏)g(𝜏) d𝜏,

(72)

u2
− 𝜎a
p − 4D

e

D

Thus,

Γ(x, t) = C0

∫0

t

where

f̂ (p) =

1
2

u
− 𝜎a
p − 4D√p,ĝ (p) = e

−

D

x

.

(73)

Noting that13

Author's personal copy
K. Amagai et al.

∫0∞e−pt eat dt =1(ℜp > a),
p−a

∫0

∞

√
b2
b
e−pt √
e− 4t dt = e−b p (ℜp > 0),2 𝜋t3∕2(74)

we obtain√��x∕ D − x2u2
e 4D𝜏 d𝜏+ 𝜎a (t − 𝜏) √∫0
4D
2 𝜋𝜏 3∕2
�
� �
√
�� 2
��
x∕ D
u
u2
x
= √ exp+ 𝜎a t exp −+ 𝜎a √
4D
4D
2 𝜋
D
2
��
⎛ ⎛ √
� ⎞ ⎞
t
⎜ ⎜ x∕ D⎟u2
1+ 𝜎a 𝜏 ⎟ ⎟ d𝜏
exp ⎜− √ −×∫0 𝜏 3∕2⎜⎟
4D⎜ ⎝ 2 𝜏
⎠ ⎟⎠
⎝
�
�
(75)
�
�� 2
��
C0
u
u2
x
=
exp
+ 𝜎a t exp −+ 𝜎a √
2
4D
4D
D
��
�⎞
�
⎛
u2
x
+ 𝜎a t ⎟−× erfc ⎜ √
⎟
⎜ 4Dt
4D
⎠
⎝
��
��
�
� ⎞�
⎛
u2
u2
2x
x
⎜
+ 𝜎a √+ 𝜎a t ⎟ ,++ exp
erfc √
⎟
⎜ 4Dt
4D
4D
D
⎠
⎝
√
2
∞
where the complementary error function is given by erfc(x) = (2∕ 𝜋) ∫x e−t dt . Finally,
� �
�
2
C0 ux
u
x
e 2D exp −+ 𝜎a √
C(x, t) =
2
4D
D
��
�⎞
�
⎛
u2
x
+ 𝜎a t ⎟−× erfc ⎜ √
(76)
⎟
⎜ 4Dt
4D
⎠
⎝
��
��
�
� ⎞�
⎛
x
u2
2x
u2
+ 𝜎a √+ 𝜎a t ⎟ .
++ exp
erfc ⎜ √
⎜ 4Dt
⎟
4D
4D
D
⎝
⎠
��

t

Γ(x, t) = C0

exp

The above solution reduces to the Ogata–Banks solution (Ogata and Banks 1961) when𝜎a = 0. In particular, we have C(x, t) → C0 as t → ∞ if 𝜎a = 0. When 𝜎a = 0, we have
�
��
�
�
�
C0
x − utx+
ut
erfc √+ eux∕D erfc √.
C(x, t) =
(77)
2
4Dt
4Dt

13

Author's personal copy
The Linear Boltzmann Equation in Column Experiments of Porous…

References
Adams, E.E., Gelhar, L.W.: Field study of dispersion in a heterogeneous aquifer, 2, spatial moments
analysis. Water Resour. Res. 28(12), 3293–3307 (1992). https​://doi.org/10.1029/92WR0​1757
Apresyan, L.A., Kravtsov, Y.A.: Radiation Transfer: Statistical and Wave Aspects. Gordon and Breach,
Amsterdam (1996)
Barichello, L.B.: Explicit formulations for radiative transfer problems. In: Orlande, H.R.B., Fudym, O.,
Maillet, D., Cotta, R.M. (eds.) Thermal Measurements and Inverse Techniques. CRS Press, Boca
Raton (2011)
Barichello, L.B., Siewert, C.E.: A new version of the discrete-ordinates method. In: Proceedings of the
2nd International Conference on Computational Heat and Mass Transfer, Rio de Janeiro, pp. 22–26(2001)
Barichello, L.B., Garcia, R.D.M., Siewert, C.E.: Particular solutions for the discrete-ordinates method.
J. Quant. Spectrosc. Radiat. Trans. 64(3), 219–226 (2000). https​://doi.org/10.1016/S0022​-4073(98)00146​-0
Benson, S.M., Orr Jr., F.M.: Carbon dioxide capture and storage. Harnessing Mater. Energy 33(4), 303–305 (2008). https​://doi.org/10.1557/mrs20​08.63
Benson, D.A., Wheatcraft, S.W., Meerschaert, M.M.: Application of a fractional advection–dispersion
equation. Water Resour. Res. 36(6), 1403–1412 (2000). https​://doi.org/10.1029/2000W​R9000​31
Benson, D.A., Schumer, R., Meerschaert, M.M., Wheatcraft, S.W.: Fractional dispersion, Lévy
motion, and the MADE tracer tests. Transp. Porous Media 42(1–2), 211–240 (2001). https​://doi.
org/10.1023/A:10067​33002​131
Berkowitz, B., Scher, H.: Theory of anomalous chemical transport in random fracture networks. Phys.
Rev. E 57(5), 5858–5869 (1998). https​://doi.org/10.1103/PhysR​evE.57.5858
Berkowitz, B., Scher, H., Silliman, S.E.: Anomalous transport in laboratory-scale, heterogeneous porous
media. Water Resour. Res. 36(1), 149–158 (2000). https​://doi.org/10.1029/1999W​R9002​95
Case, K.M.: Elementary solutions of the transport equation and their applications. Ann. Phys. 9(1), 1–23
(1960). https​://doi.org/10.1016/0003-4916(60)90060​-9Case, K.M., Zweifel, P.F.: Linear Transport Theory. Addison-Wesley, Boston (1967)
Chakraborty, P., Meerschaert, M.M., Lim, C.Y.: Parameter estimation for fractional transport: a particle-tracking approach. Water Resour. Res. 45(10), W10415 (2009). https​://doi.org/10.1029/2008W​
R0075​77
Chandrasekhar, S.: Radiative Transfer. Dover, New York (1960)
Cortis, A., Chen, Y., Scher, H., Berkowitz, B.: Quantitative characterization of pore-scale disordereffects on transport in “homogeneous” granular media. Phys. Rev. E 70(4), 041108 (2004). https​://doi.org/10.1103/PhysR​evE.70.04110​8
De Anna, P., Le Borgne, T., Dentz, M., Tartakovsky, A.M., Bolster, D., Davy, P.: Flow intermittency,
dispersion, and correlated continuous time random walks in porous media. Phys. Rev. Lett. 110(18),184502 (2013). https​://doi.org/10.1103/PhysR​evLet​t.110.18450​2Duderstadt, J.J., Martin, W.R.: Transport Theory. Wiley, New York (1979)
Furutsu, K., Yamada, Y.: Diffusion approximation for a dissipative random medium and the applications.
Phys. Rev. E 50(5), 3634–3640 (1994). https​://doi.org/10.1103/PhysR​evE.50.3634
Gelhar, L.W.: Stochastic subsurface hydrology from theory to applications. Water Resour. Res. 22(9),135S–145S (1986). https​://doi.org/10.1029/WR022​i09Sp​0135S​
Golub, G.H., Welsch, J.H.: Calculation of Gauss quadrature rules. Math. Comput. 23(106), 221–230(1969). https​://doi.org/10.2307/20044​18
Hatano, Y., Hatano, N.: Dispersive transport of ions in column experiments: an explanation of longtailed profiles. Water Resour. Res. 34(5), 1027–1033 (1998). https​://doi.org/10.1029/98WR0​0214
Ishimaru, A.: Wave Propagation and Scattering in Random Media. Academic, San Diego (1978)
Kelly, J.F., Bolster, D., Meerschaert, M.M., Drummond, J.D., Packman, A.I.: FracFit: a robust parameter
estimation tool for fractional calculus models. Water Resour. Res. 53(3), 2559–2567 (2017). https​://doi.org/10.1002/2016W​R0197​48
Kennedy, C.A., Lennox, W.C.: A stochastic interpretation of the tailing effect in solute transport. Stoch.
Environ. Res. Risk Assess. 15(4), 325–340 (2001). https​://doi.org/10.1007/s0047​70100​076
Levy, M., Berkowitz, B.: Measurement and analysis of non-Fickian dispersion in heterogeneous porous
media. J. Contam. Hydrol. 64(3–4), 203–226 (2003). https​://doi.org/10.1016/S0169​-7722(02)00204​-8
Liang, Y., Chen, W., Xu, W., Sun, H.: Distributed order Hausdorff derivative diffusion model to characterize non-Fickian diffusion in porous media. Commun. Nonlinear Sci. Numer. Sim. 70, 384–393 (2019).
https​://doi.org/10.1016/j.cnsns​.2018.10.010

13Author's personal copy
K. Amagai et al.
Liu, H., Kang, Q., Leonardi, C.R., Schmieschek, S., Narvaez, A., Jones, B.D., Williams, J.R., Valocchi,
A.J., Harting, J.: Multiphase lattice Boltzmann simulations for porous media applications. Comput.
Geosci. 20(4), 777–805 (2016). https​://doi.org/10.1007/s1059​6-015-9542-3
Matsuda, N., Mikami, S., Shimoura, S., Takahashi, J., Nakano, M., Shimada, K., Uno, K., Hagiwara, S.,
Saito, K.: Depth profiles of radioactive cesium in soil using a scraper plate over a wide area surrounding the Fukushima Dai-ichi Nuclear Power Plant, Japan. J. Environ. Radioact. 139, 427–434 (2015).
https​://doi.org/10.1016/j.jenvr​ad.2014.10.001
Meerschaert, M.M., Zhang, Y., Baeumer, B.: Tempered anomalous diffusion in heterogeneous systems.
Geophys. Res. Lett. 35(17), L17403 (2008). https​://doi.org/10.1029/2008G​L0348​99
Metzler, R., Klafter, J.: The random walk’s guide to anomalous diffusion: a fractional dynamics approach.
Phy. Rep. 339(1), 1–77 (2000). https​://doi.org/10.1016/S0370​-1573(00)00070​-3Montero, M., Masoliver, J.: Nonindependent continuous-time random walks. Phys. Rev. E 76(6), 061115(2007). https​://doi.org/10.1103/PhysR​evE.76.06111​5
Montroll, E.W., Weiss, G.H.: Random walks on lattices. II. J. Math. Phys. 6(2), 1672013181 (1965). https​://doi.org/10.1063/1.17042​69
Moroni, M., Cushman, J.H., Cenedese, A.: Application of photogrammetric 3D-PTV technique to track
particles in porous media. Transp. Porous Media 79(1), 43–65 (2009). https​://doi.org/10.1007/s1124​
2-008-9270-4
Nissan, A., Dror, I., Berkowitz, B.: Time-dependent velocity-field controls on anomalous chemical transportin porous media. Water Resour. Res. 53(5), 3760–3769 (2017). https​://doi.org/10.1002/2016W​R0201​
43
Ogata, A., Banks, R.B.: A solution of the differential equation of longitudinal dispersion in porous media.
Professional paper, 411-A (1961). https​://doi.org/10.3133/pp411​A
Rubin, Y.: Applied Stochastic Hydrogeology. Oxford University Press, London (2003)
Schumer, R., Benson, D.A., Meerschaert, M.M., Baeumer, B.: Fractal mobile/immobile solute transport.
Water Resour. Res. 39(10), 1296 (2003). https​://doi.org/10.1029/2003W​R0021​41
Sen, D., Nobes, D.S., Mitra, S.K.: Optical measurement of pore scale velocity field inside microporous
media. Microfluid. Nanofluid. 12(1–4), 189–200 (2012). https​://doi.org/10.1007/s1040​4-011-0862-x
Siewert, C.E.: A concise and accurate solution to Chandrasekhar’s basic problem in radiative transfer. J.
Quant. Spectrosc. Radiat. Trans. 64(2), 109–130 (2000). https​://doi.org/10.1016/S0022​-4073(98)00144​-7
Siewert, C.E., Wright, S.J.: Efficient eigenvalue calculations in radiative transfer. J. Quant. Spectrosc.
Radiat. Trans. 62(6), 685–688 (1999). https​://doi.org/10.1016/S0022​-4073(98)00099​-5
Sobolev, V.V.: Light Scattering in Planetary Atmospheres. Pergamon, Oxford (1976)
Succi, S.: The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Clarendon, Oxford (2001)
Sun, H., Chen, W., Chen, Y.: Variable-order fractional differential operators in anomalous diffusion
modeling. Phys. A Stat. Mech. Appl. 388(21), 4586–4592 (2009). https​://doi.org/10.1016/j.physa​.2009.07.024
Sun, H., Li, Z., Zhang, Y., Chen, W.: Fractional and fractal derivative models for transient anomalous diffusion: model comparison. Chaos Solitons Fractals 102, 346–353 (2017). https​://doi.org/10.1016/j.chaos​.2017.03.060
Thomas, S.: Enhanced oil recovery—an overview. Oil Gas Sci. Technol. 63(1), 9–19 (2008). https​://doi.
org/10.2516/ogst:20070​60
Thomas, G.E., Stamnes, K.: Radiative Transfer in the Atmosphere and Ocean. Cambridge University Press,
New York (1999)
Van Genuchten, M.T., Wierenga, P.J.: Mass transfer studies in sorbing porous media I. Analytical solutions.
Soil Sci. Soc. Am. J. 40(4), 473–480 (1976). https​://doi.org/10.2136/sssaj​1976.03615​99500​40000​
40011​x
Wei, S., Chen, W., Hon, Y.C.: Characterizing time dependent anomalous diffusion process: a survey on fractional derivative and nonlinear models. Phys. A Stat. Mech. Appl. 462(15), 1244–1251 (2016). https​://doi.org/10.1016/j.physa​.2016.06.145
Weideman, J., Trefethen, L.: Parabolic and hyperbolic contours for computing the Bromwich integral. Math.
Comput. 76(259), 1341–1356 (2007). https​://doi.org/10.1090/S0025​-5718-07-01945​-X
Williams, M.M.R.: Stochastic problems in the transport of radioactive nuclides in fractured rock. Nuclear
Sci. Eng. 112(3), 215–230 (1992a). https​://doi.org/10.13182​/NSE92​-A2907​0Williams, M.M.R.: A new model for describing the transport of radionuclides through fractured rock. Ann.
Nuclear Energy 19(10–12), 791–824 (1992b). https​://doi.org/10.1016/0306-4549(92)90018​-7
Williams, M.M.R.: A new model for describing the transport of radionuclides through fractured rock. Part
II: numerical results. Ann. Nuclear Energy 20(3), 185–202 (1993a). https​://doi.org/10.1016/03064549(93)90101​-T

13

Author's personal copy
The Linear Boltzmann Equation in Column Experiments of Porous…
Williams, M.M.R.: Radionuclide transport in fractured rock a new model: application and discussion. Ann.
Nuclear Energy 20(4), 279–297 (1993b). https​://doi.org/10.1016/0306-4549(93)90083​-2
Yoo, K.M., Liu, F., Alfano, R.R.: When does the diffusion approximation fail to describe photon transport in random media? Phys. Rev. Lett. 64(22), 2647–2650 (1990). https​://doi.org/10.1103/PhysR​evLet​
t.64.2647
Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps andinstitutional affiliations.

13