Abstract
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.
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 [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, with [17]. Many researchers have investigated the form of 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 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.
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 dead-end 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.
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 of jumps, we only investigate the waiting time distribution . 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 in the immobile zone to describe the waiting time distribution ; the relation between them can be written asor . 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–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 . Consider the case that the simulated particle is still in the dead-end pore at the time . Define . If , this particle is split in (the integer part of ) copies and each copy continues to diffuse with a new weight . If , we continue to simulate the diffusion of the particle with probability and the new weight is adjusted to . In other words, the simulation of the particle is ceased with the probability . 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 .
3. Results and Discussion
For most heavy metal ions, the order of magnitude of diffusion coefficient is about 10−10 m2/s. In this paper, the diffusion coefficient of solute is set to be m2/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 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–104 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 results in anomalous dispersion, which can also be described by time-fractional advection-dispersion equation (tFADE). If the waiting time distribution has the form , the concentration can be written as [25]where 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 and 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 readswhere and represent the concentration in dead-end pores (immobile zones) and in flow tube (mobile zones), respectively. is proportional to . If at large because of drift with water flow, decreases as , 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 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 aswhere is the length of the straight tube and . The interface between the straight tube and the mobile zone is located at . 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 can be set as . The condition of reflecting boundary at is written as . 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 and the modes with large decrease more rapidly than the mode with . Neglecting the modes with large , in approximates to . We fit the of straight tubes of different length according to and compare them with the values of calculated from (6). Table 1 shows the values of and . The two sets of values are very close. In hemisphere dead-end pores, can also be approximated by 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 in Table 1, we can find that is located in the range of time . Since is very close to (illustrated in Table 1), is approximately proportional to the square of length and inversely proportional to 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 as that in (6).
4. 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 dead-end pores, straight tube, bent tube, and hemisphere, are simulated. When a dead-end pore is large enough, follows a power-law distribution, that is, . 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 power-law 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.
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).