Effect of Roughness on Conservative Solute Transport through Synthetic Rough Single Fractures

Understanding solute transport in fractured rocks is of particular importance in many applications. Aperture values ranging from 4.7 to 8.7 mm and Reynolds number (Re) values at 9.38~1743.8 were set for investigating fluid flow through synthetic horizontal single smooth and rough fractures. The Brilliant Blue FCF dye was chosen as the tracer to visualize the transport process. This paper focuses on the dispersion process in rough single fractures under non-Darcian flow conditions. Non-Darcian flow existed in both smooth and rough single fractures and the average flow velocity–hydraulic gradient (V–J) relationships were best described by the Forchheimer equation. The main objectives were to check the existing flow and transport models and to study possible correlations between fitting parameters and heterogeneities. The classical advection dispersion equation (ADE) model failed to capture the long-tailing of breakthrough curves (BTCs). Instead, the continuous time random walk (CTRW) model was better at explaining BTCs in both smooth and rough fractures, especially in capturing the long-tailing feature. The non-Darcian coefficient βc in the Forchheimer equation and the coefficient β in the CTRW model appeared to be most relevant for characterizing the heterogeneity of the rough single fractures.


Introduction
The management of groundwater resources and control of contaminated aquifers require an understanding of the processes of flow and transport in porous or fractured rocks [1][2][3][4][5].Fractures and bedding planes or faults give rise to preferential flow paths for groundwater, pollution in dilution, and free product, and thus are of great concern [6].
A single fracture has been traditionally idealized as parallel plates to obtain a tractable description of fluid flow and solute transport, and the model used to describe fluid flow in such an idealized single fracture is the local cubic law (LCL), which is essentially the expression of Darcy's law for a single fracture [7].Also, Fickian transport is believed to be the "right" form of governing law for transport in the subsurface, where the dispersive mass flux is assumed to be proportional to the first derivative of the resident solute concentration [8].However, real fractures have rough walls with points of contact, in which the transport is found to be non-Fickian on many occasions [1,2,7].In general, there are two distinctive features of non-Fickian transport.First, the peak value of a breakthrough curve (BTC) originated from a pulse input of a contaminant or tracer arrives earlier than expected from the advection dispersion equation (ADE) (early arrival).The early arrival often suggests one or multiple preferential paths for the solute transport.Second, the tailing of the BTC lasts much longer than expected (the long tail) [2].The long tail phenomena can possibly be interpreted in terms of: a multi-rates mass exchange process between mobile and immobile domains [9,10]; fracture wall roughness [11]; the absorption of solute on fracture walls [12,13]; a permeable rock matrix [14]; and the entrapment of solute in eddies inside the fracture [15][16][17].
Bauget and Fourar [1] and Nowamooz et al. [27] studied the solute transport in a transparent replica of a real single fracture and found that the Fickian ADE was incapable of modeling the heavy (or long) tails' behavior.Alternatively, CTRW was able to describe the non-Fickian dispersion.The influence of surface roughness on nonlinear flow behaviors in three-dimensional (3D) self-affine rough single fractures was determined by Wang et al. [28] using Lattice Boltzmann simulations.Non-Fickian transport through two-dimensional (2D) rough single fractures was numerically studied by Wang and Cardenas [29], and their results clearly showed early arrival and heavy tailing in BTCs.The degree of deviation of transport from Fickian to non-Fickian is captured by the parameter β of the truncated power law used in CTRW.Specifically, β was found to be proportional to fracture heterogeneity [1,24].
On the other hand, non-Darcian flow, which can arise in a number of different ways, may have a great impact on solute transport.Non-Darcian flow is particularly prone to occuring in heterogeneous geological formations such as in a single fracture [11,[30][31][32] or in a fractured network [33] due to relatively high speeds of flow.
To investigate solute transport under non-Darcian flow conditions, Qian et al. [8] established a well-controlled physical model with a vertical smooth parallel fracture, and the non-Fickian transport process was identified with early arrival and long tails.A MIM model fitted both peak and tails of the observed BTCs better than the ADE model.On the basis of this experience, with the purpose of describing solute transport under different flow velocities and fracture apertures, Chen et al. [34] carried out a series of flow and tracer test experiments on an artificial channeled single fracture-a single fracture with contact in certain areas-constructed in the laboratory.The flow condition showed a non-Darcian feature (best described by the nonlinear Forchheimer equation) and BTCs showed a non-Fickian feature such as early arrival of the peak value, long tailing and multi-peak phenomena.The results of this study [34] showed that ADE was not adequate to describe BTCs.
In a similar way to that undertaken by the study of Chen et al. [34], laboratory flow and tracer tests have been carried out on an artificially created fractured rock sample by Cherubini et al. [7], which showed a Forchheimer type of non-Darcian flow, and profound mobile-immobile mass exchange in terms of transport.
The main goals here were to test different models on the solute transport through a rough single fracture under non-Darcian flow conditions and to study possible correlations between fitting parameters and fracture heterogeneities.With these objectives in mind, the following tasks were carried out in sequence.First, synthetic rough single fractures with different roughness elements were designed.Then a series of flow and tracer test experiments were performed with a Reynolds number (Re) ranging from 9.38 to 1743.8.An imaging process was introduced to obtain a precise 2D concentration distribution for visual inspection.Subsequently, the characteristics of BTCs were analyzed and two models (ADE and CTRW) tested for their capacity to characterize the non-Fickian transport behavior under non-Darcian flow conditions.

Flow Model
Generally the model used to describe fluid flow in fractured media is LCL, which adapts Darcy's law [35]: where J is the hydraulic gradient [dimensionless], K is the hydraulic conductivity [m/s], µ is dynamic viscosity [mPa s], ρ is water density [kg/m 3 ], g is the gravitational acceleration [m/s 2 ], V is flow velocity [m/s] and k is the Darcian coefficient.A minor note is that porosity in a single fracture is expected to be 1, meaning that there is no filling materials inside the single fracture.When a nonlinear flow feature occurs, the flow models commonly used to represent non-Darcian flow behavior are the Forchheimer law and the Izbash law.The former includes a quadratic term of velocity to represent the inertial effect [36,37]: where β c [m −1 ] is the inertial term friction coefficient or non-Darcian coefficient, and a [s/m] and b c [s 2 /m 2 ] are the linear and inertial coefficients, respectively.If b c = 0, flow becomes Darcian.If a = 0, flow becomes fully developed turbulent [30].
The Izbash law, or the power-law, equation is as follows [35]: where c [s 2 /m 2 ] and m [dimensionless] are two constant coefficients.The Forchheimer and Izbash equations are equivalent when a = 0 in Equation ( 1) and m = 2 in Equation (2).Another dimensionless parameter to analyze the inertial versus viscous forces is the Reynolds number (R e ), which is defined as [30]: where b is the fracture aperture [mm], and ν is the kinematic viscosity of water [mm 2 /s] (here ν = 1.308 mm 2 /s at 10 • C).The V-J relationships obtained by experiments were fitted by the utility Cftool in Matlab.The quality of the fitting was discussed based on the decision coefficient r 2 , the error squared SSE, and the root mean square error (RMSE), as: where J i is the ith hydraulic gradient, J i is the average hydraulic gradient, J ie is fitted hydraulic gradient, J i0 is the test hydraulic gradient, J i0 is the average test hydraulic gradient, N is the number of the test hydraulic gradients.

The Advection Dispersion Approach
The general formation of the well-known ADE for a conservative solute without any sink/sources in a one-dimensional framework is given as, where C is the average solute concentration, D is the longitudinal hydrodynamic dispersion coefficient, t is time, and x is distance from the injected position along the flow direction.An analytical solution of Equation ( 8) with an instantaneous injection of tracer is [38]: where M is the injected mass, A is the cross-sectional area of the fracture over which injection occurs.
The tracer is assumed to be injected uniformly over the entire saturated thickness of the fracture.

Continuous Time Random Walk
The continuous time random walk theory (CTRW) has been developed specifically to model the non-Fickian tracer transport [1,[39][40][41][42].The CTRW transport equation is based on the Fokker-Planck equation with a memory equation (FPME) which originally describes the temporal evolution of the probability density function of the velocity of a particle under the influence of drag and random forces.In CTRW, the Laplace transformed concentration function C(x, w), is given by the one-dimensional form of FPME, where w is the Laplace variable, and u ψ and D ψ are transport velocity and dispersion in the framework of CTRW respectively.In Equation (10), the memory function M(w) captures the anomalous or non-Fickian transport induced by local heterogeneity or process and the corresponding formulation is given by where t is some characteristic time and ψ(w) is the transition rate probability, which is the core of the CTRW model.In general, there are three different models for ψ(w), the asymptotic model, the truncated power law (TPL) model, and the modified exponential model.Details for the three different models of ψ(w) can be found in the reference Cortis and Berkowitz [18].Here, the most-used TPL model was adopted in the framework of CTRW, where t 1 represents the lower limit time when the power law behavior begins, t 2 is the cut-off time describing the Fickian behavior, and Γ() represents the incomplete Gamma function.The parameter β indicates different types of anomalous transport.Following Berkowitz et al. [40], ψ(w) may be approximated by a power-law decay function, t −1−β , where β is a parameter that characterizes dispersion.For β > 2, the process is Fickian and CTRW is equivalent to the classical ADE model.If 2 > β > 1, the process is no longer Fickian and shows moderately anomalous transport; If 1 > β > 0, the process indicates highly anomalous transport.

Experiment Setup and Test Process
Artificial horizontal single fractures were constructed in the laboratory using organic glass plates, as per the schematic shown in Figure 1a.Although such artificial single fractures may not represent all the features of natural fractures, excellent control of flow field can be imposed on the apparatus, and the influence of roughness on water flow through single fractures can be studied in great detail.A detailed setup can be found in Chen et al. [11].Three different kinds of roughness were used for the upper wall (Figure 1b): a smooth parallel plate; a rough plate with rectangular rough elements; and a rough plate with trapezoidal rough elements.For each type of rough single fractures, four types of asperity heights (∆), including 2 mm, 4 mm, 6 mm and 8 mm, were used.The corresponding relative roughness (ε), which equals the ratio of ∆ and the maximum aperture b max (10 mm), were 0.2, 0.4, 0.6 and 0.8 respectively.Therefore, the average apertures (b a ), which equal the average asperity heights of the rough single fracture (calculation shown in Figure 1b), were obtained and shown in Table 1.
Water 2017, 9, 656 5 of 17 great detail.A detailed setup can be found in Chen et al. [11].Three different kinds of roughness were used for the upper wall (Figure 1b): a smooth parallel plate; a rough plate with rectangular rough elements; and a rough plate with trapezoidal rough elements.For each type of rough single fractures, four types of asperity heights (  ), including 2 mm, 4 mm, 6 mm and 8 mm, were used.The corresponding relative roughness (  ), which equals the ratio of  and the maximum aperture bmax (10 mm), were 0.2, 0.4, 0.6 and 0.8 respectively.Therefore, the average apertures (ba), which equal the average asperity heights of the rough single fracture (calculation shown in Figure 1b), were obtained and shown in Table 1.The water levels at the recharge and discharge reservoirs were measured using piezometers with errors less than 0.5 mm.The hydraulic gradient through the single fracture was controlled through adjusted water levels between the inflow and the outflow reservoirs.The average flow rate was monitored using a calibrated flow meter LZB-25 (glass Rotameter made by Chengfeng Flowmeter Co., Ltd.Changzhou City, Jiangsu Prov.China) with a maximum permissible error of 1.5%.The average velocity was calculated by the average flow rate divided by the average aperture and width  The water levels at the recharge and discharge reservoirs were measured using piezometers with errors less than 0.5 mm.The hydraulic gradient through the single fracture was controlled through adjusted water levels between the inflow and the outflow reservoirs.The average flow rate was monitored using a calibrated flow meter LZB-25 (glass Rotameter made by Chengfeng Flowmeter Co. Ltd., Changzhou, Jiangsu, China) with a maximum permissible error of 1.5%.The average velocity was calculated by the average flow rate divided by the average aperture and width of the fracture.The flow experiments under different hydraulic gradients were carried out for a given single fracture.After that, the upper plate was replaced by a plate with a different type of roughness, and the entire flow experiment was repeated.
The dye Brilliant Blue FCF F (bis {4-(N-ethyl-N-3-sulfophenylmethyl) aminophenyl}-2-sulfophenyl methylium disodium salt) was chosen as the tracer here for visual inspection and image process.It is important to choose an appropriate tracer concentration.If the selected concentration is too high, the gravity and density of the tracer may be a concern.If the concentration is too low, it may not generate images of tracer transport that are of sufficiently high quality for visual inspection.Concentration of the dye Brilliant Blue FCF was set as 1200 mg/L after a series of preliminary tests.Instantaneous tracer was adopted using an injection syringe and 1 mL tracer was injected in each test.The tracer injection point was illustrated in Figure 1a.
The digital camera (Canon EOS 500D) was fixed on a tripod.The length of the camera lens from the main body of the single fracture was set about 650 mm.A light source was placed under the fracture and the light emitted by the lamp was evenly distributed over the transparent plexiglas plate of the main fracture.All tracer tests were carried out in a home-made digital darkroom and the only light source during the test was provided by the lamp.

Image Process
Image analysis was improved to the extent that estimation of dye concentration from soil color was possible [43].The main purpose here was to analyze the color information in the image to obtain the mathematical relationship between the color and concentration, and to further analyze the temporal and spatial variation of the concentration.The main factors that affect the color of the image during the test are the intensity of the light source, the color temperature, and the camera's parameter settings (mainly the exposure and white balance).In order to ensure the accuracy of the digital image's color, color temperature adjustment is needed for each photo.
In this study, the water-filled fracture image was used as the background of each tracer test.The background image was subtracted from the image containing the tracer using the image subtract code by Matlab software.The result is shown in Figure 2. The residual color of the fracture surface caused by the dye tracer was reduced.
The dye Brilliant Blue FCF F (bis {4-(N-ethyl-N-3-sulfophenylmethyl) aminophenyl}-2sulfophenyl methylium disodium salt) was chosen as the tracer here for visual inspection and image process.It is important to choose an appropriate tracer concentration.If the selected concentration is too high, the gravity and density of the tracer may be a concern.If the concentration is too low, it may not generate images of tracer transport that are of sufficiently high quality for visual inspection.Concentration of the dye Brilliant Blue FCF was set as 1200 mg/L after a series of preliminary tests.Instantaneous tracer was adopted using an injection syringe and 1 mL tracer was injected in each test.The tracer injection point was illustrated in Figure 1a.
The digital camera (Canon EOS 500D) was fixed on a tripod.The length of the camera lens from the main body of the single fracture was set about 650 mm.A light source was placed under the fracture and the light emitted by the lamp was evenly distributed over the transparent plexiglas plate of the main fracture.All tracer tests were carried out in a home-made digital darkroom and the only light source during the test was provided by the lamp.

Image Process
Image analysis was improved to the extent that estimation of dye concentration from soil color was possible [43].The main purpose here was to analyze the color information in the image to obtain the mathematical relationship between the color and concentration, and to further analyze the temporal and spatial variation of the concentration.The main factors that affect the color of the image during the test are the intensity of the light source, the color temperature, and the camera's parameter settings (mainly the exposure and white balance).In order to ensure the accuracy of the digital image's color, color temperature adjustment is needed for each photo.
In this study, the water-filled fracture image was used as the background of each tracer test.The background image was subtracted from the image containing the tracer using the image subtract code by Matlab software.The result is shown in Figure 2. The residual color of the fracture surface caused by the dye tracer was reduced.
The RGB color model is the most commonly used on digital cameras.We tested the relationship of red (R), green (G), and blue(B) values and the concentration and finally found that the R value in the RGB model had the strongest correlation with the concentration [44].The polynomial fitting of the relationship between the R value of the image and the concentration C was performed with a coefficient of determination of r 2 = 0.990 and RMSE = 0.757.The high r 2 value and the relatively low RMSE value suggested that the image processing method was acceptable.
Different fracture thickness and roughness for each test meant that the color of the same The RGB color model is the most commonly used on digital cameras.We tested the relationship of red (R), green (G), and blue(B) values and the concentration and finally found that the R value in the RGB model had the strongest correlation with the concentration [44].
The polynomial fitting of the relationship between the R value of the image and the concentration C was performed with a coefficient of determination of r 2 = 0.990 and RMSE = 0.757.The high r 2 value and the relatively low RMSE value suggested that the image processing method was acceptable.
Different fracture thickness and roughness for each test meant that the color of the same concentration solution in the medium was different.Therefore, the C-R relationships for each roughness type were obtained and the results were shown in Table 2, where the concentration has a unit of mg/L.

Roughness Element C-R Relation Equation r 2
Rectangular A

Flow and Tracer Experiments
The present works focused on relatively low velocity and small aperture.In order to investigate the tracer test in non-Darcian flow conditions, we chose a medium aperture ranging from 0.0047 to 0.0087 m and the Reynolds number (R e ) at 9.38~1743.8.The evidence of non-Darcian flow is reflected in the nonlinear relationship between the hydraulic gradient and the flow rate at the given range of R e .
To analyze the characteristics of solute transport through single fractures more accurately, instantaneous tracer tests using the dye Brilliant Blue FCF were performed.The horizontal fractures with smooth parallel plates, rectangular cross-section rough fractures, and trapezoidal cross-section rough fractures, were chosen.For demonstration purposes, concentrations at position x = 555 mm downstream from the point of injection, which is located at the middle elevation at x = 0 mm, were calculated under flow velocity of 9.7 mm/s and 51.7 mm/s for smooth parallel plates; 44.3 mm/s and 71.0 mm/s for rectangular cross-section rough fractures; and 20.5 mm/s and 51.2 mm/s for trapezoidal cross-section rough fractures; and the results can be found in Figures 3-5 respectively.
A few observations can be made about Figure 3.For the smooth parallel plates shown in Figure 3, when flow velocity V was as small as 9.7 mm/s, the solute plume moved forward as a lump, but the color inside the plume for each point was not the same, indicating that the solute had not yet fully mixed.The front of the plume showed a typical Gaussian distribution.Lateral dispersion was also clearly visible and irregular movement of the solute with a long tail occurred.When the flow velocity V reached a relatively higher value of 55.1 mm/s, the front of the plume also showed a Gaussian distribution but the lateral dispersion was reduced.Similarly, with that for V of 9.7 mm/s, the concentration inside the plume showed apparent difference and a distinctive long tail also appeared in this case.
Figure 4 showed the solute plume of Brilliant Blue FCF through a rough single fracture with trapezoidal roughness.The differences from that shown in Figure 3 were as follows: the front of the plume was irregular, rather than smooth, which was obviously caused by the heterogeneity of fracture Water 2017, 9, 656 8 of 17 wall roughness.Figure 4 also indicates that the solute concentration inside the plume was different and the solute migration rate was not the same for each point.In addition, retention of solute in some "dead-ends" caused by fracture wall roughness was clearly visible.Specifically, the retained solute was trapped for a long time in those dead-ends, leading to a long tail.
Figure 5 shows the solute plume through a rough single fracture with rectangular roughness.It is similar to Figure 4, with the exception that the irregularity degrees of plume were greater and the long tailing and solute-retention capacity were more pronounced.
The present works focused on relatively low velocity and small aperture.In order to investigate the tracer test in non-Darcian flow conditions, we chose a medium aperture ranging from 0.0047 to 0.0087 m and the Reynolds number (Re) at 9.38~1743.8.The evidence of non-Darcian flow is reflected in the nonlinear relationship between the hydraulic gradient and the flow rate at the given range of Re.
To analyze the characteristics of solute transport through single fractures more accurately, instantaneous tracer tests using the dye Brilliant Blue FCF were performed.The horizontal fractures with smooth parallel plates, rectangular cross-section rough fractures, and trapezoidal cross-section rough fractures, were chosen.For demonstration purposes, concentrations at position x = 555 mm downstream from the point of injection, which is located at the middle elevation at x = 0 mm, were calculated under flow velocity of 9.7 mm/s and 51.7 mm/s for smooth parallel plates; 44.3 mm/s and 71.0 mm/s for rectangular cross-section rough fractures; and 20.5 mm/s and 51.2 mm/s for trapezoidal cross-section rough fractures; and the results can be found in Figures 3-5        A few observations can be made about Figure 3.For the smooth parallel plates shown in Figure 3, when flow velocity V was as small as 9.7 mm/s, the solute plume moved forward as a lump, but the color inside the plume for each point was not the same, indicating that the solute had not yet fully mixed.The front of the plume showed a typical Gaussian distribution.Lateral dispersion was also clearly visible and irregular movement of the solute with a long tail occurred.When the flow velocity V reached a relatively higher value of 55.1 mm/s, the front of the plume also showed a Gaussian distribution but the lateral dispersion was reduced.Similarly, with that for V of 9.7 mm/s, the concentration inside the plume showed apparent difference and a distinctive long tail also appeared in this case.
Figure 4 showed the solute plume of Brilliant Blue FCF through a rough single fracture with trapezoidal roughness.The differences from that shown in Figure 3 were as follows: the front of the plume was irregular, rather than smooth, which was obviously caused by the heterogeneity of fracture wall roughness.Figure 4 also indicates that the solute concentration inside the plume was different and the solute migration rate was not the same for each point.In addition, retention of solute

Interpretation of Experiments
In the following section, we seek to interpret the obvious non-Fickian transport phenomenon observed in Figures 3-5 using an appropriate theoretical model.Before discussing the transport process, it is necessary to address the flow process first as it was the fundamental driving force for hydrodynamic dispersion.

The Non-Darcian Equations
The experimental flow velocity-hydraulic gradient (V-J) relationships for different rough single fractures under different flow conditions are illustrated in Figures 6-8 and the fitting results, using linear Equation (1) and nonlinear Equations ( 2) and (3), and goodness of fit were also shown in Tables 3-5.Fitting parameters are listed in Table 6.
Water 2017, 9, 656 9 of 17 in some "dead-ends" caused by fracture wall roughness was clearly visible.Specifically, the retained solute was trapped for a long time in those dead-ends, leading to a long tail.
Figure 5 shows the solute plume through a rough single fracture with rectangular roughness.It is similar to Figure 4, with the exception that the irregularity degrees of plume were greater and the long tailing and solute-retention capacity were more pronounced.

Interpretation of Experiments
In the following section, we seek to interpret the obvious non-Fickian transport phenomenon observed in Figures 3-5 using an appropriate theoretical model.Before discussing the transport process, it is necessary to address the flow process first as it was the fundamental driving force for hydrodynamic dispersion.

The Non-Darcian Equations
The experimental flow velocity-hydraulic gradient (V-J) relationships for different rough single fractures under different flow conditions are illustrated in Figures 6-8 and the fitting results, using linear Equation (1) and nonlinear Equations ( 2) and ( 3), and goodness of fit were also shown in Tables 3-5.Fitting parameters are listed in Table 6.Table 3. Fitting results and goodness of fit for the V-J relationship for flow through single fractures with smooth parallel plates.
The fitting parameters in Table 6 exhibited the following patterns.First, with the decrease of the asperity heights of the rough elements from 8 mm to 2 mm, the relative roughness (ε) decreased from 0.8 to 0.2 and the average aperture increased from 4.7 mm to 8.7 mm, the value of k in the Darcian equation decreased, indicating an increasing of permeability for the fracture.The value of c in the Izbash equation also decreased with the average aperture, but the value of m in the Izbash equation did not show any obvious trend.Second, fitting parameters showed different trends for smooth and rough single fractures.For smooth single fractures, the value of a in the Forchheimer equation decreased with the average aperture but the value of b c and the value of the non-Darcian coefficient β c increased with the average aperture.Third, the nonlinear phenomenon was more striking in smooth single fractures with larger average apertures.However, for rough single fractures in which flow can be described by the Forchheimer equation, a increased with the average aperture but b c and β c decreased with the average aperture.This was probably because the increase of asperity heights of rough elements in a rough single fracture would create greater fluid-solid contact surface, and enhance the frictional force for flow.Therefore, the relative importance of viscous flow versus inertial flow for a greater relative roughness is increased, reducing the flow nonlinearity (which is primarily caused by the inertial effect of flow).

The ADE and CTRW Models
The molecular diffusion coefficient D m is set as 5.0 × 10 −11 m 2 •s −1 based on the Taylor-Aris experiment [1].The Peclet numbers (P e = bV/D m ) in our experiment were greater than 5 × 10 6 .As pointed out by Bauget and Fourar [1], when the Peclet number is greater than 4000, transverse dispersion can be ignored.The P e numbers of our experiments were more than 100 times greater than 4000, which warranted the use of the one-dimensional approach.The CTRW Matlab Toolbox V3.1 [20] was employed for the CTRW TPL model.
BTCs at position x = 555 mm in single fractures with smooth parallel plates and rectangular rough single fractures were illustrated and fitted by ADE and CTRW TPL models.The results are shown in Figures 9 and 10 and Tables 7 and 8.

The ADE and CTRW Models
The molecular diffusion coefficient Dm is set as 5.0 × 10 −11 m 2 •s −1 based on the Taylor-Aris experiment [1].The Peclet numbers (Pe = bV/Dm) in our experiment were greater than 5 × 10 6 .As pointed out by Bauget and Fourar [1], when the Peclet number is greater than 4000, transverse dispersion can be ignored.The Pe numbers of our experiments were more than 100 times greater than 4000, which warranted the use of the one-dimensional approach.The CTRW Matlab Toolbox V3.1 [20] was employed for the CTRW TPL model.
BTCs at position x = 555 mm in single fractures with smooth parallel plates and rectangular rough single fractures were illustrated and fitted by ADE and CTRW TPL models.The results are shown in Figures 9 and 10 and Tables 7 and 8.
We can make several interesting observations from Figures 9 and 10.First, BTCs followed the normal distribution in smooth single fractures, especially when the flow velocity was small.Second, as the flow velocity V increased, BTCs appeared to exhibit non-normal distributions with long tails.These non-normal distributions could be found ubiquitously in rough single fractures.Third, the classical ADE model was incapable of capturing the long-tailing of BTCs either in smooth or rough single fractures.However, the CTRW TPL model could better explain BTCs in both smooth and rough single fractures, especially in capturing the long-tailing phenomenon.
Similar conclusions can be obtained from Tables 7 and 8     It can be seen from Tables 7 and 8 that 1 < β < 2, reflecting how the solute transport in both smooth and rough single fractures under non-Darcian flow conditions had not reached Fickian migration.The β value decreased with the flow velocity and the relative roughness (  ), indicating that the degree of nonlinearity was stronger with larger flow velocity; and the degree of deviation from Fickian transport was also higher.Meanwhile, an increase of the asperity height would further aggravate the irregularity of solute transport.The above observations were made based on a β value of 1.05 and a  value of 0.6 in single fractures with rectangular roughness.We can make several interesting observations from Figures 9 and 10.First, BTCs followed the normal distribution in smooth single fractures, especially when the flow velocity was small.Second, as the flow velocity V increased, BTCs appeared to exhibit non-normal distributions with long tails.These non-normal distributions could be found ubiquitously in rough single fractures.Third, the classical ADE model was incapable of capturing the long-tailing of BTCs either in smooth or rough single fractures.However, the CTRW TPL model could better explain BTCs in both smooth and rough single fractures, especially in capturing the long-tailing phenomenon.
Similar conclusions can be obtained from Tables 7 and 8, as demonstrated by relatively smaller values of R 2 and larger RMSE values for the ADE model compared with the CTRW TPL model.The R 2 values were greater than 0.95 and the RMSE values were less than 0.1 with the CTRW TPL model, indicating satisfactory goodness of fit.It could be seen that the t 1 value in the CTRW TPL model was less than the observation time of the test and the t 2 value in the same model was much larger than the observation time of the test, indicating that the solute transport through single fractures under experimental flow conditions was far from reaching Fickian transport.
It can be seen from Tables 7 and 8 that 1 < β < 2, reflecting how the solute transport in both smooth and rough single fractures under non-Darcian flow conditions had not reached Fickian migration.The β value decreased with the flow velocity and the relative roughness (ε), indicating that the degree of nonlinearity was stronger with larger flow velocity; and the degree of deviation from Fickian transport was also higher.Meanwhile, an increase of the asperity height would further aggravate the irregularity of solute transport.The above observations were made based on a β value of 1.05 and a ε value of 0.6 in single fractures with rectangular roughness.
In addition, the parameter v ψ , which is the average velocity of solute particles in the CTRW TPL model, was much larger than the actual average velocity of V.This implied that the average velocity of solute particles was greater than the average velocity of flow.
To further study the effect of relative roughness (ε) on solute transport, BTCs at position x = 355 mm downstream, V of 71.0 m/s, b of 8.7 mm and 6.0 mm, and ε of 0.2 and 0.6 in single fractures with rectangular roughness were fitted by the CTRW TPL model; while BTCs at position x = 355 mm downstream, V of 49.5 m/s, b of 8.7 mm and 7.3 mm, and ε of 0.2 and 0.4 in single fractures with trapezoidal roughness were also fitted by the CTRW TPL model.The fitting results can be seen in Table 9.As a reference for comparison, BTCs at position x = 355 mm downstream, V of 59.0 m/s, and b of 4.7 mm and 8.7 mm in single fractures with smooth parallel plates were also fitted by the CTRW TPL and the results are also listed in Table 9.From Table 9, one can see that an increase of ε resulted in a significant reduction in the value of β; that is, the degree of non-Fickian transport (or the deviation from Fickian transport) was affected significantly by ε in rough single fractures.However, the change of the b values in single smooth fractures had little effect on the value of β.Furthermore, to analyze the effect of fracture roughness type (rectangular or trapezoidal) on solute transport, BTCs from rough single fractures with rectangular and trapezoidal roughness under the same flow condition (V = 36.0mm/s) and the same ε value of 0.4 were fitted by the CTRW TPL model, and the results are shown in Table 10.As a reference for comparison, BTCs from smooth single fractures and rough single fractures with trapezoidal roughness under the same flow condition (V = 51.5 mm/s) and the same b value of 6 mm were fitted by the CTRW TPL model and the results are also listed in Table 10.We can conclude from Table 10 that the β value under the same flow condition in a rough single fracture with a rectangular-type of roughness was smaller than that with a trapezoidal-type of roughness.As shown in Chen et al. [11], recirculation of flow in the eddies and curved streamline in a rough single fracture with the rectangular-type of roughness caused greater energy loss than that in the rough single fracture with a trapezoidal-type roughness.Such additional energy loss caused a greater degree of flow nonlinearity or stronger non-Darcian flow, leading to stronger non-Fickian transport and a smaller β value.
In summary, the heterogeneity of the fractured media, which is reflected in the relative roughness and roughness structure, was primarily responsible for the non-Fickian transport through rough single fractures.Also, the increase of flow nonlinearity of non-Darcian flow can exacerbate the non-Fickian transport phenomenon.

Conclusions
Artificial smooth and rough single fractures were constructed in the laboratory using organic glass plates.Compared with previous works focused on relatively low velocity and small apertures, medium aperture values ranging from 4.7 to 8.7 mm and the Reynolds numbers of 9.38~1743.8were set to characterize the tracer transport through smooth and rough single fractures under non-Darcian flow conditions.The dye Brilliant Blue FCF was chosen as the tracer to facilitate the image process, enabling the solute transport process to be directly observed more intuitively.
It was found that non-Darcian flow existed in both smooth and rough single fractures under these experimental conditions.The nonlinear flow phenomenon was more obvious in rough single fractures, especially with large asperity heights.The non-Darcian flow velocity and hydraulic gradient (V-J) relationships were best described by the Forchheimer equation.
The classical ADE model was incapable of capturing the long-tailing of BTCs either in smooth or rough single fractures; the CTRW TPL model, with more adjustable parameters, can explain BTCs both in smooth and rough single fractures better, especially when it comes to capturing the long-tailing phenomenon.
The heterogeneity of the fractured media, which was manifested through the relative roughness and roughness structure, had a significant effect on the degree of deviation from Fickian transport (or the degree of non-Fickian transport) through rough single fractures.Also, the increase in the flow nonlinearity can further exacerbate the non-Fickian transport phenomenon.The non-Darcian coefficient β c and the coefficient β in the CTRW TPL model seem to be good parameters for characterizing the heterogeneity of the rough fracture.The results may enhance our understanding of solute transport under non-Darcian flow conditions through fractured media.The study of solute transport and reactive transport through a real single fracture under non-Darcian conditions needs to be studied further.

Figure 1 .
Figure 1.Schematic figure of the experimental setup.

Figure 1 .
Figure 1.Schematic figure of the experimental setup.

Figure 2 .
Figure 2. Difference in value between the adjusted image and the background image.

Figure 2 .
Figure 2. Difference in value between the adjusted image and the background image. respectively.

Figure 3 .
Figure 3. Solute plume of the tracer test in single fractures with smooth parallel plates.Figure 3. Solute plume of the tracer test in single fractures with smooth parallel plates.

Figure 3 .
Figure 3. Solute plume of the tracer test in single fractures with smooth parallel plates.Figure 3. Solute plume of the tracer test in single fractures with smooth parallel plates.Water 2017, 9, 656 8 of 17

Figure 4 .
Figure 4. Solute plume of the tracer test in rough single fractures with rectangular roughness.Figure 4. Solute plume of the tracer test in rough single fractures with rectangular roughness.

Figure 4 .
Figure 4. Solute plume of the tracer test in rough single fractures with rectangular roughness.Figure 4. Solute plume of the tracer test in rough single fractures with rectangular roughness.

Figure 4 .
Figure 4. Solute plume of the tracer test in rough single fractures with rectangular roughness.

Figure 5 .
Figure 5. Solute plume of the tracer test in rough single fractures with trapezoidal roughness.

Figure 5 .
Figure 5. Solute plume of the tracer test in rough single fractures with trapezoidal roughness.

Figure 6 .
Figure 6.The V-J relationship and fitting results for flow through single fractures with smooth parallel plates.

Figure 6 .
Figure 6.The V-J relationship and fitting results for flow through single fractures with smooth parallel plates.

Figure 6 .
Figure 6.The V-J relationship and fitting results for flow through single fractures with smooth parallel plates.

Figure 7 .
Figure 7.The V-J relationship and fitting results for flow through single fractures with rectangular roughness.

Figure 7 . 17 Figure 8 .
Figure 7.The V-J relationship and fitting results for flow through single fractures with rectangular roughness.Water 2017, 9, 656 10 of 17

Figure 8 .
Figure 8.The V-J relationship and fitting results for flow through single fractures with trapezoidal roughness.
, as demonstrated by relatively smaller values of R 2 and larger RMSE values for the ADE model compared with the CTRW TPL model.The R 2 values were greater than 0.95 and the RMSE values were less than 0.1 with the CTRW TPL model, indicating satisfactory goodness of fit.It could be seen that the t1 value in the CTRW TPL model was less than the observation time of the test and the t2 value in the same model was much larger than the observation time of the test, indicating that the solute transport through single fractures under experimental flow conditions was far from reaching Fickian transport.

Figure 9 .
Figure 9. BTCs and fitting results in single fractures with smooth parallel plates (x = 555 mm).

Figure 9 .
Figure 9. BTCs and fitting results in single fractures with smooth parallel plates (x = 555 mm).

Figure 10 .
Figure 10.BTCs and fitting results in single fractures with rectangular rough fractures (x = 555 mm).

Figure 10 .
Figure 10.BTCs and fitting results in single fractures with rectangular rough fractures (x = 555 mm).

Table 1 .
Average aperture and relative roughness for different single fractures.

Table 1 .
Average aperture and relative roughness for different single fractures.

Table 2 .
The concentration C-R value equations for different single fractures (unit of C is mg/L).

Table 4 .
Fitting results and goodness of fit for the V-J relationship for flow through single fractures with rectangular roughness.

Table 3 .
Fitting results and goodness of fit for the V-J relationship for flow through single fractures with smooth parallel plates.b(mm)

Table 4 .
Fitting results and goodness of fit for the V-J relationship for flow through single fractures with rectangular roughness.

Table 5 .
Fitting results and goodness of fit for the V-J relationship for flow through single fractures with trapezoidal roughness.

Table 7 .
Fitting parameters of BTC in single fractures with smooth parallel plates.

Table 8 .
Fitting parameters of BTC in single fractures with rectangular rough fractures.

Table 7 .
Fitting parameters of BTC in single fractures with smooth parallel plates.

Table 8 .
Fitting parameters of BTC in single fractures with rectangular rough fractures.

Table 9 .
Fitting results of TPL for single fractures with different ε and b.

Table 10 .
Fitting results of the CTRW TPL model for single fractures with the same b with different roughness structure.