The linear Boltzmann equation in column experiments of porous media

TL;DR

This paper demonstrates that the linear Boltzmann equation effectively models transport phenomena in porous media column experiments, and derives the advection-diffusion equation as a limiting case.

Contribution

It introduces the linear Boltzmann equation as a new modeling approach for porous media transport and connects it to classical advection-diffusion models.

Findings

Breakthrough curves are accurately reproduced by the linear Boltzmann equation.

The advection-diffusion equation is derived as an asymptotic limit.

The approach bridges microscopic transport and macroscopic modeling.

Abstract

The use of the linear Boltzmann equation is proposed for transport in porous media in a column. 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.

Equations

$$\frac{\partial }{\partial t}\psi (x,\mu ,t)+(u+v_0\mu )\frac{\partial }{\partial x}\psi +({\sigma }_a+{\sigma }_s)\psi =\frac{{\sigma }_s}{2}{\int }_{-1}^1\psi (x,\mu ,t)\mathrm{d}\mu .$$

$$\psi (x,\mu ,0)=0,$$

$$\psi (0,\mu ,t)=n_0\delta (\mu -1),\mu >-\eta ,$$

$$\eta =\frac{u}{v_0}.$$

$$n(x,t)={\int }_{-1}^1\psi (x,\mu ,t)\mathrm{d}\mu .$$

$$C(x,t)=\alpha n(x,t),$$

$$\psi (x,\mu ,t)=\frac{1}{2\alpha }C(x,t)+\frac{3}{2\alpha v_0}J(x,t)\mu ,$$

$$C(x,t)$$

$$J(x,t)$$

$$\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+\frac{\partial J}{\partial x}+{\sigma }_{a}C=0.$$

$$\frac{1}{v_0}\frac{\partial J}{\partial t}+\frac{u}{v_0}\frac{\partial J}{\partial x}+\frac{v_0}{3}\frac{\partial C}{\partial x}+\frac{{\sigma }_a+{\sigma }_s}{v_0}J=0.$$

$$J=-D\frac{\partial C}{\partial x},$$

$$D$$

$${\ell }^{\ast }$$

$$\frac{\partial C}{\partial t}-D\frac{{\partial }^{2}C}{\partial x^2}+u\frac{\partial C}{\partial x}+{\sigma }_{a}C=0,$$

$$C(0,t)-\frac{3D\eta (1-\eta )}{2}\frac{D}{u}\frac{\partial C}{\partial x}(0,t)=\frac{2\alpha n_0}{1+\eta }.$$

$$C(0,t)=C_0,$$

$$\beta =\frac{\alpha n_0}{C_0}.$$

$$\hat{\psi }(x,\mu ,p)={\int }_0^{\infty }{\mathrm{e}}^{-pt}\psi (x,\mu ,t)\mathrm{d}t.$$

$${\mu }_t=\frac{{\sigma }_a+{\sigma }_s+p}{v_0},{\mu }_s=\frac{{\sigma }_s}{v_0}.$$

$$(\eta +\mu )\frac{\partial }{\partial x}\hat{\psi }(x,\mu ,p)+{\mu }_t\hat{\psi }(x,\mu ,p)=\frac{{\mu }_s}{2}{\int }_{-1}^1\hat{\psi }(x,\mu ,p)\mathrm{d}\mu ,$$

$$\hat{\psi }(0,\mu ,p)=\frac{n_0}{p}\delta (\mu -1),\mu >-\eta .$$

$$\hat{\psi }(x,\mu ,p)=\widehat{{\psi }_b}(x,\mu ,p)+\widehat{{\psi }_s}(x,\mu ,p).$$

$$(\eta +\mu )\frac{\partial }{\partial x}\widehat{{\psi }_b}+{\mu }_t\widehat{{\psi }_b}=0,x>0,-1\le \mu \le 1,$$

$$\widehat{{\psi }_b}(0,\mu ,p)=\frac{n_0}{p}\delta (\mu -1),\mu >-\eta ,$$

$$(\eta +\mu )\frac{\partial }{\partial x}\widehat{{\psi }_s}+{\mu }_t\widehat{{\psi }_s}=\frac{{\mu }_s}{2}{\int }_{-1}^1\widehat{{\psi }_s}\mathrm{d}\mu +q,x>0,-1\le \mu \le 1,$$

$$\widehat{{\psi }_s}(0,\mu ,p)=0,\mu >-\eta ,$$

$$q(x,\mu ,p)=\frac{{\mu }_s}{2}{\int }_{-1}^1\widehat{{\psi }_b}(x,\mu ,p)\mathrm{d}\mu =\frac{n_0{\mu }_s}{2p}{\mathrm{e}}^{-x{\mu }_t/(\eta +\mu )}.$$

$$\widehat{{\psi }_b}(x,\mu ,p)=\frac{n_0}{p}{\mathrm{e}}^{-x{\mu }_t/(\eta +\mu )}\delta (\mu -1).$$

$$(\eta +{\mu }_i)\frac{\partial }{\partial x}\widehat{{\psi }_s}(x,{\mu }_i)+{\mu }_t\widehat{{\psi }_s}(x,{\mu }_i)$$

$$=\frac{{\mu }_s}{2}\sum _{j=1}^N{w}_j\left[\widehat{{\psi }_s}(x,{\mu }_j)+\widehat{{\psi }_s}(x,-{\mu }_j)\right]+q(x,{\mu }_i),$$

$$N=30.$$

Full text

1 Introduction

Mass transport in porous media is one of the central topics in hydrology. For example, mass transport 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 mass transport in porous media, is an urgent issue (Matsuda et al. 2015). The phenomenon of mass transport 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 from measured 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’s law (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 deviation from ADE was investigated by a random particle simulation (Kennedy and Lennox 2001).

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 waiting time and jump distance are suitably adjusted. Although CTRW is widely used in the study of 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 Masoliver 2007). 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; Sun et 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 reproduced by 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 the porous medium.

Although CTRW and fractional equations have achieved certain success, the physical meaning 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 Boltzmann transport picture of the transport in porous media. Tracer particles in a porous medium flow along fissures and fractures. We regard the process of a tracer particle changing its direction as a scattering. Thus for the ensemble 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 a certain probability distribution.

The rest of the paper is organized as follows. In Sec. 2, the linear Boltzmann equation is introduced. We explore the advection-diffusion approximation in Sec. 3. In Sec. 4, we solve the linear Boltzmann equation. Section 5 is devoted to column experiments. Conclusions are given in Sec. 6.

2 Linear Boltzmann transport

The linear Boltzmann equation, radiative transport equation, or transport equation governs the transport 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 intersteller medium (Apresyan and Kravtsov 1996; Chandrasekhar 1960; Ishimaru 1978; Sovolev 1976; Thomas and Stamnes 1999). Let us regard the transport of 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 in beads 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 $v_{0}/\sigma_{s}$ (see (1) for $v_{0},\sigma_{s}$). Indeed, the use of the linear Boltzmann equation for the flow of tracer particles in porous media was proposed by Williams (Williams 1992a, b, 1993a, 1993b).

A remark is necessary for our Boltzmann transport. The (nonlinear) Boltzmann equation can be simulated by the lattice Boltzmann scheme (Succi 2001), and the lattice Boltzmann method has become one of commonly used numerical methods to investigate the flow in porous media. Instead of fluid in porous media, we focus on the Boltzmann transport of tracer particles in the fluid. Hence, the proposed equation (1) below is a linear equation.

Having in mind the transport in column experiments which is described in Sec. 5, we consider the one-dimensional linear Boltzmann equation. Let $\mu\in[-1,1]$ be the cosine of the polar angle and $v_{0}>0$ be the inherent particle speed. Although the coefficient of the spatial derivative is $v_{0}\mu$ in Williams’ study (Williams 1992a, b) , in this paper, we take into account advection, which will be denoted by $u\geq 0$. Hence the velocity $v$ is written as $v=u+v_{0}\mu$. Let $\psi(x,\mu,t)$ be the angularly dependent number density of particles at position $x\geq 0$ in direction $\mu$ at time $t\geq 0$. We can write the linear Boltzmann equation as follows.

[TABLE]

where isotropic scattering is assumed in the integral on the right-hand side. Here, $\sigma_{a}\geq 0$ and $\sigma_{s}>0$ are the absorption and scattering coefficients, respectively. Let $n_{0}$ 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

[TABLE]

where $\delta(\mu-1)$ is the Dirac delta function and

[TABLE]

We note that unlike the usual half-range boundary condition of $\mu>0$ (Case and Zweifel 1967), the boundary value in (3) is specified for $\mu\in(-\eta,1]$ due to the presence of the advection. We regard the column as the half space ($0<x<\infty$) and have $\psi\to 0$ as $x\to\infty$.

The particle number density $n(x,t)$ is calculated as

[TABLE]

Then the concentration is given by

[TABLE]

where $\alpha$ is the mass per particle.

3 Advection diffusion approximation

Let us derive the advection-diffusion equation from the linear Boltzmann equation. According the usual procedure of the diffusion approximation (Duderstadt and Martin 1979; Ishimaru 1978), we assume that the angular dependence of $\psi(x,\mu,t)$ is weak and written as

[TABLE]

where,

[TABLE]

In (7), we implicitly assumed a large propagation distance and long time. If absorption is strong, $\psi$ becomes zero before the form (7) is achieved. Hence we assume $\sigma_{a}$ is small. Indeed in Sec. 5, we will see that $\sigma_{a}$ is negligibly small. Indeed, the approximation in (7) is called the $P_{1}$ approximation (Case and Zweifel 1967). Below, we will proceed further and obtain the diffusion equation. We substitute the form (7) for $\psi$ in (1). By integrating the resulting equation over $\mu$, we obtain

[TABLE]

By integrating the equation over $\mu$ after multiplying $\mu$, we obtain

[TABLE]

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

[TABLE]

where

[TABLE]

Hence we obtain the following advection diffusion equation by substituting (12) into (10).

[TABLE]

For the photon diffusion equation, $\ell^{*}$ 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 $\ell^{*}$.

The diffusion equation (15) reduces to the usual advection-diffusion equation when $\sigma_{a}=0$. Since we found $\sigma_{a}\approx 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 B, the solution to (15) is presented.

Finally, we consider the boundary condition. Let us substitute (7) for $\psi(x,\mu,t)$ in the boundary condition (3). By integrating both sides of the resulting equation over $\mu$ from $-\eta$ to $1$, we obtain

[TABLE]

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

[TABLE]

where $C_{0}>0$ denotes the boundary source. Since the diffusion approximation does not hold on the boundary, the actual proportionality constant between $\alpha n_{0}$ and $C_{0}$ is a fitting parameter. We introduce

[TABLE]

4 Concentration from the linear Boltzmann equation

Let us consider the Laplace transform.

[TABLE]

We introduce

[TABLE]

Then our transport equation is written as

[TABLE]

with the boundary condition,

[TABLE]

We will compute $\hat{\psi}$ with the analytical discrete ordinates method (ADO) (Siewert and Wright 1999; Barichello et al. 2000; Barichello and Siewert 2001). Unlike the standard ADO, the coefficient $\mu_{t}$ is a complex number since the variable $p$ is complex and moreover, the boundary condition is given for $\mu>-\eta$ instead of $\mu>0$. We begin by decomposing $\hat{\psi}$ into the ballistic term $\hat{\psi}_{b}$ and scattering term $\hat{\psi}_{s}$ as

[TABLE]

Here, $\hat{\psi}_{b}$ obeys

[TABLE]

with the boundary condition

[TABLE]

and $\hat{\psi}_{s}$ obeys

[TABLE]

with the boundary condition

[TABLE]

where

[TABLE]

We obtain

[TABLE]

Let us drop $p$ and write $\hat{\psi}_{s}(x,\mu)=\hat{\psi}_{s}(x,\mu,p)$ and $q(x,\mu)=q(x,\mu,p)$. For the computation of $\hat{\psi}_{s}$, we discretize the integral by the Gauss-Legendre quadrature and obtain

[TABLE]

where $\mu_{i},w_{i}$ ($i=1,2,\dots,2N$) are abscissas and weights, respectively. These $\mu_{i},w_{i}$ with $0<\mu_{1}<\cdots<\mu_{N}<1$ and $\mu_{N+i}=-\mu_{i}$ ($i=1,\dots,N$) can be calculated by the Golub-Welsch algorithm (Golub and Welsch 1969). In the numerical calculation in Sec. 5, we set

[TABLE]

Furthermore, we introduce $N_{\eta}$ as the largest integer such that $-\eta<\mu_{N_{\eta}}$. The scattering part $\hat{\psi}_{s}$ is obtained as

[TABLE]

where the Green’s function defined for each $p$ satisfies

[TABLE]

where $\delta_{ii_{0}}$ is the Kronecker delta. Here the boundary condition is given by

[TABLE]

Let us consider the following homogeneous equation to calculate the Green’s function with ADO (Siewert and Wright 1999; Barichello et al. 2000; Barichello and Siewert 2001). Note that the equation will be solved for each $p$.

[TABLE]

We note that $\hat{\psi}$ depends on $p$ through $\mu_{t}$. With separation of variables, we can write $\hat{\psi}$ as

[TABLE]

where $\nu$ is the separation constant. The function $\phi(\nu,\mu_{i})$ satisfies the normalization condition,

[TABLE]

For the integral equation in which the sum on the right-hand side of (35) is replaced by the integral $\frac{\mu_{s}}{2}\int_{-1}^{1}\hat{\psi}\,\mathrm{d}\mu$, this $\phi$ is called the singular eigenfunction (Case 1960). However, by the discretization, we obtain

[TABLE]

assuming $\nu\neq(\mu_{i}+\eta)/\mu_{t}$. It is known that $\mu\neq\mu_{i}/\mu_{t}$ if $\eta=0$ and $p$ is real (Siewert and Wright 1999). We can show the following orthogonality relation.

[TABLE]

where

[TABLE]

We can find $2N$ eigenvalues $\nu=\nu_{n}$ ($n=1,2,\dots,2N$). Moreover there are $N_{\eta}$ eigenvalues with positive real parts. See Appendix A for the eigenvalues.

Taking the fact that $G$ vanishes as $x\to\infty$, we can write

[TABLE]

where coefficients $B(\nu;x_{0},\mu_{i_{0}})$ will be determined later. The free-space Green’s function $G_{0}$ is calculated as

[TABLE]

where upper signs are chosen for $x>x_{0}$ and lower signs are used for $x<x_{0}$. Since $G$ vanishes at $x=0$ for $\mu_{i}>\mu_{N_{\eta}}$, we have from (4),

[TABLE]

Let us multiply (43) by $\exp(-x_{0}\mu_{t}/(\eta+1))$, then integrate both sides with $x_{0}$ and take the sum with respect to $i$. We obtain

[TABLE]

where

[TABLE]

We can numerically obtain $E(\nu_{n})$ from the linear system (44). Finally, we obtain

[TABLE]

In this way, $\hat{\psi}_{s}$ is computed using ADO.

The Laplace transform of the particle number density $\hat{n}(x,p)$ is obtained as

[TABLE]

We note that the particle number density $n(x,t)$ is given by

[TABLE]

The Bromwich integral in the inverse Laplace transform is numerically evaluated by the trapezoidal rule. In general, deformation of the contour can be considered for the inversion (Weideman and Trefethen 2007). We found $\gamma=0.04$ is suitable. Hence, the concentration is obtained as

[TABLE]

5 Column experiments

Column experiments are often used for the study of transport in random media (Cortis 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 or more until steady flow field is achieved. Next, tracer solution is injected to displace the fresh water. The discharged solution is collected by the fraction collector (CHF161RA, Advantec) at regular intervals. Column experiments are performed at room temperature of $25$–$26\,{}^{\circ}{\rm C}$.

Tracer experiments are conducted until the discharge concentration becomes equal to the influent concentration. For the sake of the accuracy of measurements, we correct the measurement 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 with Soda glass; (Run 1) smaller beads with diameter ranging from $0.177\,{\rm mm}$ to $0.250\,{\rm mm}$ and (Run 2) larger beads with diameter from $0.500\,{\rm mm}$ to $0.710\,{\rm mm}$. The column length is $20.0\,{\rm cm}$ with the internal diameter of the column $2.9\,{\rm cm}$. In order to prevent entrainment of 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 stirrer beads while filling 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\,{\rm\mu 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 $\pm 2\%$. The temperature compensation circuit is installed.

We prepare the NaCl solution of $0.1\%$ with the salt of purity $99.5\%$ (Wako). Injection of 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 $\eta=0.386$. The bed height is $18.0\,{\rm cm}$ with section area $6.61\,{\rm cm}^{2}$. The flow rate is $0.958\,{\rm cm}^{3}/{\rm min}$.

Figure 2 shows the chloride breakthrough curve from Run 1. The vertical axis shows the relative concentration $C/C_{0}$ and the horizontal axis is time $t$. The measured values from the column experiment (black dots) are compared to $C/C_{0}$ 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 squares of the residuals between the experimental and theoretical curves); they are obtained as $\sigma_{a}=1.0\times 10^{-8}\,{\rm min}^{-1}$, $\sigma_{s}=5.1645\,{\rm min}^{-1}$, $v_{0}=5.3073\,{\rm cm}/{\rm min}$, $u=1.6445\,{\rm cm}/{\rm min}$, $\beta=0.09130$ for LBE, and $u=1.2886\,{\rm cm}/{\rm min}$, $D=1.8379\,{\rm cm}^{2}/{\rm min}$ for ADE. In Fig. 2, experimental results agree well with both LBE and ADE. Figure 3 shows semi-log plots of $1-C/C_{0}$ as a function of $t$.

5.1.2 Run 2: large beads

The column is filled with large beads. The measured porosity is $\eta=0.367$. (This value of the porosity is smaller than that of Run 1. This may be due to small bubbles attached on the beads surface during Run 1.) The bed height is $18.5\,{\rm cm}$ with section area $6.61\,{\rm cm}^{2}$. The flow rate is $0.875\,{\rm cm}^{3}/{\rm 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 $\sigma_{a}=1.0\times 10^{-8}\,{\rm min}^{-1}$, $\sigma_{s}=4.6778\,{\rm min}^{-1}$, $v_{0}=8.1753\,{\rm cm}/{\rm min}$, $u=2.4999\,{\rm cm}/{\rm min}$, $\beta=0.08579$ for LBE, and $u=1.6189\,{\rm cm}/{\rm min}$, $D=3.8519\,{\rm cm}^{2}/{\rm min}$ for ADE. In Fig. 4, experimental results are described well by LBE and ADE. The long-time behavior is shown in Fig. 5, in which $1-C/C_{0}$ 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$\mu$m. As a preparation, we eliminate the organic matter that may have been contained in the sand by soaking it in ${\rm HNO}_{3}$ solution. The zinc solution is 2 ppm. We set a filter on the top of the sand bed made with glass wool. We also put the same filter at the bottom of the column.

A short column of length $12.0\,{\rm cm}$ with internal diameter $3.1\,{\rm 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\,{\rm cm}^{3}/{\rm h}$. The bed height is $10.7\,{\rm cm}$ and the section area is $7.54\,{\rm cm}^{2}$.

Figure 6 shows the relative concentration $C/C_{0}$ of the zinc breakthrough curve for Run 3 (black dots) with calculated values from the linear Boltzmann equation (LBE) (red solid line) and the Ogata-Banks ADE model (blue dashed line). By nonlinear least-squares fitting, we obtain $\sigma_{a}=1.0\times 10^{-7}\,{\rm h}^{-1}$, $\sigma_{s}=2.8134\,{\rm h}^{-1}$, $v_{0}=5.0663\,{\rm cm}/{\rm h}$, $u=1.9876\,{\rm cm}/{\rm h}$, $\beta=0.1739$ for LBE, and $u=1.1406\,{\rm cm}/{\rm h}$, $D=4.3864\,{\rm cm}^{2}/{\rm h}$. The dashed curve for ADE partially deviates from the measured values. The semi-log plot in Fig. 7 shows a clear discrepancy between the experimental result and ADE.

5.2.2 Run 4

Figure 8 for Run 4 is the same as Fig. 6 but we repeat another run for the standard sand and adsorbed solute. The measured porosity is $0.289$ and the flow rate is $11.25\,{\rm cm}^{3}/{\rm h}$. The bed height is $10.7\,{\rm cm}$ and the section area is $7.54\,{\rm cm}^{2}$. By nonlinear least-squares fitting, we obtain $\sigma_{a}=1.0\times 10^{-7}\,{\rm h}^{-1}$, $\sigma_{s}=2.6773\,{\rm h}^{-1}$, $v_{0}=4.1166\,{\rm cm}/{\rm h}$, $u=1.7999\,{\rm cm}/{\rm h}$, $\beta=0.2091$ for LBE, and $u=1.1384\,{\rm cm}/{\rm h}$, $D=3.9814\,{\rm cm}^{2}/{\rm h}$ for ADE. The discrepancy of ADE is even more apparent. The semi-log plot in Fig. 9 shows a discrepancy between the experimental result and ADE.

6 Conclusions

Not only for field observation reported in studies such as Adams and Gelhar (Adams and Gelhar 1992) and Benson and his collaborators (Benson et al. 2001, 2000), anomalous diffusion also appears in lab-scale experiments (Cortis 2004). In Sec. 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, 1993b). 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/C_{0}$ in the tail.

The diffusion approximation holds when the propagation distance of tracer particles is sufficiently larger than the transport mean free path $\ell^{*}$, which is given in (14). In the case of light transport, the photon diffusion equation starts to work when the propagation distance becomes $10$ 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^{\prime}$, from parameters 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 Sec. 5.1.1, $\ell^{*}=1.03$ and $D^{\prime}=1.818$, whereas $D=1.8379$. The relative difference is $|D-D^{\prime}|/D=0.0108$. For Run 2 in Sec. 5.1.2, $\ell^{*}=1.75$ and $D^{\prime}=4.762$, whereas $D=3.8519$. The relative difference is $|D-D^{\prime}|/D=0.236$. In either case, the propagation distance of $18\,{\rm cm}$ or $18.5\,{\rm cm}$ is much larger than $\ell^{*}$ and the transport is expected to be in the diffusion regime. Indeed, Fig. 2 and Fig. 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 Sec. 5.2.1, $\ell^{*}=1.80$ and $D^{\prime}=3.041$, whereas $D=4.3864$. The relative difference is $|D-D^{\prime}|/D=0.307$. For Run 4 in Sec. 5.2.2, $\ell^{*}=1.54$ and $D^{\prime}=2.110$, whereas $D=3.9814$. The relative difference is $|D-D^{\prime}|/D=0.470$. The propagation distance of $10.7\,{\rm cm}$ is not as long as the distance in Run 1 and Run 2 in terms of the transport mean free path. The discrepancy of the curve for ADE in Fig. 6 and Fig. 8 might also attribute to 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 the fact that the surface of the sand is rough and the particle speed cannot be modeled by a constant $v_{0}$.

Although we assumed isotropic scattering in (1), it is possible to take anisotropic scattering into account and write the integral term in a more general form as $\sigma_{s}\int_{-1}^{1}p(\mu,\mu^{\prime})\psi(x,\mu^{\prime},t)\,\mathrm{d}\mu^{\prime}$, where $p(\mu,\mu^{\prime})$ is the probability that tracer particles moving in the direction $\mu^{\prime}$ change their direction to the direction $\mu$ by scattering. We can similarly develop the analytical discrete ordinates method for anisotropic scattering by writing $p(\mu,\mu^{\prime})=\frac{1}{2}\sum_{l}\beta_{l}P_{l}(\mu)P_{l}(\mu^{\prime})$ with some constant coefficients $\beta_{l}$ and Legendre polynomials $P_{l}$ ($l=0,1,\dots$) (Barichello 2011; Siewert 2000).

In this paper we have shown that the mass transport in column experiments, which show 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 $\ell^{*}$ 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 the Interdisciplinary 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) 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 A Eigenvalues

Let us write the homogeneous equation as

[TABLE]

where $I$ is the identity, $\Xi_{\pm}$ and $W$ are matrices, and $\Phi_{\pm}$ are vectors defined as

[TABLE]

We obtain

[TABLE]

where

[TABLE]

We label the eigenvalues as $\nu_{n}$ ($n=1,2,\dots,2N$).

By deforming the Bromwich contour, we obtain a contour $(-\infty,\gamma)$, $\gamma>0$, on which $s$ is real. In this case we have

[TABLE]

where $\lambda_{0}$ is real and $w_{0}$ a $2N$-dimensional real vector. We note that $M=\Re{M}$ in this case and moreover we can write

[TABLE]

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,

[TABLE]

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_{\eta}$ positive eigenvalues. Therefore there are $N_{\eta}$ eigenvalues $\lambda_{0}$.

Next we assume that $s$ has a small imaginary part. Then we can write

[TABLE]

We treat the imaginary part as perturbation and express the matrix-vector equation as

[TABLE]

By collecting terms of order $O(\epsilon^{0})$, we have $S(\Re{M})w_{0}=\lambda_{0}w_{0}$. From terms of $O(\epsilon^{1})$, we have

[TABLE]

Let us multiply $w_{0}^{T}$ on both sides of the above equation. We obtain

[TABLE]

Using $w_{0}^{T}(S(\Re{M}))^{T}=w_{0}^{T}S\Re{M}=\lambda_{0}w_{0}^{T}$, we obtain

[TABLE]

That is, $\lambda_{1}$ is pure imaginary and the number of positive $\Re{\lambda}$ does not change. This fact implies that there are $N_{\eta}$ eigenvalues $\nu_{n}$ ($n=1,\dots,2N$) such that $\Re{\nu_{n}}>0$.

Appendix B Advection-diffusion equation with absorption

Let us consider the following advection-diffusion equation with the absorption term.

[TABLE]

Let us introduce $\Gamma(x,t)$ as

[TABLE]

We have

[TABLE]

Let us consider the Laplace transform:

[TABLE]

Then we obtain

[TABLE]

We obtain

[TABLE]

Thus,

[TABLE]

where

[TABLE]

Noting that

[TABLE]

we obtain

[TABLE]

where the complementary error function is given by $\mathop{\mathrm{erfc}}(x)=(2/\sqrt{\pi})\int_{x}^{\infty}\mathrm{e}^{-t^{2}}\,\mathrm{d}t$. Finally,

[TABLE]

The above solution reduces to the Ogata-Banks solution (Ogata and Banks 1961) when $\sigma_{a}=0$. In particular, we have $C(x,t)\to C_{0}$ as $t\to\infty$ if $\sigma_{a}=0$. When $\sigma_{a}=0$, we have

[TABLE]