The understanding of solute transport in rock fractures is of great importance in many engineering fields. In this study, two groups of experiments on artificial and natural single fractures with different fracture apertures and roughness were conducted to investigate the changes of solute transport regimes. The time fractional advection–dispersion equation (tFADE) as a promising model was applied to describe the anomalous transport. The performance of the classical advection-dispersion equation (ADE) and tFADE was compared according to the fitting precision of breakthrough curves (BTCs). The responses of the fitting parameters in the tFADE to the experimental conditions were also discussed. The results indicated that the non-Fickian transport more likely occurred in the short distance transport, and the larger Peclet number () led to the increase of the exponent of the power-law function in the phase of concentration decline. The tFADE was superior to ADE in capturing the non-Fickian transport especially the tailing behavior. The fractional order of time in the tFADE was the key parameter to describe the anomalous transport process, and its responding mechanisms of were revealed: the best-fit decreased with the increase of flow velocity and the decrease of the fracture aperture. The roughness of the single fracture which leads to a complex flow field had a significant effect on the best-fit . The findings of this study can help for better understanding the effectiveness and physical significance of the tFADE.

1. Introduction

Solute transport in fractured rocks is of great significance to the underground engineering such as tunnel excavation, shale gas development, and CO2 geological storage. Due to the variable geometry and hydraulic conditions in the geological environment, the accurate quantitative prediction of solute transport in the fracture is an important research with great challenges [1, 2].

In order to predict the solute transport process in fractured media, mathematical models have been developed and applied in the laboratory investigation and numerical simulation [35]. The acknowledged advection–dispersion equation (ADE), derived from Fick’s law, has been used to describe solute transport in fractures for years [6]. However, many laboratory researches claimed that medium heterogeneity deeply affects the transport of solute in the fracture, resulting in the limitation of the ADE model in describing the early arrival or tailing behavior which is typical non-Fickian behavior [79]. In single fractures, the geometry structure of walls and inertial effects of flow are mainly considered as the causes of the non-Fickian phenomena [10, 11]. That is, the rougher fracture and higher velocity of the flow lead to the differentiation of the main flow paths and recirculation zones which affect the transport process deeply and cause the failures of the ADE (Figure 1) [12, 13]. Moreover, non-Fickian phenomena are common in nature [14]. Thus, quantitative capture of non-Fickian phenomenon at a laboratory scale contributes to the study of solute transport in a natural scale [15].

To improve the precision accuracy of pollutant transport and capture the non-Fickian transport, several mathematical approaches have been developed [16, 17]. Based on a transition rate approach, continuous time random walk (CTRW) which includes the Fokker-Panck equation with a memory equation has been widely used to capture solute transport in both the fracture system and heterogeneous porous media [18, 19]. Nowamooz et al. [20] interpreted the tracer transport in the original fractures with the equivalent-stratified model and the continuous time random walk (CTRW) and analyzed the advantages and disadvantages of the two models. Though CTRW can capture the transport process well, the physical relationship between the transport velocity and the average fluid velocity is not clear enough and the parameters of CTRW are relatively more than those of other models which result in the difficulty in simulation. Improved on the ADE, several mathematical approaches can also capture the solute transport process with a certain degree of accuracy like the mobile-immobile (MIM) model and fractional advection-dispersion equation (FADE) [21, 22]. Qian et al. [23] used the single-rate MIM model, ADE, and advection–dispersion equation with retardation (ADE-R) to fit the sodium chloride transport in the filled fracture and test the precision of the MIM. The results showed that MIM can capture the anomalous transport process well including early arrival and tailing behaviors. The equivalent-stratified method can quantify the degree of fracture heterogeneity, but the goodness of fit for the solute transport process of short-distance transportation is low. By adding a fractional-order operator to time in the traditional ADE, the time fractional advection-dispersion equation (tFADE) can describe the solute remaining in an immobile domain [24, 25]. Hence, the tFADE has the potential to describe the tailing behavior caused by the geometry structure of the single fracture. On the other hand, the spatial fractional advection-dispersion equation (sFADE) has been used in describing solute transport for many years which can describe the solute transport at different transport distances [3, 4, 9]. Benson et al. [26] regarded FADE as a predictive tool with certain precision by fitting the transport in a relatively homogeneous sandbox. Huang et al. [27] simulated the atrazine transport in a sand column by numerical approach and found that the FADE can capture well the non-Fickian transport especially the tailing behaviors. Hence, the FADE is a promising approach to describe the solute transport process in the fracture system with a complex interference factor. However, the physical significance of the fractional order is still unclear.

In summary, the capture of non-Fickian transport still remains to be further researched. The effectiveness of the FADE to the simulated solute transport process under different conditions at the laboratory scale needs to be proved. Researches on the physical significance and influence factors of the parameters in the FADE are necessary. Hence, systematical studies of the FADE for describing the solute transport in a single fracture are still a requisite [28]. It can provide strong theoretical support for the wider application of the FADE.

This study is aimed at analyzing the cause of the non-Fickian transport process in single fractures, testing the fitting precision of the tFADE under different experimental conditions, and investigating the effect of the flow rate, fracture aperture, and roughness in granite fractures and plexiglass fractures on the parameters of the tFADE.

The research steps are organized as follows. First, the theoretical background of solute transport is introduced in Section 2. Then, 2 groups of solute transport experiments are performed in Section 3. Third, the experimental results and fitting results with ADE and tFADE are described, as well as the comparison of the fitting of the two models and discussion of the influence of experimental conditions on the fitting accuracy and the parameters of the tFADE. Finally, the conclusion and limitation of this study are summarized in Section 5.

2. Theoretical Background of Solute Transport

The fractional operator has been developed in many fields for years [29]. Mehdinejadiani and Fathi [30] simulated water table profiles using the space fractional Boussinesq equation. Jafari et al. [31] developed a fractional Glover-Dumm equation to simulate the groundwater flow under heterogeneous unconfined aquifers. The fractional-derivative operator of the FADE can be used to describe anomalous diffusion for its nonlocality [32]. Many studies have shown that the FADE can describe the abnormal transport of particles [4, 33, 34]. The time fractional advection-dispersion equation (tFADE) is an upgrade based on the ADE by replacing a fractional derivative from 0 to 1 to the first time derivative, arising from the power law of the time distribution of particles [35, 36]. The basic form of the tFADE is as follow: where is the concentration of the solute (mol/L), is the dispersion coefficient (m2/s), is the advective velocity (m/s), is the time (s), is the spatial coordinate of the solute along the direction of flow (m), and (dimensionless) is the order of fractional differentiation. The smaller value of explains the more remarkable solute retention in the fracture. In addition, when the value is 1, the tFADE recedes to the ADE. The dispersion coefficient in the tFADE is different from that in the ADE, because the fractional operator describes the partially local velocity rather than the mean velocity.

The numerical solutions of equation (2) are obtained by finite difference approximation, and the initial and the boundary conditions of solute transport are shown as follows: where represents the duration of solution injection and represents the injected concentration and, in this study, is set to 1 for the normalized concentration.

To evaluate the fitting accuracy of the tFADE, the coefficient of determination () and root mean square error (RMSE) were calculated. The coefficient of determination can express the goodness of fit by calculating the proportion of variability in the dependent variable that can be explained in the FADE, and the range of is from 0 to 1. RMSE is a non-negative number which reflects the measure of dispersion of an observation [37]. It should be mentioned that the normalized concentration is used in this study, so both parameters are dimensionless. In the experiments, the closer near to 1 and the smaller RMSE represent the better fitting result. The mathematical expressions of and RMSE are as follows: in which is the observed concentration, is the average value of the observed concentration, and is the simulated concentration.

In this study, two important dimensionless numbers are used. The Reynolds number () is a classic quantity to describe the ratio of inertial effect to viscous effects on the field flow, and the Peclet number () is defined as the ratio of the convection, and the diffusion can be used to describe the importance of convective flux. The larger indicates the more significant convective transport [38]. The computational expressions are shown as follows: where is the fluid density (103 kg/m3), is the average flow velocity, is the apparent fracture aperture, is the dynamic fluid viscosity (10−3 Pa·s), and is the molecular diffusion coefficient.

The residence-time distributions (RTDs), effective metric for evaluate the features of non-Fickian transport, are developed by the nonlocal dispersion theory [39, 40]. In this study, RTDs were calculated for the solute transport in the granite fracture and plexiglass fracture to examine the features of tailing behavior. RTDs were derived from the time derivative of BTCs of the step injection for its resident concentration time series, and the residence time distribution can be expressed as follows [41]: where is a dimensionless time called pore volume, is the flow rat e(), and is the volume of the fracture (). It should be noted that a single instantaneous injection was simulated in our study and the BTCs can be converted from the single instantaneous injection to the step injection [42]. where the and represent the concentration obtained in response to the step injection and pulse injection, respectively. According to equation (8), equation (9) can be transformed into

3. Experimental Apparatus and Procedure

3.1. Experimental Apparatus

In this study, two groups of experiments in a single fracture under different experimental conditions were conducted and each group involved two patterns.

The schematic diagram of the laboratory setup are shown in Figure 2. Three main parts including a recharge flume, a discharge flume, and a single fracture are assembled in group 1. Both the recharge flume and the discharge flume consist of a water pipe for the control device and an overflow launder. Three hydraulic gradients for each pattern in group 1 are designed by adjusting the height of the recharge flume with the constant position of the discharge flume. The steady flow rates are calculated by a chronograph and a measuring cylinder after the water head is stable. Sampling analyses are taken at certain intervals to monitor the change of solute concentration. Each experiment was repeated three times.

Brilliant blue FCF (BBF), an ionic molecule, was used in several solute transport experiments as a conservative tracer recently for its nontoxicity, mobility, and low adsorption [4345]. To measure the concentration of BBF in the fracture, the standard curve is calculated by configuring a series of concentrations of BBF solutions and measuring the absorbance. The value of absorbance of the sample is measured at 630 nm by a spectrophotometer SP-754 (Shanghai Precision and Scientific Instrument) in group 1 and 752 N (Shanghai Precision & Scientific Instrument Co. Ltd., China) in group 2.

Two patterns of fractures were contained in group 1. The natural and rough fracture was made up of two granite plates 59 cm long and 12 cm wide. The average aperture of the natural fracture is calculated by measuring the width of 110 points with a vernier caliper because of unobvious depressions, and projections are present on the granite plate, and the apertures of 1.71 and 2.60 mm are set. Another pattern is made up of two parallel plexiglass plates, which were 69 cm long, 13 cm wide, and 1 cm thick. 1 of which is smooth, and the other is rough, evenly inlaid with 0.3 cm blockages spaced 1 cm apart to improve roughness. The fracture aperture was adjustable, and three sets of planes were made to obtain different fracture apertures with 2 and 3 mm. The sampling point was set at the end of the fracture.

Group 2 consists of a peristaltic pump, a single fracture, and a discharge flume. A flow distribution device was provided at the entrance of the fracture. A peristaltic pump was set as 55 and 100 rpm before each experiment of group 2. The hydraulic gradient was made constant by the overflow flume in every single experimental process. The sampling point was set at the center point of the fracture. The cross-section of the fracture was narrow, and the ratio of the length to the width of the fractures was large both in groups 1 and 2. Hence, vertical flow can be neglected and the one-dimensional transport was considered in this study.

Three smooth fractures with average apertures of 2, 3, and 4 mm were prepared by mounting two plexiglass plates horizontally. Each plate was set to be 60 cm long, 15 cm wide, and 1 cm thick. The single rough fracture was composed of a smooth plate and a rough plate formed by cutting a slot with 4 mm wide and 4 mm deep spaced every 4 mm on the side of the plexiglass plate. The space between the two plates is the same as that of the smooth fracture but the average apertures of rough fractures were 4, 5, and 6 mm.

In group 1, pattern 1 represents the granite fracture and pattern 2 represents the plexiglass fracture. And in group 2, pattern 1 represents the rough fracture and pattern 2 represents the fracture with a smooth surface. The summary of experimental models parameters can be seen in Table 1.

3.2. Experimental Procedure

The tracer experiments and simulations included the following steps which can also be seen briefly in Figure 3. The experimental apparatus was installed to ensure a watertight seal according to the design. The glass cement was used to prevent water leakage and stabilize the physical model. Water was continuously pumped into the fracture, and ensure that the air was expelled completely. The height of the recharge flume was adjusted by maintaining the height of the discharge flume to control the hydraulic head before injecting the BBF solution in group 1. When the flow rate stabilized, the BBF solution was injected through the injection syringe in an instant pulse while the chronograph begins to count. 5 ml of BBF solution with the concentration of 0.1 g/l was injected every single experiment. In group 2, the speed of the peristaltic pump was adjusted before adding the BBF solution and the dosage and concentration of the BBF solution were 2 ml and 0.15 g/l every time. Samples were taken at sampling points regularly till the BBF solution passed through the device completely. Each sampling volume was 3 ml to minimize quality loss effect, and the absorbance was determined to obtain the concentration. After each test, the plates were detached from the fracture and rinsed thoroughly with water. ADE and tFADE were conducted to simulate the transport process by adjusting the parameters of the models till the accuracy could not be improved.

4. Results and Discussion

4.1. Flow Fields and Transport Processes in the Single Fractures

The flow fields can deeply affect advection and mechanical dispersion. To quantify the flow regime and transport pattern, and are calculated by equations (4) and (5) and listed in Table 2. Noted that the for BBF is set as cm2/s according to Kone et al. [46]. The flow regime in this study is considered as the laminar flow because ranged from 5.0 to 11.6 which is obviously less than 100 [47]. The for the experiments ranges from 7693 to 19446 which are over 4000, indicating that molecular diffusion and transverse dispersion in the laboratory experiments can be negligible, and advection and mechanical dispersion in the flow direction made the main contribution [18, 48, 49].

To analyze the transport processes of the solute, the RTDs are calculated under different as equation (6). Figure 4 shows the RTDs for BBF transport in the granite fracture with the aperture of 1.71 mm (Figure 4(a)) and the RTDs for that in the smooth fracture with the aperture of 2 mm (Figure 4(b)).

When the solute transport follows Fick’s law, the BTC should present Gaussian distribution and the corresponding RTD will tend to be symmetric as shown in Figure 4(a). On the contrary, the tailing behaviors can be observed as the power-law drop in the RTDs if the non-Fickian transport occurs as shown in Figure 4(b).

If the solute transport process follows Fick’s law, the value of RTD after the peak will drop sharply rather than following the power law. To evaluate the degree of tailing behavior, the RTD value after the peak is fitted with power function of time and the is calculated. As seen in Figure 4(a), the three RTDs all present nearly inverted “U”-shaped curves indicating that the RTDs follow Gaussian distribution and the power-law tailing behavior does not happen. With the increase of Re and Pe, the of the power-law function increases. It indicates that the transport pattern is changing into non-Fickian transport. Fitting results with power function in Figure 4(b) show the obvious power-law decay of RTDs. The larger and , controlled by adjusting the flow velocity, lead to the larger fitted exponent increase which indicates the heavier tail. This phenomenon is consistent with the study of Dou et al. [50]. The increasing exponents may be due to the growth of the recirculation zones and eddy which lead to multirate exchange processes with the increasing [13]. The increase of and will lead to the heterogeneity of the flow field, thus affecting the migration of solute and making the tailing behavior more obvious.

The convection plays a dominant role in the transport process in this study, which is controlled by the velocity field. It has been revealed that the RTDs tend to be more symmetric and follow a Gaussian distribution with the increase of travel distance which leads to the scaling effect [37, 48]. The traveling distance and the structure of the velocity field may be the main cause which leads to the difference of BTCs in the two groups of experiments. It should be mentioned that time is taken as the -coordinate instead of the analyzed pore volume because the sampling method is used to measure the concentration of BBF in this study and the early arrival behavior was hard to be observed accurately. The tailing behavior as the typical characteristic of the subdiffusion was mainly discussed in this paper.

4.2. Model Comparison in the Weak Non-Fickian Transport Case

In order to evaluate the precision of the tFADE, the ADE and the tFADE are employed to fit the measured BTCs of BBF transport. Figure 5 shows the measured and fitted BTCs in group 1. To better observe the capacity of capturing the tailing behavior of the tFADE, logarithmic coordinates are set. A more gradual drop can be noticed in the BTC fitted by the tFADE than those fitted by the ADE. The relatively slow drop is controlled by the fractional-order . The declining delays the arrival of the peak at the same time. In the case of Fickian transport, the advantages of the tFADE are not obvious.

The associated values of RMSE and for the ADE and the tFADE are summarized in Table 3. The coefficient value is greater than 0.9, and RMSE is less than 0.1, illustrating the satisfactory precision of both the tFADE and the ADE in fitting BBF transport. The fitting errors of the tFADE are close to those of the ADE.

4.3. Model Comparison in the Anomalous Transport Case

To compare the performance of the tFADE in capturing the BBF transport in a single smooth or rough fracture, the BTCs for transport are fitted by the two models. The fitted BTCs for the fracture with an average aperture of 4 mm are shown in Figure 6, and the log charts are shown in the inserts. Tailing behaviors can be observed in the BTCs, and a more obvious phenomenon can be found in the rough fracture. The roughness elements enhance the retardation of BBF transport, which leads to subdiffusion. The difference between the fitting results of the tFADE and ADE is obvious. Though the ADE can capture the peak of the measured BTCs approximately, it cannot be used to capture the latter half of BTCs. In contrast, the tFADE can capture the trends of the measured BTCs overall and the tailing behaviors can be accurately described.

Besides the global error, the errors of the measured BTC are also calculated after the peak value using the two models. As shown in Table 4, the tFADE is more effective than the ADE in the rough fractures by comparing the and RMSE values. The tFADE is a precise model to capture solute transport in the fracture when non-Fickian transport happens especially in the fracture with complex morphology. Comparing the error in the phase of concentration decline, the precision of the tFADE is obviously improved for the rough fracture, which indicates the advantage of the tFADE in capturing the tailing behaviors.

The errors of the tFADE and ADE in the experiments of group 2 are presented in Figure 7. The circular and square markers represented the errors for the rough and smooth case, respectively. The dotted line is drawn as a reference for the accuracy comparison between the two models. Most of the points are distributed in the upper left part in the dotted line, indicating that the fitting accuracy of the tFADE is better than that of the ADE. As shown in Figure 7, the circular markers distribute on the left of the square markers indicating that the superiority of the tFADE is more obvious in capturing the transport process in the rough fracture. The subplot of Figure 7 presents RMSE comparison between the two models, and a similar phenomenon can be summarized as the major. The roughness elements affect the BBF transport by causing recirculation zones which enhance the heterogeneity of the flow field dictating the mass transfer [12, 51, 52]. The fractional operator has obvious advantages in capturing BBF transport in heterogeneous flow fields in the single fracture.

4.4. Fitting Parameters of the tFADE

The responses of fitting parameters of the tFADE to experimental conditions are summarized in Table 5. is the measured average flow velocity. is the best-fit flow velocity, and is the diffusion coefficient. In Table 5, is much less than in both two groups of experiments. It should be noted that there represents the average flow velocity and the transport velocity of the BBF can be affected by the local velocity variability which may be caused by fracture heterogeneity [53]. It may be the reason for the deviation of and . With the increase of measured flow velocity, also presents an increasing trend. As the pattern two of group 2 shown in Table 6, decreases with the increase of average aperture which results in the decreasing flow velocity with the increasing mean aperture. It indicates that decreases with the increasing aperture in the smooth single fracture. The same trend can be found in the rough fracture. It is because the flow rate is constant in each experiment of group 2. The increasing aperture causes a wider flow section which leads to the decrease of flow velocity.

Dispersion coefficient (), as a physical parameter to consider the heterogeneity of the fracture, is affected by the hydraulic aperture in the previous studies. However, the variation in this study was negligible (Tables 5 and 6). This is because is much lower than the flow velocity, and in our experiments is over in which case the transverse dispersion is negligible [8, 23]. Thus, it can be considered as a constant value for every single group of experiments.

The value of the fractional derivative order () is an important parameter to describe the anomalous diffusion. In Table 5, is insensitive to the variation of the measured flow velocity in both the granite and the plexiglass fractures in the range of 0.898 to 0.930. The inconspicuous range of the order indicated the unremarkable variation of the solute retention time. This is because the increasing is the cause of the growth of recirculation zones while the was small in these experiments [54]. In addition, the geometry structure of the fracture also affects the streamline and leads to the anomalous transport. As seen in pattern two of group 1, the value of decreased with the increase of the flow velocity which indicated the longer retention time of solute transport with the enhancing of the heterogeneity of the flow field.

The aperture is another factor affecting the value of . As seen in Table 6, the value of increased with the average aperture in pattern 2 of group 2, while the trend is not obvious in the pattern 1. It should be noted that the roughness of the fracture in pattern 2 is constant with the increase of the average aperture, while it decreases in pattern 1. When the aperture of the fractures increased, heterogeneity of the flow field attenuates and anomalous transport tends to approximate normal transport. The roughness of fractures contributes to secondary flows in the flow field such as eddies which partly interferes with the solute transport process [13, 55]. Hence, comparing to in the smooth case, the value of in the rough fracture in the same average aperture is relatively smaller overall indicating the complex transport process. The roughness elements were set with sharp corners in this work which caused resistance to flow and generated eddies to affect the transport process heavily [56].

5. Summary and Conclusions

In this work, two groups of experiments on BBF transport with different fracture media, relative roughness, flow rates, and fracture apertures were conducted. The potential cause of anomalous transport case was discussed briefly. The classical ADE and the tFADE were applied to capture the measured breakthrough curves (BTCs), and both Fickian and non-Fickian transports were conducted to evaluate the performance of the tFADE. The response of the fitting parameters and fitting errors of the tFADE was qualitatively analyzed. The main findings drawn from the experimental and fitting results are summarized as follows: (1)The experimental results showed that the more obvious tailing behavior occurred under the relatively short solute transport distance. Relatively larger enhanced the non-Fickian effect, and the measured RTDs demonstrated a notable power-law drop(2)The fitting results show that when non-Fickian transport is inconspicuous, both the ADE and tFADE can capture the transport process with satisfactory precision. However, when tailing behaviors which represent the non-Fickian transport occur, the classical ADE fails to capture the transport process while the tFADE can describe it well, especially in the reduction stage after the peak(3)An analysis of the physical meaning and response of the parameters of the tFADE indicates that the best-fit velocity presented a similar trend as the flow velocity increased. The best-fit velocity is notably smaller than the measured value. With the aperture increasing, the value of the fractional derivative order increased indicating the relatively low heterogeneity. The rough elements caused the eddy in the flow field and enhanced the heterogeneity, which lead to a smaller

The tFADE conducted in this study provides a precise approach for capturing the solute transport process in single fractures. The contribution of this study is to provide understanding of the applicability and parameter responses of the tFADE which help better understand the physical meaning of the model. Since artificial roughness was considered in this study, the performance of the tFADE in rougher fractures and the sensibility of the fractional order to the roughness should be further investigated. The limitation of this study is that the flow velocity was considered as stable, while in reality, it might be variable temporally.


:The order of fractional differentiation (−)
:The concentration of the solute ()
:The dimensionless injected concentration (−)
:The dimensionless concentration under pulse injection condition (−)
:The dimensionless concentration under step injection condition (−)
:The apparent fracture aperture ()
:The dispersion coefficient ()
:The molecular diffusion coefficient ()
:The dynamic fluid viscosity ()
:Peclet number (−)
:The flow rate ()
:The coefficient of determination (−)
:Reynolds number (−)
:The root mean square error (−)
:The fluid density ()
:The time ()
:The duration of solution injection ()
:The dimensionless time (−)
:The average flow velocity ()
:The advective velocity ()
:The volume of the fracture ()
:The spatial coordinate of the solute along the direction of flow ().

Data Availability

The data that support the findings of this study are available from the corresponding author, Yong Liu, upon reasonable request.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


This work was supported by the National Natural Science Foundation of China (nos. 41831289, 41877191, and 42072276) and the Shandong Province Science and Technology Support Pilot Program for Achieving Water Ecological Civilization under Contract no. SSTWMZCJHSD04.