# The Form of Waiting Time Distributions of Continuous Time Random Walk in Dead-End Pores.

1. Introduction

Anomalous dispersion widely exists in solute transport in ground water or in other porous media [1-7], which makes it difficult to simulate solute transport by traditional advection-dispersion equation. In recent years, scholars have developed a few models to describe the anomalous transport such as continuous time random walk (CTRW) [8-12] and fractional advection-dispersion equation [13-16]. The simulations by CTRW can agree well with the experimental data by fitting the transfer probability density function [phi](x, t) [2, 17]. Assume [phi](x, t) has the form [phi](x, t) = [phi](x)[psi](t), where [psi](t) denotes the waiting time distribution between two successive jumps. Corresponding to anomalous transport, [psi](t) should follow a power-law decline, that is, [psi](t) ~ [t.sup.-1-[alpha]] with 0 < [alpha] < 2 . Many researchers have investigated the form of f(t) by simulating the transport of solute in porous media [18-21], which may show power-law decline. However, fitting parameters from transport data in complex fluid fields are not suitable to uncover the physical mechanism of the power-law decline. An alternative way for the investigation of the form of [psi](t) is to directly simulate the motion of solute particles in immobile zones. In this paper, we construct the dead-end pores as immobile zones and simulate the waiting time distributions to explore the modes of solute transport.

[psi](t) is more likely to show an exponential decline than a power-law decline because particles in stagnant zones are determined by diffusion equation. As many experiments show, an anomalous transport eventually converts to a Fickian process [2, 19]. Consequently, in the mathematical models, the power-law decline is often truncated arbitrarily by an exponential decline in the assumption of [psi](t) [4, 9]. However, why does there exist a power-law decline of waiting time distribution? When does a power-law decline begin to convert to an exponential decline? In addition, it is not clear which factors determine the value of the power index [alpha]. To explore the general form of waiting time distribution function, we directly simulate [psi](t) in dead-end pores. In this paper, we select three types of dead-end pores. By these dead-end pores with different sizes and shapes, we aim to get a universal form of [psi](t), which may show both the power-law part (corresponding to anomalous transport) and the exponential decline part (Fickian transport). In addition, from our simulations, we expect to get the index a which describes the anomalous degree.

2. Method of Simulation

In the framework of CTRW model for solute transport, the distribution of waiting time between two successive jumps determines whether the dispersion is anomalous or not. For simplicity, skipping over the displacement distribution [phi](x) of jumps, we only investigate the waiting time distribution [psi](t). In our physical picture, waiting events happen when nonreactive particles are trapped in immobile zones. Dead-end pores are typical immobile zones. The flow velocity in a dead-end pore is zero and solute particles escape only by diffusion. The waiting time is the period from entering a dead-end pore to moving out of it. We consider three types of dead-end pores: straight tube, bent tube, and hemisphere as illustrated in Figure 1. These three dead-end pores lie above the flow tube which represents the mobile zone. Tracer particles on the interface between the flow tube and dead-end pores may enter the dead-end pore by diffusion. We begin to count waiting time when the particles pass through interfaces to dead-end pores and end counting when particles pass through interfaces to the flow tube.

For simplicity, we use the surviving probability P(t) in the immobile zone to describe the waiting time distribution [psi](t); the relation between them can be written as

P (t) = 1 - [[integral].sup.t.sub.0] [psi]([tau]) d[tau] = [[integral].sup.[infinity].sub.t] [psi](t) d[tau], (1)

or [PSI](t) = -dP(t)/dt. Compared to P(t), [PSI](t) fluctuates more violently as the derivative enlarges the fluctuations of simulated values of P(t).

P(t) can be directly calculated by counting the number of particles in the dead-end pores and the particles will not be tracked when they move out of immobile zones. Then, [psi](t) is obtained from the derivative of P(t). However, the probability P(t) is very small when the time is large, especially if P(t) decreases in an exponential form. By conventional particle tracing method with each particle of the same weight, P(t) is proportional to the number of particles staying in the immobile zone. The number corresponding to small P(t) is zero or fluctuates violently. To simulate P(t) in a wide range, variable weight for tracing particles is adopted in this paper. In recent years the pruned-enriched method has been developed to simulate the random walk with variable weight [22-24]. We use the pruned-enriched method to adjust the weight according to P(t). At the beginning of simulation, each particle is released at the interface between the flow tube and dead-end pores with the same weight W(0) = 1. Consider the case that the simulated particle is still in the dead-end pore at the time t = n[DELTA]t. Define r = W(t)/P(t). If r > 1, this particle is split in Int(r) (the integer part of r) copies and each copy continues to diffuse with a new weight [W.sub.new] = W(t)/Int(r). If r [less than or equal to] 1, we continue to simulate the diffusion of the particle with probability r' = 1/Int(1/r) and the new weight [W.sub.new] is adjusted to W(t)/r'. In other words, the simulation of the particle is ceased with the probability 1 - r'. P(t) is updated in every step. By the scheme of variable weight, the smallest probability we could simulate can be reduced by many orders of magnitude. In the whole range of simulated time, the number of samples for each t has a nearly uniform distribution. Without the scheme of variable weights, the number of visited samples is proportional to P(t) and decreases sharply with t, which leads to violent oscillation of simulated P(t).

3. Results and Discussion

For most heavy metal ions, the order of magnitude of diffusion coefficient is about [10.sup.-10] [m.sup.2]/s. In this paper, the diffusion coefficient of solute is set to be 4 x [10.sup.-10] [m.sup.2]/s. When tracer particles enter the dead-end pores through the interface, we start the time. When particles escape from the dead-end pores, the tracing process of this particle will be terminated. For simplicity, particles after escaping from a dead-end pore cannot be allowed to enter the same dead-end pore again because they drift in the flow tube (mobile zone). Firstly, we consider the surviving probability P(t) in large dead-end pores. The lengths of straight tube and bent tube are both taken 10 mm and the radius of hemisphere is taken 15 mm. Figure 2 shows the power-law decrease of simulated P(t) and f(t). The P(t) and [psi](t) in bent tube and in hemisphere are reduced by the factor 1/2 and 1/4, respectively, which makes the comparison in Figure 2 clear without altering the power law. We take the values of power from the data of P(t) because [psi](t) has more violent fluctuations. During the time interval 10 s-[10.sup.4] s, P(t) decays in a power law, that is, P(t) ~ [t.sup.~[alpha]]. Corresponding to straight tube, bent tube, and hemisphere, the value of [alpha] is 0.50, 0.49, and 0.52, respectively. They are close to 0.5. Thus, it is reasonable to suppose that the power [alpha] is approximately equal to 0.5 in these three types of dead-end pores.

Let us discuss the physical meaning of the simulated results. This power-law form of P(t) with [alpha] = 0.5 results in anomalous dispersion, which can also be described by time-fractional advection-dispersion equation (tFADE). If the waiting time distribution has the form [psi](t) ~ [t.sup.-1-[alpha]], the concentration C(x, t) can be written as 

[partial derivative]C (x, t)/[partial derivative]t = [sub.0][D.sup.1-[alpha].sub.t] L [C (x, t)], (2)

where [sub.0][D.sup.1-[alpha].sub.t] is the fractional Riemann-Liouville operator; that is,

[sub.0][D.sup.1-[alpha].sub.t] f (t) = 1/[GAMMA] ([alpha]) [partial derivative]/[partial derivative]t [[integral].sup.t.sub.0] d[tau] f ([tau])/[(t - [tau]).sup.1-[alpha]]. (3)

Here, L[C(x, t)] is the operator representing dispersion. Therefore, the fractional order of time derivation is about 0.5 if anomalous dispersion is caused by solute trapping in long dead-end pores.

The sizes of above dead-end pores are set large in the simulations. Let us consider the dead-end pores with short length. Tracer particles can jump out of them easily. Which form does surviving probability density P(t) show? To investigate it, we set the length of straight tubes to be 1 mm, 2 mm, 3 mm, and 4 mm, respectively. The diffusion coefficient is the same as above. The simulated P(t) and [psi](t) are illustrated in Figure 3. At the early time, P(t) and [psi](t) in each straight tube decrease with the power law. As time increases, P(t) and [psi](t) deviate from the power-law decline. The shorter the length is, the earlier the deviation happens. Therefore, anomalous dispersion strongly depends on the sizes of dead-end pores. The duration of anomalous dispersion becomes longer as the length of dead-end pore increases.

Now let us consider the form of P(t) at large time in small dead-end pores. We also investigate three types of dead-end pores, that is, straight tube, bent tube, and hemisphere, but their sizes are shortened. The length of straight tube and bent tube is set to be 2 mm and 3 mm, respectively. The radius of hemisphere is set to be 3 mm. Using the particle tracing method with viable weight, we can simulate the surviving probability P(t) with a high resolution. Figure 4 shows that log(P(t)) and log([psi](t)) are linear with time t in every type of dead-end pore, which indicates the exponential decline when t is large. The exponential decline can be easily explained by immobile-mobile model  by neglecting the variability of concentration in immobile zones, which reads

[partial derivative][C.sub.im] (x, t)/[partial derivative]t = -[beta] ([C.sub.im] - [C.sub.m]), (4)

where [C.sub.im] and [C.sub.m] represent the concentration in dead-end pores (immobile zones) and in flow tube (mobile zones), respectively. [C.sub.im] (x, t) is proportional to P(t). If [C.sub.m] = 0 at large t because of drift with water flow, [C.sub.im] decreases as [C.sub.im] ~ [e.sup.-[beta]t] that is, an exponential decline.

From simulations above, we can see that surviving probability P(t) decreases with time in a power law at early time and in an exponential form at large time. The power law with [alpha] [approximately equal to] 0.5 corresponds to an anomalous dispersion, while the exponential decline corresponds to a Fickian dispersion. Therefore, it is important to determine the transition time [T.sub.t] at which the decline transits from a power law to an exponential form. To get the transition time [T.sub.t], we first calculate the exponent [beta] in (4). [T.sub.t] has a close relation with the reciprocal of [beta]. For simplicity, we use a straight tube as the dead-end pore. The diffusion equation about the concentration C(x, t) in the tube can be written as

[partial derivative]C (x, t)/[partial derivative]t = [D.sub.m] [[partial derivative].sup.2]C (x, t)/[partial derivative][x.sup.2] (0 [less than or equal to] x [less than or equal to] L), (5)

where L is the length of the straight tube and P(t) = [[integral].sup.L.sub.0] C(x, t) dx. The interface between the straight tube and the mobile zone is located at x = 0. As tracer particles move out of the straight tube, they are not allowed to enter again because of drift with flow. Consequently, the boundary condition at x = 0 can be set as C(0, t) = 0. The condition of reflecting boundary at x = L is written as [partial derivative]C(x, t)/[partial derivative]x[|.sub.x=L] = 0. By the method of variable separation, C(x, t) has the form C(x, t) = X(x)[e-.sup.[beta]t] and we can get the minimum of [beta]; that is,

[[beta].sub.min] = [[pi].sup.2][D.sub.m]/4[L.sup.2]. (6)

Other eigenvalues of [beta] are much larger than [[beta].sub.min] and the modes with large [beta] decrease more rapidly than the mode with [[beta].sub.min]. Neglecting the modes with large [beta], [beta] in C(x, t) = X(x)[e.sup.-[beta]t] approximates to [[beta].sub.min]. We fit the [beta] of straight tubes of different length according to P(t) ~ [e.sup.[beta]t] and compare them with the values of [[beta].sub.min] calculated from (6). Table 1 shows the values of [beta] and [[beta].sub.min]. The two sets of values are very close. In hemisphere dead-end pores, [beta] can also be approximated by [[beta].sub.min] as the tube despite replacing L with R and adding a coefficient.

Now let us discuss the transition time [T.sub.t] at which the decline of P(t) transits from a power-law form to an exponential decline form. In Figure 3, [T.sub.t] is the time when the curve deviates from the straight line. Comparing the transition time [T.sub.t] to l/[beta] in Table 1, we can find that [T.sub.t] is located in the range of time [1/2[beta], 1/[beta]]. Since [beta] is very close to [[beta].sub.min] (illustrated in Table 1), [T.sub.t] is approximately proportional to the square of length L and inversely proportional to [D.sub.m] as shown in (6). As [T.sub.t] represents the transition time from power-law decline to exponential decline, the enduring time of anomalous dispersion has the same relation with the length L and [D.sub.m] as that in (6).

4. Conclusions

In this paper, we investigate the waiting time distribution [psi](t) of nonreactive solute particles in immobile zones, which are represented by dead-end pores. Three types of dead-end pores, straight tube, bent tube, and hemisphere, are simulated. When a dead-end pore is large enough, [psi](t) follows a power-law distribution, that is, [psi](t) ~ [t.sup.-1-[alpha]]. All the simulated dead-end pores in this paper show that the values of a are very close to 0.5. According to CTRW theory, this power-law distribution results in anomalous dispersion [17, 25]. If dead-end pore is finite, [psi](t) follows an exponential decrease finally. Therefore, the transport is anomalous at early time and is Fickian at later time. The transition time [T.sub.t] is proportional to the square of the length of dead-end pores and inversely proportional to the diffusion coefficient of solute. From CTRW theory, we could use [T.sub.t] to sign the transition time from anomalous transport to Fickian transport. Thus, we can draw the conclusion that anomalous transport phenomena will be more easily observed in the porous media with large immobile zones.

In practice, dead-end pores may be very long, such as the fractures in fracture porous media. In this case, the waiting time distribution can seem as a power-law distribution during the whole observation period. Therefore, anomalous dispersion lasts for the whole observation period. The power-law distribution of [psi](t) has a close relation with fractional dispersion equation . According to our simulations, the fractional order of time derivative approximates to 0.5 when using tFADE (i.e., (2)) to describe this class of anomalous dispersion.

https://doi.org/10.1155/2018/8329406

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported by the National Key Research and Development Program of China (Grant no. 2016YFC040280801), the National Natural Science Foundation of China- (NSFC-) Xinjiang project (Grant no. U1503282), and NSFC projects (Grant no. 41672233).

References

 E. E. Adams and L. W. Gelhar, "Field study of dispersion in a heterogeneous aquifer: 2. spatial moments analysis," Water Resources Research, vol. 28, no. 12, pp. 3293-3307, 1992.

 M. Levy and B. Berkowitz, "Measurement and analysis of non-fickian dispersion in heterogeneous porous media," Journal of Contaminant Hydrology, vol. 64, no. 3-4, pp. 203-226, 2003.

 U. M. Scheven, D. Verganelakis, R. Harris, M. L. Johns, and L. F. Gladden, "Quantitative nuclear magnetic resonance measurements of preasymptotic dispersion in flow through porous media," Physics of Fluids, vol. 17, no. 11, Article ID 117107, pp. 1-7, 2005.

 Y. Zhang, H. Xu, X. Lv, and J. Wu, "Diffusion in relatively-homogeneous sand columns: a scale-dependent or scale-independent process?" Entropy, vol. 15, pp. 4376-4391, 2013.

 E. Crevacore, T. Tosco, R. Sethi, G. Boccardo, and D. L. Marchisio, "Recirculation zones induce non-fickian transport in three-dimensional periodic porous media," Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, vol. 94, no. 5, Article ID 053118, 2016.

 A. Tyukhova, M. Dentz, W. Kinzelbach, and M. Willmann, "Mechanisms of anomalous dispersion in flow through heterogeneous porous media," Physical Review Fluids, vol. 1, Article ID 074002, 2016.

 A. Nissan, I. Dror, and B. Berkowitz, "Time-dependent velocity-field controls on anomalous chemical transport in porous media," Water Resources Research, vol. 53, pp. 3760-3769, 2017

 B. Berkowitz, G. Kosakowski, G. Margolin, and H. Scher, "Application of continuous time random walk theory to tracer test measurements in fractured and heterogeneous porous media," Groundwater, vol. 39, no. 4, pp. 593-604, 2001.

 M. Dentz, A. Cortis, H. Scher, and B. Berkowitz, "Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport," Advances in Water Resources, vol. 27, no. 2, pp. 155-173, 2004.

 R. Poloczanski, A. Wylomanska, M. Maciejewska, A. Szczurek, and J. Gajda, "Modified cumulative distribution function in application to waiting time analysis in the continuous time random walk scenario," Journal of Physics A: Mathematical and Theoretical, vol. 50, no. 3, Article ID 034002, 2017

 A. Puyguiraud, M. Dentz, and P. Gouze, "Transport upscaling from pore-to Darcy-scale: incorporating pore-scale berea sandstone lagrangian velocity statistics into a Darcy-scale transport CTRW model," in EGU General Assembly Conference Abstracts, vol. 19, p. 10943, 2017

 M. F. Shlesinger, "Origins and applications of the montroll-weiss continuous time random walk," The European Physical Journal B, vol. 90, no. 5, 93 pages, 2017

 R. Metzler and J. Klafter, "The random walk's guide to anomalous diffusion: a fractional dynamics approach," Physics Reports, vol. 339, no. 1, pp. 1-77, 2000.

 Y. Zhang, D. A. Benson, and D. M. Reeves, "Time and space nonlocalities underlying fractional-derivative models: distinction and literature review of field applications," Advances in Water Resources, vol. 32, no. 4, pp. 561-581, 2009.

 J. F. Kelly and M. M. Meerschaert, "Space-time duality for the fractional advection-dispersion equation," Water Resources Research, vol. 53, pp. 3464-3475, 2017

 R. Schumer, D. A. Benson, M. M. Meerschaert, and S. W. Wheatcraft, "Eulerian derivation of the fractional advection-dispersion equation," Journal of Contaminant Hydrology, vol. 48, no. 1-2, pp. 69-88, 2001.

 B. Berkowitz, A. Cortis, M. Dentz, and H. Scher, "Modeling Non-fickian transport in geological formations as a continuous time random walk," Reviews of Geophysics, vol. 44, no. 2, Article ID RG2003, 2006.

 B. Bijeljic and M. J. Blunt, "Pore-scale modeling and continuous time random walk analysis of dispersion in porous media," Water Resources Research, vol. 42, no. 1, Article ID W01202, 2006.

 B. Berkowitz and H. Scher, "Exploring the nature of non-fickian transport in laboratory experiments," Advances in Water Resources, vol. 32, no. 5, pp. 750-755, 2009.

 B. Bijeljic, A. Raeini, P. Mostaghimi, and M. J. Blunt, "Predictions of non-fickian solute transport in different classes of porous media using direct simulation on pore-scale images," Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, vol. 87, no. 1, Article ID 013011, 2013.

 J. Jiang and J. Wu, "Continuous time random walk in homogeneous porous media," Journal of Contaminant Hydrology, vol. 155, pp. 82-86, 2013.

 T. Prellberg and J. Krawczyk, "Flat histogram version of the pruned and enriched rosenbluth method," Physical Review Letters, vol. 92, no. 12, pp. 120602-1, 2004.

 J. Jiang, Y. Huang, and J. Wu, "Can the pruned-enriched method be used for the simulation of fluids?" Journal of Statistical Physics, vol. 136, no. 5, pp. 984-988, 2009.

 J. Jiang and J. Wu, "Simulation of continuous-time random walks by the pruned-enriched method," Physical Review E, vol. 84, Article ID 036710, 2011.

 E. Barkai, R. Metzler, and J. Klafter, "From continuous time random walks to the fractional Fokker-Planck equation," Physical Review E: Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics, vol. 61, no. 1, pp. 132-138, 2000.

 P. P. Wang, C. Zheng, and S. M. Gorelick, "A general approach to advective-dispersive transport with multirate mass transfer," Advances in Water Resources, vol. 28, no. 1, pp. 33-42, 2005.

Yusong Hou, Jianguo Jiang (iD), and J. Wu

School of Earth Sciences and Engineering, Nanjing University, Nanjing, China

Correspondence should be addressed to Jianguo Jiang; jianguo.jiang@nju.edu.cn

Received 22 July 2017; Revised 16 November 2017; Accepted 12 December 2017; Published 8 January 2018

Caption: Figure 1: The schematic diagram of mobile zone and immobile zone. Three types of dead-end pores, that is, straight tube, bent tube, and hemisphere, represent the immobile zones.

Caption: Figure 2: The power-law form of the surviving probability P(t) and the waiting time distribution W(t) in dead-end pores at early time. The power [alpha] approximates to 0.5 for each type of dead-end pore.

Caption: Figure 3: The deviation from the power-law decline of P(t) and [PSI](t) in straight tubes of different length L as t increases. The shorter the length is, the faster the deviation happens.

Caption: Figure 4: The exponential decrease of P(t) and [PSI](t) at large time for each type of dead-end pore.
```Table 1: The comparison between the [beta] fitted by
P(t) ~ [e.sup.-[beta]t] and the [[beta].sub.min] of (6).

[[beta].sub.min]
L (mm)      [beta] ([s.sup.-1])          ([s.sup.-1])

1            9.72 x [10.sup.-4]       9.87 x [10.sup.-4]
2            2.44 x [10.sup.-4]       2.47 x [10.sup.-4]
3            1.08 x [10.sup.-4]       1.10 x [10.sup.-4]
4            6.10 x [10.sup.-4]       6.20 x [10.sup.-4]

[beta]/
L (mm)      [[beta].sub.min]

1                  0.98
2                  0.99
3                  0.98
4                  0.99
```