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

Anomalous dispersion of solute in porous media can be explained by the power-law distribution of waiting time of solute particles. In this paper, we simulate the diffusion of nonreactive tracer in dead-end pores to explore the waiting time distributions. The distributions of waiting time in different dead-end pores show similar power-law decline at early time and transit to an exponential decline in the end. The transition time between these two decline modes increases with the lengths of dead-end pores. It is well known that power-law distributions of waiting time may lead to anomalous (non-Fickian) dispersion. Therefore, anomalous dispersion is highly dependent on the sizes of immobile zones. According to the power-law decline, we can directly get the power index from the structure of dead-end pores, which can be used to judge the anomalous degree of solute transport in advance.


Introduction
Anomalous dispersion widely exists in solute transport in ground water or in other porous media [1][2][3][4][5][6][7], which makes it difficult to simulate solute transport by traditional advectiondispersion equation.In recent years, scholars have developed a few models to describe the anomalous transport such as continuous time random walk (CTRW) [8][9][10][11][12] and fractional advection-dispersion equation [13][14][15][16].The simulations by CTRW can agree well with the experimental data by fitting the transfer probability density function (, ) [2,17].Assume (, ) has the form (, ) = ()(), where () denotes the waiting time distribution between two successive jumps.Corresponding to anomalous transport, () should follow a power-law decline, that is, () ∼  −1− with 0 <  < 2 [17].Many researchers have investigated the form of () by simulating the transport of solute in porous media [18][19][20][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 () is to directly simulate the motion of solute particles in immobile zones.In this paper, we construct the deadend pores as immobile zones and simulate the waiting time distributions to explore the modes of solute transport.() 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 () [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 .To explore the general form of waiting time distribution function, we directly simulate () in dead-end pores.In this paper, we select three types of dead-end pores.By these deadend pores with different sizes and shapes, we aim to get a universal form of (), 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  which describes the anomalous degree.

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 () of jumps, we only investigate the waiting time distribution ().In our physical picture, waiting events happen when nonreactive particles are trapped in immobile zones.Deadend 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 deadend 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 () in the immobile zone to describe the waiting time distribution (); the relation between them can be written as or Ψ() = −()/.Compared to (), Ψ() fluctuates more violently as the derivative enlarges the fluctuations of simulated values of ().() 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, () is obtained from the derivative of ().However, the probability () is very small when the time is large, especially if () decreases in an exponential form.By conventional particle tracing method with each particle of the same weight, () is proportional to the number of particles staying in the immobile zone.The number corresponding to small () is zero or fluctuates violently.To simulate () 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][23][24].We use the pruned-enriched method to adjust the weight according to ().At the beginning of simulation, each particle is released at the interface between the flow tube and dead-end pores with the same weight (0) = 1.Consider the case that the simulated particle is still in the dead-end pore at the time  = Δ.Define  = ()/().If  > 1, this particle is split in Int() (the integer part of ) copies and each copy continues to diffuse with a new weight  new = ()/Int().If  ≤ 1, we continue to simulate the diffusion of the particle with probability   = 1/Int(1/) and the new weight  new is adjusted to ()/  .In other words, the simulation of the particle is ceased with the probability 1 −   .() 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  has a nearly uniform distribution.Without the scheme of variable weights, the number of visited samples is proportional to () and decreases sharply with , which leads to violent oscillation of simulated ().

Results and Discussion
For most heavy metal ions, the order of magnitude of diffusion coefficient is about 10 −10 m 2 /s.In this paper, the diffusion coefficient of solute is set to be 4 × 10 −10 m 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 deadend pore again because they drift in the flow tube (mobile zone).Firstly, we consider the surviving probability () 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 () and ().The () and () 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 () because () has more violent fluctuations.
During the time interval 10 s-10 4 s, () decays in a power law, that is, () ∼  − .Corresponding to straight tube, bent tube, and hemisphere, the value of  is 0.50, 0.49, and 0.52, respectively.They are close to 0.5.Thus, it is reasonable to suppose that the power  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 () with  = 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 () ∼  −1− , the concentration (, ) can be written as [25]  (, ) where 0  1−  is the fractional Riemann-Liouville operator; that is, Here, [(, )] 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 () 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 () and () are illustrated in Figure 3.At the early time, () and () in each straight tube decrease with the power law.As time increases, () and () 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 () 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 () with a high resolution.Figure 4 shows that log(()) and log(()) are linear with time  in every type of dead-end pore, which indicates the exponential decline when  is large.The exponential decline can be easily explained by immobile-mobile model [26] by neglecting the variability of concentration in immobile zones, which reads where  im and  m represent the concentration in dead-end pores (immobile zones) and in flow tube (mobile zones), respectively. im (, ) is proportional to ().If  m = 0 at large  because of drift with water flow,  im decreases as  im ∼  − , that is, an exponential decline.
From simulations above, we can see that surviving probability () decreases with time in a power law at early time and in an exponential form at large time.The power law with  ≈ 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   at which the decline transits from a power law to an exponential form.To get the transition time   , we first calculate the exponent  in (4).  has a close relation with the reciprocal of .For simplicity, we use a straight tube as the dead-end pore.The diffusion equation about the concentration (, ) in the tube can be written as where  is the length of the straight tube and () = ∫  0 (, ) .The interface between the straight tube and the mobile zone is located at  = 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  = 0 can be set as (0, ) = 0.The condition of reflecting boundary at  =  is written as (, )/| = = 0.By the method of variable separation, (, ) has the form (, ) = () − and we can get the minimum of ; that is, Other eigenvalues of  are much larger than  min and the modes with large  decrease more rapidly than the mode with  min .Neglecting the modes with large ,  in (, ) = () − approximates to  min .We fit the  of straight tubes of different length according to () ∼  − and compare them with the values of  min calculated from (6).Table 1 shows the values of  and  min .The two sets of values are very close.
In hemisphere dead-end pores,  can also be approximated by  min as the tube despite replacing  with  and adding a coefficient.Now let us discuss the transition time   at which the decline of () transits from a power-law form to an exponential decline form.In Figure 3,   is the time when the curve deviates from the straight line.Comparing the transition time   to 1/ in Table 1, we can find that   is located in the range of time [1/2, 1/].Since  is very close to  min (illustrated in Table 1),   is approximately proportional to the square of length  and inversely proportional to  m as shown in (6).As   represents the transition time from power-law decline to exponential decline, the enduring time of anomalous dispersion has the same relation with the length  and  m as that in (6).

Conclusions
In this paper, we investigate the waiting time distribution () of nonreactive solute particles in immobile zones, which are represented by dead-end pores.Three types of deadend pores, straight tube, bent tube, and hemisphere, are simulated.When a dead-end pore is large enough, () follows a power-law distribution, that is, () ∼  −1− .All the simulated dead-end pores in this paper show that the values of  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, () follows an exponential decrease finally.Therefore, the transport is anomalous at early time and is Fickian at later time.The transition time   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   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 powerlaw distribution of () has a close relation with fractional dispersion equation [25].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.

Figure 1 :
Figure1: 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.

Figure 2 :
Figure2: The power-law form of the surviving probability () and the waiting time distribution Ψ() in dead-end pores at early time.The power  approximates to 0.5 for each type of dead-end pore.

Figure 4 :
Figure 4: The exponential decrease of () and Ψ() at large time for each type of dead-end pore.

Table 1 :
The comparison between the  fitted by () ∼  − and the  min of (6).