Lingcho Zou, Vldimir Cvetkovic

a Division of Resources,Energy and Infrastructure,Department of Sustainable Development,Environmental Science and Engineering,Royal Institute of Technology,Stockholm,10044, Sweden

b Department of Physical Geography, Stockholm University, Stockholm,10691, Sweden

ABSTRACT The impact of normal stress-induced closure on fluid flow and solute transport in a single rock fracture is demonstrated in this study.The fracture is created from a measured surface of a granite rock sample.The Bandis model is used to calculate the fracture closure due to normal stress,and the fluid flow is simulated by solving the Reynold equation. The Lagrangian particle tracking method is applied to modeling the advective transport in the fracture.The results show that the normal stress significantly affects fluid flow and solute transport in rock fractures.It causes fracture closure and creates asperity contact areas,which significantly reduces the effective hydraulic aperture and enhances flow channeling. Consequently, the reduced aperture and enhanced channeling affect travel time distributions. In particular, the enhanced channeling results in enhanced first arriving and tailing behaviors for solute transport. The fracture normal stiffness correlates linearly with the 5th and 95th percentiles of the normalized travel time. The finding from this study may help to better understand the stress-dependent solute transport processes in natural rock fractures.

Keywords:Normal stress Fluid flow Solute transport Stiffness Particle tracking

1. Introduction

Modeling of fluid flow and solute transport processes in fractured rocks is an important topic for many rock engineering projects, such as mining, groundwater management and geological disposal of nuclear wastes (Neuman, 2005). In particular, solute transport in rock fractures, quantified as the residence time distribution, is critical in safety assessment of geological disposal of nuclear wastes(e.g.Cvetkovic et al.,1999;Bodin et al.,2003;Tsang et al., 2015). A fundamental understanding of solute transport in single/discrete fractures is essential for modeling of flow and transport in fracture networks at field application scales(e.g.Zhao et al., 2011; Cvetkovic and Frampton, 2012). Numerous experimental and numerical studies have shown that the fluid flow and solute transport processes in rock fractures are largely affected by mechanical behaviors (e.g. Gangi, 1978; Kranz et al., 1979;Witherspoon et al., 1980; Tsang and Witherspoon, 1981; Barton et al., 1985; Zimmerman and Bodvarsson, 1996; Yeo et al., 1998;Rutqvist and Stephansson, 2003; Selvadurai, 2015; Kang et al.,2016; Zou et al., 2017a; Vogler et al., 2018). For instance, Li et al.(2008), Zhou et al. (2015) and Chen et al. (2019) conducted a series of laboratory tests on fluid flow in rough fractures coupling with normal or shear stress, and they found that the effective transmissivity or hydraulic apertures were linear or nonlinear functions of applied stresses. Selvadurai (2015) proposed that the stress relief of a fracture could result in a significant permeability increase by radial flow hydraulic pulse testing of a fracture in a granite cylinder. Vogler et al. (2018) presented experimental test and numerical simulation results of hydro-mechanical behaviors of rough fractures exposed to high-pressure fluid injection,and found a nonlinear increase in fluid injection pressure with loading. The discovered stress-dependent flow would consequently affect the solute transport processes (e.g. Zhao et al., 2013). This paper is an attempt to improve the basic understanding of solute transport in a single rough-walled rock fracture affected by normal stress.

The normal stress affects the flow and transport in rock fracture because it causes fracture asperities to deform and close,which can largely alter the geometrical structures for fluid flow and solute transport (e.g. Pyrak-Nolte and Morris, 2000; Kang et al., 2016).Over the past decades, many studies have sought to quantify the constitutive equation for the fracture closure due to the normal stress. Several classical empirical relationships between normal stress and normal displacement have been established based on a large number of experimental data. For example, Goodman (1976)developed an empirical hyperbolic function to relate the fracture closure with applied normal stress.Bandis et al.(1983)introduced the initial normal stiffness of the fracture interface and the maximum closure to uniquely define the hyperbolic stress-closure function for rock fractures.Considering the elastic contact of rough rock surfaces during closure,a class of analytical models applied or extended the Greenwood—Williamson (GW) model based on the Hertzian contact theory for rock fractures (Greenwood and Williamson,1966; Brown and Scholz,1985,1986; Xia et al., 2003;Tang et al., 2016). In addition, many numerical models have been developed to simulate the contact of rough fractures under normal stresses, commonly using direct finite element simulations (e.g.Lavrov,2017),and boundary element modeling based on half-space approaches (e.g. Hopkins, 1991; Li et al., 2015; Kang et al., 2016;Kling et al., 2018). These mechanical models provided convenient approaches for quantifying fracture closure. In particular, the classical empirical models provided simple constitutive equations to estimate the fracture closure directly and have been widely used in rock engineering.

Based on the quantified constitutive equations or numerical simulations for fracture closure under normal stress, many studies investigated the correlations between the fracture specific stiffness (defined as the ratio between the increment of stress and the increment of displacement caused by the deformation of the fracture) and effective transmissivity. For instance, Tsang and Witherspoon (1981), Hopkins (1991),Pyrak-Nolte and Morris (2000) and Pyrak-Nolte and Nolte(2016) studied the impact of normal stress on fracture transmissivity. An intrinsic correlation between fracture specific stiffness and fracture effective transmissivity has been presented by Pyrak-Nolte and Morris (2000) and Pyrak-Nolte and Nolte (2016). Using Hopkins method, Wang and Cardenas(2016) studied the impact of surface roughness on the correlation between fracture specific stiffness and fracture effective transmissivity. These studies generally revealed that the fracture transmissivity is intrinsically dependent on the fracture normal stiffness. However, most of these studies only focused on the hydro-mechanical processes, whereas the correlation between fracture stiffness and solute transport properties remains a comparatively less studied topic.

Solute transport in rock fractures involve complex processes,including advection, diffusion, decay, chemical reaction and mass exchanges with rock matrix(Bodin et al.,2003).Numerous studies have shown that the transport behaviors in rock fractures are largely controlled by the flow properties due to the inherent physical process of advection (e.g. Cvetkovic et al., 1999; Bodin et al., 2003; Zou et al., 2016). The stress-dependent permeability of rock fractures would consequently affect the transport processes.Moreover, solute transport processes, e.g. matrix diffusion, is also dependent on geometrical structures that can be altered by the normal stress(Cvetkovic et al.,1999; Neretnieks, 2014). Therefore,solute transport processes in rock fractures are most likely also intrinsically dependent on the normal stress. A few studies attempted to quantify the impact of normal stress on the solute transport in rock fractures. For example, Kang et al. (2016) presented a generic study on the emergence of anomalous transport in randomly generated rough-walled fractures under increasing normal stresses,showing that the normal stress results in anomalous transport in rock fractures. However, due to the complicated mechanical and transport processes combined with complex structure, the impact of normal stress on the solute transport in natural rough-walled rock fractures is still not fully understood.

The objective of this work is to quantify the impact of fracture closure induced by normal stress on fluid flow and solute transport in a rough-walled fracture through numerical modeling. The primary questions we attempt to address in this work are:(1)what are the mechanisms for normal stress affecting solute transports in natural rough-walled rock fractures? and (2) what is the possible correlation between the fracture specific stiffness and advective solute transport? To answer these two questions, a rough-walled fracture model is created based on a digital surface measured from a granite rock sample. At this stage, we only use the most convenient mechanical and fluid flow models and only consider advective transport.The Bandis model is employed to calculate the fracture closure and stiffness subjected to normal stress. The fluid flow is simulated by solving the Reynolds equation using finite element method (FEM). The advective transport is simulated through a Lagrangian particle tracking approach. The results from this study improve our understanding of the impact of normal stress on solute transport in a natural crystalline fracture on the laboratory scale.

2. Materials and methods

2.1. Rough fracture model

In order to represent realistic feature roughness, a digital fracture surface scanned from a granite rock sample is used as the parent surface to build a rough-walled fracture geometry model for analyzing fracture closure,fluid flow and solute transport.The size of original rock specimen is 0.1 m in width and 0.2 m in length.This digital surface is obtained with an interval of 0.2 mm in bothx-andy-direction with an accuracy of ±20 μm and a resolution of 10 μm(Koyama et al.,2008).We use an ‘a(chǎn)rtificial shear’step inx-direction by 1 mm and then a ‘lift up’ step by 0.8 mm, to create a more realistic rough-walled fracture with variable apertures (Zou et al.,2017a), where the original surface is the lower surface and the sheared and lifted surface is the upper surface. The digital surface and the approach of creating the rough fracture model are illustrated in Fig.1.

The ‘a(chǎn)rtificial shear’and ‘lift up’steps create a numerical fracture model with variable apertures, which performs like a simplified shear test conducted in the laboratory,under small or zero normal stress. Although this ‘a(chǎn)rtificial shear’ is largely simplified in comparison to actual laboratory tests, it provides a fracture geometry model with realistic rough surfaces. A 0.05 m × 0.05 m section selected from the created fracture will be the rough-walled fracture geometry model for analyzing the fracture closure, fluid flow and solute transport (Fig. 1). The selected section avoids the effect of sample edges as well as the extreme asperities(i.e.the highest and the lowest surface heights near the edges).

The generated fracture model with variable apertures and the probability density function (PDF) of its fracture aperture are presented in Fig. 2. The local apertures are discretely distributed in space, with different correlation lengths inx- andy-direction(Fig. 2a) because of the artificial shear inx-direction. The mean aperture(0.804 mm)roughly equals the lift up distance,i.e.0.8 mm and the standard deviation is around 0.156 mm.The distribution of the aperture is well approximated by the normal distribution (the red curve in Fig.2b):this feature of the aperture has been reported in previous study for artificial shearing of self-affine surfaces by Mallikamas and Rajaram (2005).

Fig.1. The original digital surface and illustration for creating the rough fracture model using an ‘a(chǎn)rtificial shear’ approach.

2.2. Normal stress-induced fracture closure

Many analytical or empirical closure models for rock fractures under normal stresses have been proposed based on theoretical analyses or experimental data. One of the most widely used empirical closure model is the Bandis model, which is used in the present study to model the closure. Bandis model is a hyperbolic formula that describes the relationship between the effective normal stress and the closure displacement, expressed as (Bandis et al.,1983):

where σnis the effective normal stress(Pa);kn0is the initial normal stiffness coefficient(Pa/m);umaxis the maximum closure(m);and Δuis the closure displacement (m) at a specific normal stress,which can be rewritten as (Bandis et al.,1983):

Accordingly, the stiffnesskncan be expressed as (Bandis et al.,1983):

The parameters of the maximum closureumaxand initial normal stiffness coefficientkn0for a specific fracture are normally determined by compression experiment. In this work, we assume that the maximum closureumaxis equal to the initial mean aperture valuee0, i.e.umax=e0, for a generic study. The underlying assumption is that the hydraulic aperture is equal to the mechanical aperture and the fracture can be closed completely.Therefore,for a single fracture, with known initial mean aperture valuee0and initial normal stiffness coefficientkn0, the apertureeavailable for flow at a normal stress can be obtained by

2.3. Fluid flow

The Reynolds equation has been widely applied in modeling of flow and transport in rock fractures,due to its simplicity and ability to account for surface roughness by assuming that the cubic law is applicable locally (e.g. Brown, 1987). The Reynolds equation is expressed as

wherePis the hydraulic pressure(Pa);μ is the dynamic viscosity of water(Pa s);andTis the transmissivity obtained by the cubic law,written as

wheree(x，y)is the local aperture.Note that the Reynolds equation is valid only for laminar flow conditions when the Reynolds number is less than 1.The Reynolds number is defined asRe=ρQ/(Wμ),where ρ is the water density,Qis the specific flow rate,andWis the width of fracture.

Fig. 2. Aperture distributions of the fracture model: (a) Spatial distribution; and (b)Probability density function, where the red curve is a normal probability density function. Std denotes the standard deviation.

A Galerkin FEM code is developed to solve the Reynolds equation and to obtain the flow fields.The structured mesh with bilinear quadrilateral elements is adopted for finite element simulations.The structured mesh with the size of 0.2 mm×0.2 mm is consistent with the grid of the scanned surfaces, so that it can represent the full geometrical features of the surface roughness from the measured topography data. In addition, the bilinear quadrilateral elements instead of the linear triangular elements are used to avoid constant gradient at each triangle for the linear triangular elements. The Gaussian quadrature method using four Gaussian integration points are adopted in the FEM code.The local aperturee(x，y) at each Gaussian integration point used to calculate the transmissivityT(see Eq.(6))is interpolated from the four nodes of an element using its shape functions. By solving the Reynolds equation (Eq. (5)), we can obtain the pressure field. The flow velocity field can be calculated by

whereuis the velocity vector (m/s).

2.4. Advective transport

In this study, we only consider the advective transport because of its dominant impact on contaminant transport in natural (heterogeneous)rock fractures(Cvetkovic et al.,1999). The Lagrangian particle tracking approach is used to simulate the advective transport in the fracture model. Compared to the Eulerian gridbased methods, such as finite difference, finite element and finite volume methods, the Lagrangian particle tracking method avoids errors caused by numerical dispersion, so that it has unique advantages to model the advective transport process (e.g. Cvetkovic et al.,1999; Salamon et al., 2006).

In this Lagrangian particle tracking approach,the flow trajectory is described by the Lagrangian position vector as a function of time:X(t) = [x(t)，y(t)].In this study,we adoptx-direction as the main flow direction. The travel/residence time τ for a particle from the release point (x= 0) to the control plane (x=L) is defined by

whereuxis the velocity inx-direction. In practice, the residence time τ of each particle is calculated by its discrete form:

where the residence time τ is the sum of travel time at each stepiwhen a particle is advected from the release point (x= 0) to the control plane (x=L). The classical fourth-order Runge-Kutta method is adopted for the particle tracking. The details of this particle tracking method are available in Zou(2016).

3. Simulation results of fracture closure,fluid flow and solute transport

3.1. Fracture closure and specific stiffness

According to the single fracture model created from the digital surface, the initial mean aperture is 0.8 mm. A typical value of the initial normal stiffness for granite,i.e.kn0= 25 GPa/m,is adopted to parameterize the mechanical model (e.g. Bandis et al., 1983).Only one value for the initial normal stiffness is used, since the impact and sensitivity of the initial normal stiffness are obvious according to the Bandis model that has been adopted in this study for illustration.In order to quantify the impact of normal stress on fluid flow and solute transport, eight cases with varying normal stresses from 0 MPa to 35 MPa will be considered in this study.The model parameters and simulated cases are summarized in Table 1.

The fracture closure and stiffness under the applied normal stresses can be obtained directly from the Bandis model,i.e.by Eqs.(2)and(3).The calculated fracture closure and stiffness as functions of the normal stress are presented in Fig.3,showing that both the fracture closure and stiffness increase with the increasing normal stress.However,the fracture closure increases with a reducing rate of change, whereas the normal stiffness increases with an increasing rate of change. The fracture closure is used to calculate the aperture distributions for each case of normal stress, and the stiffness is used to explore its correlations with contact areas and transport properties.

Table 1Summary of model parameters and simulated cases.

Fig. 3. Fracture closure and specific stiffness as functions of normal stress.

Fig.4presentsthePDFsof thefractureaperturefordifferentnormal stresses, i.e. σn= 0 MPa， 20 MPa and 35 MPa. It shows that the mean aperture decreaseswith theincreasingnormalstress.Compared with the case σn= 0 MPa， when σn= 20 MPa and 35 MPa， the probability density curves are truncated on the left-hand side due to the contact areas caused by the fracture closure. These increasing contact areas also slightly reduce the standard derivation of the aperture with the increasing normal stresses.

The asperity contact areas directly influence the normal stiffness,because the applied normal stresses are held by the asperity contact areas. The asperity contact area is defined as the overlapping zones of the upper and lower fracture surfaces after closure. We use the mean aperture (i.e. subtract height of lower surface from the height of upper surface) of each element to determine the contact areas. Once the mean aperture of an element is below zero, this element is considered as contacted.The relative contact areas is the ratio between the total contact areas and the entire vertical projected fracture area(i.e.0.05 m×0.05 m). Fig. 5 presents the correlation between the stiffness and the relative contact areas, showing that the stiffness increases with increasing relative contact areas.Note that the contact area is also affected by the geometry condition of surface asperities,especially for the initial contacting stage under relatively small normal stress;there are negligible contact areas when the normal stress is less than 5 MPa for instance. This part of the results is affected by the initial mean aperture. When the relative contact area is larger than 1%, the stiffness is essentially linear with the relative contact area.

Fig. 4. Evolution of aperture distributions with the increasing normal stress.

Overall, increasing normal stress causes nonlinear increases of fracture closure and stiffness. The increasing fracture closure gradually reduces the mean aperture and increases the contact area.These changes of aperture structure and stiffness are the main reason accounting for changes in fluid flow and solute transport.

3.2. Flow field and effective transmissivity

To obtain the flow fields by solving the Reynolds equation, a constant pressure gradient is applied inx-direction,where the lefthand side is the inlet boundary and the right-hand side is the outlet boundary. We only consider laminar flow in the fracture, and a relatively small pressure gradient, i.e. 200 Pa/m, is applied to ensuring that the Reynolds number is less than 1 for all simulation cases. To calculate flow streamlines for visualization of the flow field and analyzing advective transport,5000 particles are injected at uniform spacing along the inlet boundary (see Fig. 1) for all simulation cases. The resident injection method rather than the flux-weighted method is used to capture the overall features of flow and transport in the whole fracture.

Fig. 6 presents the simulated velocity fields and streamlines under three exemplified cases with normal stresses of σn= 0 MPa, 20 MPa and 35 MPa. Initially, when σn= 0 MPa, the flow shows slight channeling due to variable apertures caused by heterogeneous transmissivity. In contrast, for the cases with normal stresses σn= 20—35 MPa, there are increasing areas of irregular contacts caused by the fracture closure. These contact areas block the flow and push the flow streamlines into a few channels around the contact areas. When σnincreases from 20 MPa to 35 MPa, the streamlines indicate that the flow becomes more channelized due to increased contact area.

Fig. 5. Correlation between stiffness and relative contact areas.

Fig. 6. Velocity fields and streamlines for different normal stresses σn: (a, c, e) Velocity fields, and (b, d, f) Streamlines for σn = 0 MPa, 20 MPa and 35 MPa, respectively.

In order to quantify the dependency of flow properties on the normal stress, the effective hydraulic aperture is back-calculated from the simulated flow rates using the cubic law(see Eq.(6)).The relationship between the effective hydraulic aperture and the applied normal stress is presented in Fig.7,showing that the effective hydraulic aperture decreases with increase of applied normal stress.This feature of stress-dependent hydraulic aperture has been extensively reported in previous studies (e.g. Swan, 1983;Zimmerman et al.,1992; Renshaw,1995; Jeong, 2017). The contact area as a function of normal stress is also presented in Fig.7,showing that it largely increases when σn>10 MPa. The main reason for decreasing effective hydraulic aperture is that the mean aperture is reduced due to the increasing fracture closure. Another reason is that the increasing contact area enhances channeling flow and increases the flow tortuosity that results in a larger pressure drop.

Fig. 7. Hydraulic aperture and relative contact area as functions of the normal stress.

3.3. Travel time distributions

The travel time for each particle is recorded as part of the particle tracking algorithm. The distribution of particle travel time represents the overall transport behavior.Note that the mean travel time is different for various cases of normal stress because of the effective hydraulic aperture and the mean velocity decreases with the increasing normal stress. This feature can be easily extracted and estimated by the cubic law once the effective hydraulic aperture is obtained. In practice, it is more important to quantify the deviation between a realistic travel time and the mean travel time estimated from the simplified homogeneous model.Therefore,the travel time is normalized by the mean travel time obtained from the cubic law based on an effective hydraulic aperture for the corresponding case of normal stress.

Fig. 8 presents the cumulative distribution function (CDF) and complementary CDF (CCDF) of the normalized travel time for the eight cases of normal stress.The CDF is used to show the first arrival features and the CCDF is used to show the tailing behavior.Generally, the CDFs and CCDFs show multi-rate transport features even for the case without normal stress because of the channeling(see Fig.6).It is expected that this multi-rate features are gradually enhanced with increasing normal stress because of enhanced channeling. More importantly, the large difference between the first arrivals and tailings for different cases of normal stress indicates that the normal stress significantly influences the breakthrough and tailing.Specifically,the solute arrives relatively earlier than the corresponding mean travel time when the normal stress is high. Meanwhile, the solute has relatively longer travel time for higher normal stress cases.

To quantify the influence of normal stress on the relative first arrival, median and tailings, the 5th, 50th and 95th percentiles of the normalized travel time as a function of normal stress are presented in Fig.9.The median travel time for all cases of normal stress is around 5%smaller than the corresponding mean travel time,and it slightly reduces when the normal stress increases from 0 MPa to 30 MPa. The 5th percentile of travel time is around 20% smaller than the mean travel time when the normal stress is 0 MPa, and then it gradually reduces to 45% smaller when the normal stress increases from 0 MPa to 35 MPa.In contrast,the 95th percentile of the travel time is 28% higher than the mean travel time when the normal stress is 0 MPa, and then it gradually increases to 93%higher when the normal stress is 35 MPa. Comparing the impacts between the first arrival represented by the 5th percentile of travel time and the tailings represented by the 95th percentile, it shows that the normal stress has a greater influence on asymptotic travel time.The reason is that the increasing contact areas caused by the elevated normal stresses induce more impact on the tortuosity of flow streamlines than the enhanced channeling flow (see Fig. 6),resulting in the higher magnitude of influence on asymptotic travel time.

Fig.8. (a)Cumulative distribution function(CDF)and(b) Complementary CDF(CCDF)of travel time for different normal stress cases.

Fig. 9. The 5th, 50th and 95th percentiles of the normalized travel time for different normal stresses.

3.4. Correlation between fracture specific stiffness and travel time

Here we explore the correlation between fracture normal stiffness and typical travel time percentiles,i.e.the 5th,50th and 95th,of the normalized travel times, respectively, as shown in Fig. 10.Since the stiffness is a function of the normal stress,the correlation between fracture stiffness and travel time percentiles exhibits a similar pattern as the relationship between travel time percentiles and normal stress.The 5th percentiles of the normalized travel time decrease linearly with increasing stiffness, whereas the 95th percentiles of the normalized travel time increases linearly with increasing stiffness. This indicates that potentially strong linear correlation can be established between the fracture normal stiffness and key arrival time percentiles for advective solute transport in heterogeneous rock fractures.In practice,such correlation can be applied to predicting the transport properties from the stiffness information, and vice versa.

4. Discussions

We demonstrated the impact of normal stress on advective solute transport in a rough-walled rock fracture created from a measured rock surface. The simulation of fluid flow and particle advection confirmed that the normal stress could have a significant impact on solute transport. Both the relative first arrivals and tailings were enhanced by increasing normal stress. The main mechanism is that the normal stress causes fracture closure and creates asperity contact areas.These contact areas result in channeling flow that consequently enhances the dispersion processes in the fracture. The first arrival and tailing are important features in risk assessment or safety evaluation for contaminant transport in fractured rocks (e.g. Frampton and Cvetkovic, 2011). The revealed stress-dependent transport behavior in rock fractures confirmed that under certain conditions of stress change, it might be important to consider normal stress impact on fluid flow and solute transport. Applying simplified homogeneous models to simulate fluid flow and solute transport in natural fractures may introduce higher degree of errors or uncertainties for the cases with relatively high normal stress.

Fig.10. Correlation between fracture stiffness and normalized travel time percentiles.

The correlation between the normal stiffness and travel time percentiles had been explored in this study for the first time. The results indicated a potentially linear correlation between the stiffness and normalized first arrival and asymptotic arrival time in the fracture. It is similar to the correlation between the fracture stiffness and flow properties,i.e.flow rates or permeability,addressed in previous studies (e.g. Hopkins, 1991; Pyrak-Nolte and Morris,2000;Pyrak-Nolte and Nolte,2016;Wang and Cardenas,2016).The correlation between stiffness and transport properties provided a potentially useful relationship to link hydro-mechanical and chemical transport properties,e.g.to predict transport parameters from the measured stiffness or flow properties. Similarly, once reliable correlations between stiffness and transport properties were determined, the hydro-mechanical properties of rock fractures,e.g.compressibility and permeability,could also be estimated from measured transport information, e.g. the measured breakthrough curves by tracer tests. However, the present study only explored such potential correlations using relatively simple mechanical and flow models. More comprehensive modeling and validation studies are needed to obtain more reliable correlations that can be used in applications. Nevertheless, the findings from this study demonstrated the importance and possibilities of such correlations between hydro-mechanical and chemical transport properties.

Several remaining issues are summarized for future studies as follows:

(1)Dependency of geometry structures.The complex surface roughness that yields heterogeneous structures of rock fractures is the fundamental reason for the uncertainty in the correlations between hydro-mechanical and chemical transport properties.Although a measured surface was used to create the fracture model to represent realistic surface roughness features of rock fractures in the present study,the dependency of geometry structures on the correlation between hydro-mechanical and chemical transport properties has not been considered. This is limited by the relatively simple mechanical model (Eq. (2)), and the simplified flow model,i.e.the Reynolds equation(Eq.(5))used in our study.Such models that use simplified roughness parameters, e.g.joint roughness coefficient (JRC), may have limitations in representing the impact of geometry structures(Bandis et al.,1983; Barton et al., 1985). Therefore, more sophisticated mechanical models,such as direct finite element(e.g.Lavrov,2017) or boundary element simulations (e.g. Li et al., 2015;Kang et al., 2016; Kling et al., 2018), and rigorous flow equations, such as Navier—Stokes equations (Zou et al.,2017a, b) are needed to study the dependency of geometrical structures on the correlation between hydromechanical and chemical transport properties.

(2)Dependency of stress history.Natural rock fractures or faults have often undergone a complex stress history of loading and unloading process. The cyclic loading process will cause irreversible changes of the fracture geometric structures,which has been observed in laboratory tests of single rock fracture (e.g. Bandis et al., 1983; Brown and Scholz, 1985,1986; Boulon et al.,1993; Selvadurai, 2015). These irreversible changes will consequently affect the hydro-mechanical and chemical transport behaviors of the fractures or faults(e.g.Sathar et al.,2012).In this work,we did not consider the cyclic loading conditions and the irreversible changes of the geometric structures. However, it is an important issue for quantifying the realistic hydro-mechanical and chemical transport properties of fractured rocks or faults in rock engineering practices (Bandis et al., 1983; Sathar et al., 2012;Selvadurai, 2015).

(3)Coupling with shear deformation.In reality, natural rock fractures are always under a three-dimensional stress state that combines both normal and shear stresses. Extensive experimental and numerical studies based on various assumptions generally demonstrated that the shear deformation significantly affects flow and transport in rock fractures, because the shear deformation typically changes fracture structure (e.g. Li et al., 2008; Vilarrasa et al., 2011;Chen et al., 2019). However, due to challenges in geometric structure characterization and difficulties in conducting experiments in consideration of hydro-mechanical and transport coupling, the correlation between the shear stiffness and transport properties remains an open issue. A better understanding of the impact of normal deformation on the transport processes is essential for further considering the coupling with shear deformation(Selvadurai et al.,2018).

(4)Validation by experimental data.Numerous experimental studies have sought to test the hydro-mechanical couplings or solute transport processes in single rock fracture or its replica(e.g.Li et al.,2008;Zhou et al.,2015;Chen et al.,2019).High-quality and high-resolution experimental data for hydro-mechanical coupling and solute transport in realistic rock fractures are essential for validation of mechanical,flow and transport models, as well as for validation of the relationships among stress, flow and transport properties.However,reliable and comprehensive experimental data are still mostly unavailable. Nevertheless, the demonstrated important impacts of normal stress on solute transport in rock fractures presented in this study can be helpful for design of such validation experiments.

(5)Scale effects.Similar to many previous studies, the present study focused on a single fracture of limited size. An important challenging issue is scale-dependency of the stress impact on solute transport. For instance, Rutqvist(2015) found that it is important to calibrate the stress—permeability relationship against field data of stress-induced permeability changes in fractured rocks. Similarly, a better understanding of the scale effects and calibration with field data is important for application of such correlations between the hydro-mechanical and chemical transport properties in engineering practices.

5. Conclusions

In the present study, we demonstrated the potential impact of normal stress-induced fracture closure on fluid flow and solute transport in a single rock fracture, illustrating the potential correlation between the fracture normal stiffness and percentiles of the solute travel time. The findings are summarized as follows:

(1) The normal stress-induced fracture closure may have a significant impact on solute transport processes in single roughwalled rock fracture. The fracture closure reduces effective hydraulic aperture and increases the mean travel time.Relatively, the increasing normal stress enhances both the first and asymptotic arrival times for advective solute transport.

(2) The main mechanism for the stress-enhanced early arrival and tailing behavior for solute transport is that the fracture closure results in asperity contact areas that enhance channeling flow in the fracture, which in turn affects early and late arrival times.

(3) The fracture stiffness is correlated with the 5th and 95th percentiles of the normalized travel time,possibly following a linear function. This correlation between the stiffness and transport properties can be potentially useful to predict solute transport parameters from stiffness information in engineering applications, and vice versa.

(4) Without consideration of normal stresses, using simplified homogeneous models, i.e. parallel plates model, may introduce significant uncertainty in modeling fluid flow and solute transport in rock fractures,especially for the cases with higher normal stresses.

Several issues, including the impact of geometric structures,coupling with shear deformation, validation by experiments and scale effects, are essential for comprehensively understanding the hydro-mechanical and chemical coupled processes in natural rock fractures, which are important topics for future studies.

Declaration of Competing Interest

The authors wish to confirm that there are no known conflicts of interests associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

We would like to acknowledge the funding provided by the Swedish Nuclear Fuel and Waste Management Co. (SKB). We are very grateful to the reviewers for their insightful comments and valuable suggestions that are very helpful in improving this paper.