Influence of Shear Displacement on Fluid Flow and Solute Transport in a 3D Rough Fracture

Fractured rocks in the subsurface are ubiquitous, and the dynamics of mass transfer in fractured rocks plays an important role in understanding the problem in engineering geology and environmental geology. In this study, the influence of shear displacement on fluid flow and solute transport in a 3D rough fracture was investigated. A 3D self-affine rough fracture was generated using the modified successive random addition (SRA) technology, and three sheared fractures with different shear displacements were constructed based on the mechanistic model. A direct numerical model based on the Navier-Stokes equation and the advection-diffusion equation was developed to solve the fluid flow and the solute transport. The results showed that shear displacement had a significant influence not only on the fluid flow but also on the solute transport. A global measure of the spatial variability of the flow velocity showed that the heterogeneity became weaker with decreasing shear displacement. All measured BTCs deviated from the Gaussian profile and exhibited the typical anomalous behaviors, such as the long tail and the early arrival. Although the best-fitted results of the advection-dispersion equation (ADE) model and mobile-immobile model (MIM) were generally consistent with those of the BTCs, the MIM was more capable than the ADE model for characterizing the shear-induced anomalous behavior of the BTCs. It was found that the mass exchange process between the immobile and mobile domains was enhanced in the sheared fractures while the fraction of the advection-dominant mobile domain decreased as the shear displacement increased. Furthermore, the deviation of the Taylor dispersion coefficient from the fitted dispersion coefficient by the ADE model and MIM in the sheared fractures was confirmed due to the influence of shear displacement.


Introduction
Fractures are ubiquitous, exist in near-surface bedrock outcrops, and provide dominant pathways for fluid flow in low-permeability rock [1,2]. Accurate characterization of fluid flow and solute transport through fractures plays a crucial role in the contaminant transport, enhancement of oil recovery, nuclear waste disposal, and carbon capture and storage [3][4][5].
Fractures that occur ubiquitously in geological formations are usually simplified as a pair of smooth and parallel plates. In this simplification, the solute transport in the parallel plates is considered to follow Fick's law, in which the dispersion coefficient is spatially and temporally constant. However, for natural fractures, heterogeneity exists at all scales [6,7]. The heterogeneity of fractures is a major source of uncertainty for solute transport through fractures. It has been found that many factors (e.g., roughness [8], contact area [9,10], and topology [11]) significantly affect the solute transport process and cause non-Fickian solute transport. Non-Fickian solute transport is manifested by anomalously early arrival and power-law scaling of late arrival times.
The traditional Fickian transport model, the advectiondispersion equation (ADE), has been widely used to describe solute transport for past decades [12][13][14][15][16]. However, it has been found incapable of capturing the non-Fickian features (i.e., scale-dependent spreading, early arrival, and long tails in BTCs). Non-Fickian solute transport in fracture systems has been widely observed and discussed [17][18][19]. Several models have been developed to overcome the shortcomings of the ADE model and adopted to explain the non-Fickian transport such as the mobile-immobile model (MIM) [20], the continuous-time random walk (CTRW) model [21], and the multirate mass transfer (MRMT) model [22]. Detailed descriptions and comparisons of these transport models can be found in Gao et al. [23].
Previous studies [9,24,25] have used a "single rough fracture" where the lower wall of the fracture is smooth and the upper wall is rough. However, this model has considered the roughness and shown non-Fickian solute transport. It is still far from reality and has little contribution to computation time in numerical simulations. However, it is still a challenge to directly measure fluid flow and solute transport processes within fractures with small apertures and rough surfaces. Benefiting from the development of computer computing power, numerical simulations that can provide pore-scale details are widely adopted to investigate fluid flow in fractures or porous media [10,[26][27][28][29].
The stress fields applied on the fractured rock masses can be altered during earthquakes, tectonic movements, and/or underground excavations, resulting in shear displacements of the fractures. In the past decades, many scholars [30][31][32] have paid attention to the mechanisms of how shear displacements affect the fluid flow process in the inner fracture and its hydraulic conductivity. It has been found that shear displacements can open fractures due to dilation and further increase the conductivity of fractures. Since the solute transport process in fractures highly depends on the fluid flow in pore spaces (i.e., aperture), the solute transport process can be varied in fractures with different shear displacements as the fluid flow paths change significantly [33]. Vilarrasa et al. [34] found that shear displacements greatly increase the flow rate in directions perpendicular to the shearing direction and promote the solute transport in the fracture. Liu et al. [35] calculated the Peclet number during the shearing and found that the dispersion becomes much significant in the increasingly concentrated contact areas and flow channels. However, the studies mentioned before are mostly conducted on two-dimensional (2D) fractures or on single rough fractures. The effects of shear displacement on the solute transport process in threedimensional (3D) rough fractures have not been fully investigated.
The main objective of this study is to investigate the influence of shear displacement on fluid flow and solute transport processes in fractures. For this, three 3D sheared fractures with different shear displacements were constructed based on the modified successive random additions (SRA) and the mechanistic model. A direct numerical model coupled with the Navier-Stokes equation and advection-diffusion equation was implemented to simulate the fluid flow and solute transport in the 3D sheared factures. The characteristics of the flow field and concentration field were analyzed, and the corresponding BTCs were characterized by inversely fitting the ADE model and MIM.

Construction of the 3D Sheared Fracture
There are three steps to construct a 3D sheared fracture with a certain shear displacement. First, a 3D self-affine rough fracture surface was generated. Then, the generated 3D fracture surface was considered to be the lower fracture surface, and the upper fracture surface was translated a vertical distance without any shear displacement. In this step, although both the lower and upper surfaces were rough, the aperture in the constructed 3D fracture was constant. Lastly, the shear process was simulated by moving the upper surface to a certain displacement and keeping the lower surface fixed. Since the shear process was completely determined by the displacement of the upper surface, U s is defined as the shear displacement of the upper surface in the 3D fracture. A more detailed application for constructing the 3D sheared fracture can be found in previous work [9].
Previous studies [8,9,36,37] showed that the aperture in the rough fracture naturally follows a self-affine fractal distribution, implying that the roughness of the rough fracture surface exhibits a repeating pattern at different scales. In this study, a 3D self-affine fracture was generated with the technique of successive random additions [38]. As an efficient and fast algorithm, the fractional Brownian motion was employed based on the technique of successive random additions. It should be noted that as a key characteristic parameter, the Hurst exponent, H, is related to the spatial correlation and the roughness of the fracture surface. Here, a 3D self-affine rough fracture surface with the size of 14 mm length and 10 mm width was generated by selecting H = 0:5, since the H varies between 0.45 and 0.85 for natural rock fractures [39].
Once the 3D self-affine rough fracture surface was generated, the generated rough fracture surface was constructed to a 3D rough fracture with a constant aperture. The shear process was implemented with three different shear displacements (e.g., U s = 2:2, 2.6, and 3.0 mm). It should be mentioned that to ensure the consistency of the constructed 3D sheared fracture, a smaller zone with the size of 10 mm length and 10 mm width in the constructed 3D rough fracture was selected since the shear displacement follows the longitude direction of the generated rough fracture surface. Figure 1 shows the constructed sheared fracture with three different shear displacements. From Figure 1, it can be seen that as the shear displacement increases, not only the void spaces in the sheared fracture change but also the contact area increases. The aperture distribution of the sheared fracture with three different shear displacements is shown in Figure 1. The average aperture is 0.29, 0.48, and 0.67 mm in the sheared fracture with the shear displacements U s = 2:2, 2.6, and 3.0 mm, respectively.

Direct Numerical Model and Inverse Model
3.1. Flow and Solute Transport Equations. In this study, the velocity fields u within the 3D sheared fractures were obtained by solving the continuity and Navier-Stokes equations. The fluid in the flow field was assumed to be incompressible and not affected by gravity.
where u = ½u, v, w is the 3D velocity vector [LT -1 ] with three components, ρ is the density of the fluid [ML -3 ], p is the fluid pressure [ML -1 T -2 ], and μ is the fluid dynamic viscosity [ML -1 T -1 ]. The standard water properties at 293.15 K were used, for example, ρ = 998:2 kg/m 3 and μ = 1:002 × 10 −3 Pa · s. The conservative solute transport within 3D sheared fractures was obtained by solving the advection-diffusion equation.
where c is the solute concentration [ML -3 ], t is time [T], and D m is the molecular diffusion coefficient [L 2 T -1 ]. The velocity vector u is from the flow field based on the solution of Equations (1) and (2). The molecular diffusion coefficient was set to 2:0 × 10 −9 m 2 /s for three sheared fractures. The value of D m relates to the common dissolved organic contaminant (e.g., dichloroethene, toluene, MTBE ethylbenzene, and xylene) in groundwater quality problems [40].
To investigate fluid flow and solute transport processes within the sheared fractures, a pressure gradient was applied at inlet and outlet boundaries (parallel to the y-z plane) to drive fluid flow in the x direction. The side walls (parallel to the x-z plane) and the top and bottom surfaces were set as nonslip boundaries. A dimensionless parameter Pv (pore volume) related to the time is introduced as follows: where Q is the flow rate [L 3 T -1 ] per length and V is the total volume of the void space within the fractures [L 3 ]. In this study, the solute was introduced by a pulse injection method with a flux weight distribution along the inlet boundary. The void space in the fractures was initially saturated with pure water. A total of 1 Pv of a solute, with an initial concentration c 0 of 1 mol/m 3 , was injected at the inlet boundary. Flow and solute transport simulations are performed in the void space by applying pressure gradients that lead to a certain dimensionless Peclet number (Pe) defined as follows: which is the ratio of the characteristic diffusion time ðτ D = ð bÞ 2 /D m Þ to the characteristic advection time where u is the average velocity in the fracture and b is the mean aperture. In this study, solute transport in three sheared fractures was under the same condition Pe = 50. Pe = 50 means that advection dominates, and the velocity is large enough to neglect the effects of molecular diffusion. The dimensionless Reynolds number (Re) is given as follows: where W denotes the width [L] of sheared fractures in the x direction and is equal to 10 mm for all cases. The flow was confined to the laminar regime to avoid the influence of eddies, which occur at high flow velocity on the transport process. Under the condition Pe = 50, the corresponding Re are 0.99, 1.02, and 1.21 for three sheared fractures. The Galerkin finite element software package COMSOL Multiphysics package (COMSOL Inc., Burlington, MA, USA) was used to solve Equations (1)- (3). Since the void spaces in the sheared fracture vary with the shear displacement, three sheared fracture models were meshed into 1.95, 191, and 320 million tetrahedral elements for fractures with the shear displacements U s = 2:2, 2.6, and 3.0 mm, respectively. For all sheared fracture models, the mesh inde-pendence analysis was performed to ensure numerical stability and accuracy (see Figure 2(c)). The steady-state flow rate as a metric of mesh independence was calculated, and the corresponding results indicate that the number of meshed elements was sufficient to provide stable and accurate numerical results.

Inverse Model for Solute Transport
3.2.1. The ADE Model. The classical ADE model, which is based on Fick's law, can be used to quantify the breakthrough caves (BTCs) for solute transport through sheared fractures. The ADE model was obtained from the transformation of the advection-diffusion equation (Equation (3)): replacing the velocity vector u with the average velocity u and changing the molecular diffusion term into the dispersion term. The analytical solution for the pulse injection condition is given by  Lithosphere where C m and C im are defined as the solute concentrations in the mobile and immobile domains, respectively; θ m and θ im are the water contents in the mobile and immobile domains, and the sum of them is equal to one for a saturated fracture; D MIM is the fitted dispersion coefficient; and α is the rate coefficient [T -1 ] associated with the mass transfer between the mobile and immobile domains. A partitioning coefficient β = θ m /ðθ m + θ im Þ is introduced as the mobile water fraction. The STANMOD V2.08 based on the CXTFIT code [41] was used for the inversion of both models. Since the initial estimation of the parameters significantly affects the accuracy of the inverse results. Accurate initial estimation can improve the convergence to a minimum error, especially for the MIM with more inverse parameters. Therefore, the fitting results of the ADE model are used as the initial estimate for the MIM. For both models, the R 2 value and the global error E i were calculated to assess the goodness of fit. The E i is given by where c fit,i is the fitted result and c s,i are the BTCs from numerical simulations, and N is the total amount of the data in BTCs.

Flow Field in Sheared
Factures. On the one hand, the fluid flow as an important background information highly controls the solute transport process in the sheared fracture.
On the other hand, the changed void spaces in the sheared fracture have a significant influence on the characteristics of the flow field. Although the solute transport in all three sheared fractures was under the same Pe, the average flow velocity decreases as the shear displacement increases due to the increasing average aperture as listed in Table 1. The basic parameters of the flow field in three sheared fractures are listed in Table 1. From Figure 3, it can be seen that the shear displacement has a significant influence on the flow field. Generally, the spatial distribution of the flow velocity in the flow field is highly heterogeneous in all sheared fractures. As the shear displacement increases, the originally closed fractured zone (e.g., the aperture is zero) becomes dilated. Comparing  Figures 3(a) and 3(c), this dilated void space not only provides more lager void spaces but also decreases the permeability of the fracture. For all sheared fractures, the relatively high local flow velocity is captured at both the small-aperture zone and the intersecting fractured zone, while the relatively low local flow velocity is observed at the large-aperture zone. The streamlines become smoother in the sheared fracture with the shear displacement U s = 3:0 mm than with the shear displacement U s = 2:2 mm. This indicates that the tortuosity of the flow path in the sheared fracture is significantly influenced by the shear displacement. As a result, as both shear displacements and dilated void spaces increased, the preferential flow paths disappeared.
Since the preferential flow path development in the flow field is the result of the complex fractured space, the presence of the preferential flow is related to the heterogeneity of the flow field. To quantify the shear-induced flow field where the preferential flow occurs, the coefficient of flow velocity variation, a global measure of the velocity spatial variability [27,42], was calculated for all sheared fractures: where σ u is the standard variance of the flow velocity, u av is the average velocity magnitude, and V is the volume of the pore space occupied by fluid in the sheared fracture. In Table 1, it is found that the calculated CV u decreases as the shear displacement increases. This indicates that the velocity spatial variability is enhanced by the decreasing shear displacement. These results of the calculated CV u are also consistent with the characteristics of the spatial distribution of the flow velocity in the flow field.

Concentration Field and BTCs.
In this study, to better capture the characteristics of solute transport in the sheared fracture, the solute was injected by the step injection boundary condition. The solute transport in all sheared fractures was carried out under the asymptotic transport regime in which the Pe is less than the proposed critical Pe a . This could ensure the characteristics of non-Fickian transport as a  6 Lithosphere result of the shear displacement influence rather than the fracture length scale. Figure 4 shows the spatial-temporal evolution of solute transport in the 3D sheared fracture with different shear displacements. From Figure 4, it can be seen that the shear displacement has a significant influence on solute transport. The evolution of the solute in the sheared fracture is sensitive to the local aperture distribution. In the sheared fracture with the shear displacement U s = 2:2 mm, the solute transports faster in the preferential path along the flow direction than other largeaperture zones. Since the flow velocity is relatively high in the preferential path across the entire sheared fracture, the solute could break through early the entire sheared fracture. As the shear displacement increases, not only the area of the closed fractured zone decreases but also the flow velocity in the preferential path decreases. In Figures 4(f)-4(j), it can be seen that the solute transports in the vicinity of the closed fractured zone. It shows clearly that there are two preferential paths in the sheared fracture with the shear displacement U s = 2:6 mm. Due to the heterogeneous aperture distribution, the solute is transported slower in the newly formed preferential path. This nonuniform transport of the solute could be attributed to the spatial variability of flow velocity in the flow field. It should be noted that although the Pe was the same and set to Pe = 50 for all sheared fractures, the local flow velocity was varied depending on the local aperture. This indicates that the solute transport was advection-dominant at the scale of the entire sheared fracture but was diffusion-dominant in the small-aperture zone where the flow velocity is relatively low.
Once the solute was transported into the small-aperture zone, the slow flow and diffusion dominate the solute transport. Consequently, the transport of the solute could be significantly delayed. In fact, compared to the preferential path that is advection-dominant, the small-aperture zone can be considered an "immobile" region and provide strong resistance against the spreading and mixing of the solute. In addition, since the nonslip boundary was implemented in 7 Lithosphere the direct numerical model, the distribution of the local flow velocity is nonuniform. The local flow velocity at the crosssection of the sheared fracture shows a parabolic distribution. This leads to slower solute transport near the fracture surface than in the middle of the sheared fracture. Figure 5 shows the BTCs for three sheared fractures. The flux-averaged BTCs were determined to investigate the influence of shear displacement on the BTCs in the sheared fractures. The flux-averaged concentration was normalized as From Figure 5, it can be seen that the shear displacement has a significant influence on the characteristics of BTCs. All determined BTCs show the deviation from the Gaussian profile and the typical anomalous behaviors, such as the long tail and the early arrival. Although the average velocity in all the sheared fractures was the same, the breakthrough time for the BTCs was influenced by the shear displacement. As the shear displacement increases, the breakthrough time for the BTCs increases, associated with the increasing residence time of the solute. However, it was found that the breakthrough time is the same in the sheared fractures with the shear displacements U s = 2:6 and 3.0 mm. This is because the spatial variability of the shear-induced flow in the sheared fractures with the shear displacements U s = 2:6 and 3.0 mm is the same (see Table 1), which is in accordance with the results of the flow field in Section 4.1. Since the flow field is highly heterogeneous in all the sheared fractures, the long tail of BTCs was observed. From Figure 5, it can be seen that the long tail of BTCs was stronger in the sheared fracture with the shear displacements U s = 2:6 and 3.0 mm than with the shear displacement U s = 2:2 mm, indicating that the shear-induced flow increases solute residence time. Another interesting finding is the influence of shear displacement on the peak concentration of BTCs. As shown in Figure 5, the peak concentration of BTCs is lower for the sheared fracture with the shear displacements U s = 2:6 and 3.0 mm than with the shear displacement U s = 2:2 mm, indicating that the shear displacement increases as the peak concentration of BTCs decreases. From the perspective of the solute mixing process in the sense that the greater magnitude of mixing leads to a lower peak concentration, it can be concluded that the magnitude of solute mixing is enhanced by the shear displacement. It should be noted that all determined BTCs were analyzed under the same Pe. Therefore, the difference of characteristics of BTCs for all cases is attributed to the influence of the shear displacement.

Application of the ADE Model and MIM.
To further investigate the shear displacement influence on the solute transport behavior, both the ADE model and the MIM were applied inversely to fit the determined BTCs. From Figure 6, it can be seen that the best-fitted lines of both the ADE model and the MIM agreed generally with the results of the determined BTCs from the direct numerical simulations for all sheared fractures. The determination of errors in both the ADE model and the MIM is more than 0.99 for all sheared fractures as listed in Table 2. However, the ADE model poorly fits the typical anomalous characteristics of BTCs (e.g., the long tail and the early arrival). This also can be seen from the global error in both the ADE model and the MIM where the global error is smaller in the MIM than the ADE model for all sheared fractures. Based on the results from Tables 2 and 3 and Figure 6, although the results from both the ADE model and the MIM tend to yield an increasing global error as the shear displacement increases, the best-fitted results of BTCs are better from the MIM than from the ADE model, indicating that the MIM is more capable than the ADE model for characterizing the shear-induced anomalous behaviors of BTCs.
In Tables 2 and 3, both the dispersion and the average flow velocity in three sheared fractures are estimated inversely from the ADE model and MIM. Compared to the results in Table 1, it is found that the best-fitted average flow velocities are lower than the average flow velocities calculated from the direct numerical model in the sheared fractures with the shear displacements U s = 2:2 and 2.6 mm for both the ADE model and the MIM. This is because, on the one hand, the transport velocity in the small-aperture zone is slower than the average flow velocity; on the other hand, the presence of the closed fractured zone leads to the strong preferential flow path where the flow velocity is larger than the average flow velocity. Moreover, the deviation of the fitted velocity from the average flow velocity increases as listed in Tables 2 and 3 for both the ADE model and the MIM. This also indicates that the closed fractured zone has a significant influence on the best-fitted average flow velocities from the ADE model and MIM. It should be noted that the deviation of the fitted velocity from the average flow velocity from the MIM is larger than that from the ADE model. This is attributed to the fact that the MIM characterized the velocity of the solute transport only in the mobile domain.
Unlike the ADE model, the estimated parameters α and β from the MIM provide insight into the mechanism of the influence of shear displacement on solute transport. In   Table 3, it is found that the α increases as the shear displacement increases while the β decreases as the shear displacement increases. The increasing α indicates that the mass exchange process between the immobile and mobile domains is enhanced. The decreasing β means that the fraction of the diffusion-dominant immobile domain increases. Although the area of the closed fractured zone decreases as shear displacement increases, the fraction of the mobile domain increases due to the decreasing spatial variability of velocity (see Table 1). It also can be seen from Figure 4 that as the area of the closed fractured zone decreases, the patterns of the solute transport are more uniform in the sheared fracture with the shear displacement U s = 3:0 mm than with the shear displacements U s = 2:2 and 2.6 mm.
In addition, the dispersion in a fracture can be directly estimated by the Taylor analytical solution. The seminal Taylor's work [43] demonstrates that for solute transport through parallel plates where the fluid follows the stratified flow, the effective solute transport eventually becomes one dimensional and the shear flow-induced dispersion can be measured by the effective longitudinal dispersion (i.e., Taylor dispersion) and average flow velocity. The Taylor dispersion coefficient for an ideal paralleled fracture is given as follows: Depending on Equation (12), the Taylor dispersion coefficient can be directly estimated by the Q and D m without conducting solute transport simulations. The Peclet number corresponding to the critical advection time, Pe a , is determined depending on our previous work [8]. The Pe a is equal to 138.48, 82.59, and 60.02 for the sheared fractures with the shear displacements U s = 2:2, 2.6, and 3.0 mm, respectively. It should be mentioned that all Pe in three sheared fractures are equal to 50 and smaller than the determined Pe a , which indicates that all of the simulated solute transport in this study is asymptotic and the anomalous solute transport behavior could  9 Lithosphere be considered a result of the influence of shear displacement. The Taylor dispersion coefficient in three sheared fractures was calculated depending on Equation (12). The D Taylor was 4:68 × 10 −8 , 3:27 × 10 −8 , and 2:58 × 10 −8 m 2 /s for the sheared fractures with the shear displacements U s = 2:2, 2.6, and 3.0 mm, respectively. It is found that all calculated D Taylor in three sheared fractures are smaller than the D ADE and D MIM . This deviation of the D Taylor from the D ADE and D MIM is due to the fact that the Taylor analytical solution for D Taylor is under the assumption of an ideal paralleled fracture, while the sheared fractures in this study are rough fractures. The magnitude of this deviation indicates the influence of the shear displacement on the dispersion of the solute.

Summary and Conclusions
In this study, we focused on the characteristics of fluid flow and solute transport through the sheared fracture and the influence of shear displacement on fluid flow and solute transport in the 3D sheared fractures was numerically investigated.
It was found that as both shear displacements and dilated void spaces increase, the spatial variability of flow velocity decreases. The presence of the preferential flow is found to be related to the heterogeneity of the flow field. Non-Fickian behavior (i.e., early arrival and long tails) was observed in all BTCs in the sheared fractures. It was found that as the shear displacement increased, the breakthrough time for the BTCs increased. Inversely applying the ADE model and MIM to characterize the BTCs showed that the MIM is more capable than the ADE model for characterizing the shear-induced anomalous behaviors of BTCs. It was also confirmed that the mass exchange process between the immobile and mobile domains was enhanced as the shear displacement increased.
Although our work highlighted the influence of shear displacement on fluid flow and solute transport in the sheared fractures, future studies would focus on the different directions of shear displacement that influences the fluid flow and solute transport. Furthermore, in this study, it was assumed that the rock matrix was impermeable, and thus, the fracture surface is a nonslip boundary. A lowpermeability rock matrix should be considered because this complex fracture system is of key importance to hydraulic fracturing, back-diffusion processes, oil recovery efficiency, and groundwater remediation.

Data Availability
The data used to support the findings of this study are included within the paper.

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.  10 Lithosphere