An Improved Model for Evaluating the Hydraulic Behaviour of a Single Rock Joint Considering Contact Area Evolution During Shearing

The contact of rock joints during shearing induced by underground excavation significantly impacts the hydro-mechanical behaviour of fractured rock mass, since fluid tends to flow through a rough-walled fracture along connected channels while bypassing the contact areas with tortuosity. However, previous research mostly considered joint roughness or aperture changes based on 2D joint profiles, while the contact and tortuosity using 3D surfaces were often overlooked. This paper considers the evolution of contact area and aperture distribution during shearing. The concept of the critical inclination angle is introduced and correlated with the contact area during shearing based on Grasselli’s criterion. The standard deviation of the mean mechanical aperture is calculated using the modified algorithm. An improved model for estimating the hydraulic aperture with the mechanical aperture is then developed, applying an aperture correction term and a contact correction term. A number of shear-flow tests on artificial joint samples under constant normal loads are conducted. The accuracy and reliability of the proposed model are verified by comparing it against the experimental results and by comparing the prediction performance with other published models. Results show that the proposed model can improve the prediction of the hydraulic aperture and hydraulic conductivity of a single rock joint during shearing. An improved model is proposed to estimate hydraulic aperture with mechanical aperture, incorporating aperture and contact correction terms. The evolutions of geometric morphologies during shear can be computed, and the effects on flow behaviour are well described. A series of coupled shear-flow tests are conducted on two artificial joints with different surface characteristics under constant normal loads. The proposed model can predict the shear-flow coupled behaviour of rock joints with good accuracy. An improved model is proposed to estimate hydraulic aperture with mechanical aperture, incorporating aperture and contact correction terms. The evolutions of geometric morphologies during shear can be computed, and the effects on flow behaviour are well described. A series of coupled shear-flow tests are conducted on two artificial joints with different surface characteristics under constant normal loads. The proposed model can predict the shear-flow coupled behaviour of rock joints with good accuracy.


Introduction
The presence of rock joints has a considerable impact on the strength, deformability, and permeability of fractured rock masses by providing surfaces of weakness on which further deformation are prone to occur, and by serving as the major water flow channels due to low matrix permeability (10 -20 -10 -19 m 2 ).Therefore, the understanding of water flow through rock joints has aroused more and more concerns (Dang et al. 2019;Li et al. 2019;Gui et al. 2020) for ensuring the safety and economic performance of engineering applications.
In earlier studies, water flow through rock joints is often assumed as laminar flow between two ideal smooth parallel plates, and the well-known cubic law was derived and widely used to characterize the hydraulic behaviour of rock joints.However, smooth planar fractures are unrealistic in actual geological conditions.Instead, natural rock fractures are irregular and of spatially varied aperture (Hakami 1995;Zhang and Chai 2020).According to this, apertures can generally be defined as mechanical ( e m ) and hydrau- lic ( e h ).The mechanical aperture is defined as the average vertical distance between two rock fracture surfaces and is geometrically measured such as with epoxy injection.The hydraulic aperture signifies the aperture of equivalent idealized smooth parallel plates generating the same pressure drop for a given flow rate as its original rough fractures and is measured by fluid flow analysis (Olsson and Barton 2001;Zhang and Chai 2020).Thus, the concept of the equivalent hydraulic aperture was introduced to modify the cubic law.Numerous works have been conducted on the relationship between hydraulic and mechanical aperture incorporating various geometrical parameters, as summarized in Table 1.
It has been proved by previous studies that surface roughness is one of the most sensitive factors for water flow through rock joints.The joint roughness coefficient ( JRC ) (Barton and Choubey 1977) is a widely adopted parameter to quantify the surface roughness of rock joints.Barton et al. (1985) incorporated JRC in evaluat- ing hydraulic aperture and proposed an empirical equation mainly based on normal deformation fluid flow tests.Olsson and Barton (2001) modified Barton's equation by considering the coupled shear-flow process.However, their equations do not conform to the dimensional consistency.Li and Jiang (2013) and Li et al. (2019) considered the roughness by incorporating the root mean square slope Z 2 .Renshaw (1995) and Zimmerman and Bodvarsson (1996) studied the irregularity of the joint surface by considering the effect of different aperture distributions on hydraulic conductivity.However, the effects of the aperture distribution evolution on flow behaviour during shear were not considered.Xiong et al. (2011) then improved Zimmerman's model by incorporating the evolution of fracture aperture distribution during shear and its effects on the fluid flow behaviour.However, they did not account for surface damage in the model, as the coupled shearflow tests were conducted under relatively small normal stresses (1.0 MPa and 1.5 MPa) to avoid large damage on asperities, although asperity degradation through shear is inevitable to some extent that cannot be overlooked.
The contact of rock joints during shearing induced by engineering practices also has a significant impact on the hydro-mechanical behaviour of fractured rock mass, since fluid tends to flow through a rough joint along connected channels while bypassing the contact obstacles with tortuosity.Nonetheless, previous research mostly considered joint roughness or aperture changes based on 2D joint profiles, while the contact and tortuosity based on 3D surfaces were often overlooked.Only a few investigations were made on quantifying the effect of contact area on flow behaviour.Walsh (1981) applied Maxwell's effective medium method to a fracture with randomly located circular obstructions.Zimmerman et al. (1992) extended Walsh's equation to cases where the obstacles are elliptical in shape, with random orientations, and the shapes of the irregular obstacles were described by the aspect ratio.However, the aspect ratio of irregularly shaped obstacles is hard to determine.Thus, Zimmerman and Bodvarsson (1996) incorporated the contact correction term (1 − 2c) (Kirkpatrick 1973) instead of the previous one.They also took into account the aperture variation and produced an analytical expression as the aperture correction term.Yeo (2001) conducted finite element simulations through rock fractures with conductive areas and impermeable obstacles, and modified the contact correction term to (1 − 2.4c) .Even though it yielded better estimations, this coefficient has no clear physical meaning.Furthermore, some studies (Yeo 2001;Li et al. 2008) are conducted by pre-arranging contact points between parallel joint surfaces, but these ignored contact area changes induced by shearing and the associated impact on flow behaviour.
Even though the influences of various geometrical parameters on fluid flow through rock fractures have been wellstudied by previous works, as listed in Table 1, only a few works considered the interactions of the shear-flow process.In a coupled shear-flow system, the geometrical characteristic of a fracture could be very complex since it varies with the change of normal constraints and shear displacement.Thus, there is still a lack of models to accurately describe this complicated process.c denotes the contact ratio Barton et al. (1985) e h = e 2 m JRC 2.5 JRC denotes the joint roughness coefficient.The units of e h and e m are microns Brown (1987) e Z 2 is the root mean square of the slope of the fracture profile A 0 ∕A and II represents the effect of primary and secondary roughness, respectively, C denotes the degree of heterogeneity Summarily, the quantification of fracture geometry effects on water flow in single rock joints still needs further investigation, especially when considering coupled shear-flow processes.This paper, therefore, aims to study the impacts of contact area and aperture distribution of the single rock joint on the flow behaviour.The evolutions of contact area and aperture distribution induced by progressive shearing are considered.Then, an improved model for estimating the hydraulic aperture with the mechanical aperture incorporating the above factors is developed.Finally, a series of shearflow tests are conducted on artificial joint samples, and the model performance is verified by comparing it against the results of experimental tests.

Specimen Preparation and Joint Surface Characterization
Two kinds of artificial fractures (labelled as J1 and J2) with different surface characteristics were used in tests.Both fracture surfaces were generated by splitting intact Gosford sandstone blocks, which are common construction materials in Australia, obtained from a quarry in the Sydney Basin.The sizes of two intact sandstone blocks are both 100 × 100 × 100 mm 3 .The artificial fracture J1 has a rough surface with no obvious major but plenty of small asperities, while J2 has two obvious major asperities on the surface as well as plenty of small ones.The morphology of each fracture surface is then scanned by a high-resolution 3D optical scanner system at UNSW ENG Makerspace, which can achieve precision up  to 0.04 mm, as illustrated in Fig. 1a.In addition, to avoid any tilting during shearing (Nguyen et al. 2014), and to keep a joint length of 100 mm involved in shearing, the lengths of the lower fractures are extended by 12 mm (total shear displacement), as shown in Fig. 1b, c.The scanned data are then used to digitalize the fracture surfaces and designed for reproducing the replicas.
A mixture of grout and water was used to manufacture joint replicas at a mass ratio of 1 ∶ 0.2 .The fact that replicas made from the same cast have identical fracture surfaces permits one to perform several shear-flow tests on joints with the same morphology under different normal stress conditions.The artificial fracture specimens are composed of upper and lower halves.Table 2 shows the mechanical properties of these rock-like specimens.

Asperity Degradation and Dilation
The mechanical behaviour of a single rock joint has been numerously studied over the decades.Patton (1966) assumed that asperities on the fracture surface were regular saw-tooth, and the dilation angle was described as the inclination of the teeth in an idealized model.They proposed the simplest yet widely used form of shear models for rough rock joints: = n ⋅ tan b + i .Barton and  Choubey (1977) proposed JRC and incorporated it in their shear model to describe the shear behaviour of natural rock joints with different roughness.Li et al. (2016) indicated that the shear behaviour of a joint was dominated by critical asperities with the steepest waviness facing the shear direction and the unevenness of the largest base length.
In this paper, the commonly used joint constitutive model-Barton-Bandis model (Barton and Choubey 1977;Barton 1982) is adopted for describing the asperity degradation since it is simple and can be easily upscaled to field scale.After peak shear strength, a mobilized JRC ( JRC mob ) is used to describe the asperity degradation.The joint model is expressed as: (1) where n is the normal stress, is the shear stress, r is a residual friction angle, JCS is joint wall compressive strength.Barton (1982) suggested a table for estimating the ratio JRC mob ∕JRC from the ratio s ∕ peak , where s is the shear displacement, and peak is the peak shear displacement (Asadollahi and Tonon 2010).
The dilation can be calculated accordingly with the Barton-Bandis model.However, since the main objective of this work is to develop a new prediction model for evaluating hydraulic conductivity during shearing, the dilation of the fractures will directly use experimental results.By doing this, the calibration of the shear model for predicting dilation during shearing can be omitted.

Evolution of Contact Area
The evolution of the contact area during shearing is calculated based on Grasselli's criterion (Grasselli et al. 2002;Grasselli 2006), which proposed a three-dimensional morphology characterization approach and expressed the variation of the actual contact area A * as a function of the apparent dip angle * of the surface along the shear direction.The equation is expressed as: where A 0 and * max are the maximum possible contact area ratio and the maximum apparent dip angle in the shear direction, respectively.C is a fitting parameter.The parameters A 0 , C and * max depend on the various fracture surfaces, the specified shear direction, as well as on the three-dimensional surface representation (i.e., triangulation algorithm and measurement resolution).
The concept of the threshold inclination angle is introduced to estimate the effective area involved in the flow process in a single rock joint, which is equivalent to the threshold apparent dip angle * cr for Grasselli's criterion.It is assumed that only those zones facing the shear direction with the areas of the surface inclined greater than the threshold value are potentially in contact and involved in the flow process.The maximum possible contact area ratio A 0 is calculated at a threshold angle * cr of 0°, indicating that (2) The threshold inclination angle evolution model the sum of all the areas of the surface facing the shear direction normalized with respect to the total area of the fracture.
The threshold inclination angle corresponds to the applied normal load and will be mobilized as asperity degradation during shearing.Initially, as illustrated in Fig. 2, only those zones of the surface steeper than i 0 (areas of blue in Fig. 2) are involved in shearing.As shearing progresses, the threshold inclination angle will decrease, and then all those areas with a surface steeper than i mob (areas of red in Fig. 2) will be involved in shearing.
Based on the Barton-Bandis model, the threshold inclination angle, which would be mobilized during shearing, is defined as: Thus, the relation between the mobilized threshold inclination angle i mob and mobilized potential contact area ratio c mob is calculated as: In this work, the 3D scanning point cloud data of two rough joint surfaces are processed using MATLAB code (Heinze et al. 2021) to obtain the three Grasselli's parameters A 0 , * max and C .The JRC value can be obtained from a statistical parameter formula or back-calculated from the direct shear test.Here, we adopted the back-calculation analysis from shear tests since this provides more accurate JRC values for 3D fractures.The geometrical properties of the two fracture surfaces are listed in Table 3.

Determination of Mechanical Aperture and Variable Aperture Distribution
Both normal and shear stresses change the fracture void geometry that serves as the spaces for water flow.The mechanical aperture of a rock joint, e m , can be calculated from: where e 0 represents the initial aperture at the given stress environment.Δe n is the variation of aperture induced by nor- mal loading, which is equal to 0 under constant normal load conditions.Δe s is the variation of aperture due to shearing, which is the change in the normal displacement Δ n .The initial mechanical aperture e 0 can either be obtained through the cyclic loading-unloading test, or directly using the measured hydraulic aperture before shearing (Xiong et al. 2011).
Many previous works have stated that the flow behaviour is highly sensitive to the aperture distribution (Renshaw 1995; Zimmerman and Bodvarsson 1996;Xiong et al. 2011;Xie et al. 2015;Finenko and Konietzky 2021), which is always described by the mean mechanical aperture and its standard deviation.However, the difficulty is the measurement of the variable joint aperture distribution, especially when considering the evolution of the aperture distribution during shearing.Hakami (1995) introduced several physical techniques to measure fracture aperture, including measuring the surface topography, injecting resin into a fracture and measuring the resin thickness, and casting to make a replica of the void space.But it is generally difficult to measure the aperture distribution changes during shearing using a physical measurement method.Tan et al. (2020) pasted the marking points on the surface of two fracture halves, then scanned and digitized the fracture surfaces and calculated (5) e m = e 0 − Δe n + Δe s , Fig. 3 The measurement of the variable aperture distribution the variable aperture distribution.However, these methods may have limitations to be commonly applied to other cases.Some researchers (Li et al. 2008;Xiong et al. 2011;Huang et al. 2017;Wang et al. 2022) utilized a computational method for calculating the aperture distributions with assumptions that no large damage occurred on asperities under relatively small normal stresses and gouge materials developed during shearing have negligible influence on the fluid flow.It is evident that this method is simplified since the surface damage and gauge materials could exist and have an influence on the fluid flow to a certain extent.However, the computational method can be easily applicable in various contexts.Furthermore, in our work, the influence of asperity degradation is primarily considered within the above contact treatment method.
Here, a modified computational algorithm was adopted based on an open-access source code FSAT in MATLAB (Heinze et al. 2021), which provides a method for measuring the variable aperture distribution with the input mean mechanical aperture.The algorithm was improved by considering the evolution of the mechanical aperture during shearing and calculating the standard deviation of the aperture for merely the open area without contact, as the following steps (Fig. 3): First, align the sides of the upper and lower surfaces according to the specific shear displacement s .Then, shift the upper surface relative to the lower to match the given mean mechanical aperture, which is calculated from Eq. ( 5).When the height of the lower surface is higher than the upper, the height of the lower surface will be corrected to match the upper, as illustrated in Fig. 3 a to b.Hence, the apertures of these corrected zones are zero.Finally, the standard deviation of the mean mechanical aperture can be calculated from the fractures after being shifted.It should be noted that only those non-zero apertures were considered in calculating the mean mechanical aperture and the standard deviation in our work.
The aperture distributions at three specific shear displacements ( s = 2 mm, 4 mm, and 8 mm) of two fracture surfaces (J1 and J2) measured by the modified algorithm are shown in Fig. 4. The dark blue areas in it represent the areas of zero apertures due to damage or contact.
Figure 5 shows the comparison of the aperture distribution of the fracture surface J2 at 8 mm shear displacement under n =1.5 MPa obtained from the modified computational procedure and experimental works.Details of the shear-flow experiments will be described in the later section.A thin white paint was sprayed on the surface in this case.It can be seen from Fig. 5b that the areas where the paint is worn off represent the areas of zero apertures due to damage or contact.The zero aperture areas are circled with red dashed lines in both images from the computational measurement (Fig. 5a) and experimental work (Fig. 5b).The position and shape of these areas from the two images are found to be roughly the same, indicating that the modified computational algorithm we adopted is effective for evaluating the aperture distributions.

An Improved Equation for Evaluating Hydraulic Behaviour During Shearing
When water flows through rock joints, the complexity mainly comes from the surface irregularity and the tortuosity of the flow path caused by contact areas, as illustrated in Fig. 6.
As summarized above in Table 1, Zimmerman and Bodvarsson (1996) took into account both aperture variation in the void area and contact obstacles in their equation, which adopted 1 − 1.5 2 apert ∕e 2 m as the aperture correction term, thereafter (1 − 2c) as the contact correction term.However, this equation indicates that when apert ∕e m values over 0.816, the ratio e 3 h ∕e 3 m will be less than 0 .This is unreasonable, since even in a fracture with large roughness that apert ∕e m is larger than 0.816, the flow may still be present in the fracture.Xiong et al. (2011) improved this term with 1 − 1.0 apert ∕e m through curve fitting from simulation results.They did not include the contact correction term in their equation since they stated that the contact areas with zero aperture had been taken into account.However, this is still a 2D description of the contact of the fractures, which is limited to fully capturing the complexities of 3D flow behaviour.Also, this term is presented based on a simplified Fig. 6 Schematic representation of water flow through a natural rough rock joint computational method that neglects surface damage.Thus, the contact correction term is still necessary for evaluating hydraulic behaviour, which incorporates a 3D representation of fracture contact and considers the effects of asperity degradation during shear.Yeo (2001) modified the contact correction term to (1 − 2.4c) through finite element simulations.But this coef- ficient still had no clear physical meaning.This equation indicates that when c value approaches 1/2.4, the ratio e 3 h ∕e 3 m will approach 0, and the flow will be entirely blocked off.Thus, the coefficient before c should be calculated from the maximum possible contact area in the shear direction, i.e., A 0 as mentioned above in Sect.3.2.Then the contact correc- tion term can be modified as 1 − 1 Therefore, an improved equation is proposed by incorporating two correction terms, the aperture correction term and contact correction term, which explain the reduction of flowrate by roughness and contact obstacles, respectively: where e h is the hydraulic aperture.e m is the mean mechani- cal aperture, obtained from Eq. ( 5), and e_mob is the standard deviation of the mean mechanical aperture, which will mobilize during shearing and be calculated using the computational procedure described above in Sect.3.3.A 0 is the maximum possible contact area ratio in the shear direction.Here, 1∕A 0 is given values of 2.29 and 2.22 for J1 and J2, respectively.c mob is the mobilized contact area ratio, calcu- lated from Eq. (4).c mob ∕A 0 can be defined as the normalized contact area ratio, denoted by c * mob .Thus, the proposed equation can be expressed as: where the normalized contact area ratio c * mob can be calculated from * max −i mob * max C by substituting Eq. (4) into Eq.( 6).
Furthermore, Xiong's equation is adopted as the aperture correction term.However, the improvement is that the contact areas with zero aperture were not considered in the proposed equation when measuring the variable aperture distribution, while the reduction of flow rate by contact obstacles will be corrected by the specific contact correction term.

Dimensional Analysis and Parametric Sensitivity Analysis
The proposed equation conforms to the dimensional consistency, since the hydraulic aperture e h , the mean mechanical aperture e m and its standard deviation e_mob are all in the 0.0 0 .5 1.0 1 .5 0.0 , JRC mob and c mob are dimensionless, c mob is a dimensionless ratio, and 1∕A 0 is a constant dependent on the geometric characteristic of different fracture surfaces.Furthermore, the proposed equation possesses a clear physical significance.The relationships between two correction terms with mobilized aperture distribution and mobilized contact area ratio are illustrated in Fig. 7, where f aper represents aperture correction term and f cont represents con- tact correction term.It predicts e h = e m at e = 0 and c = 0 , which recovers the case assumed in the smooth parallel plane surface.For rough joint surfaces with larger e and c , e h will decrease, since the surface roughness produces extra flow resistance and contact obstacles block off water flow, both resulting in the decrease of flow rate.When c value approaches the maximum possible contact area in the shear direction A 0 of Grasselli's criterion, the water flow will be completely blocked off, and the ratio e 3 h ∕e 3 m will approach 0. Figure 8 plots the predicted variation of e h ∕e m with the aperture parameter e ∕e m at different values of the contact area ratio c .The result implies that e h ∕e m decreases with both the increase of the aperture parameter e ∕e m and the contact area ratio c .The variation in e h ∕e m corresponding to an increase of 0.1 in c when e ∕e m = 0 is more pronounced compared with that corresponding to an increase of 0.1 in e ∕e m when c = 0 , indicating that the predicted aperture ratio is slightly more sensitive to the contact area ratio c.

Model Implementation
Figure 9 shows a flow chart for implementing the proposed model.Geometrical parameters, including three Grasselli's parameters A 0 , * max , C and the JRC value, can be readily obtained from 3D scanning point cloud data of joint surfaces by processing with the FSAT toolbox in MATLAB.It should be noted that the JRC value uti- lized in this paper is determined through back-calculation analysis based on direct shear tests.The Barton-Bandis model is adopted for describing the mechanical behaviour of rock fractures.In the post-peak stage, JRC mob is used to describe asperity degradation.The dilation can also be calculated accordingly with the Barton-Bandis model.However, in this paper, dilation of the fractures will be directly based on experimental results.
The critical inclination angle is predicted by Eq. ( 3), which will be mobilized due to asperity degradation through shear.Based on Grasselli's criterion, Eq. ( 4) continuously updates the contact area ratio in the shear direction during shearing, acting as the governing parameter in the contact correction term.Additionally, the mean mechanical aperture is updated with dilation as Eq. ( 5), and its standard deviation is updated as well using the computational method coded in MATLAB, serving as the governing parameters in the aperture correction term.The incorporation of these parameters, which are easily obtained following the treatment method explicated above, ensures the capability of the model to investigate the hydraulic behaviour of rock joints, considering the evolution of the aperture distribution and contact during shearing.Thus, the hydraulic aperture evolution during shearing can be obtained from Eq. ( 6).

Shear-Flow Experimental Apparatus
A laboratory shear-flow apparatus was adopted, which consists of a shear-flow box, a control and data acquisition system, and normal and shear actuators that are servocontrolled to apply normal and shear forces with a capacity of 500 kN and 300 kN, respectively (Fig. 10a, b).Four linear variable differential transducers (LVDTs) are placed at four corners of the upper box, and the dilations are measured by the mean values of four LVDTs.An all-around gas-pressure silicon rubber seal was adopted for sealing the water during the experiments, and two rubber side seals were used to ensure water flow through the fracture rather than through the sides of the samples, as shown in Fig. 10c-e.The applied gas pressure was 100 kPa, which was relatively small compared with the applied normal stresses (ranging from 1 to 2 MPa).Hence, the sealing system is considered to have negligible influence on the shear-flow tests.Water inlet and outlet pressure were monitored using two pressure transducers at a precision of 0.25 kPa.The pressure drop was calculated from the differential pressure of the water inlet and outlet.In 12 that curves of d= 0 mm in some cases do not conform to the general trend.This is because the fractures were not ideally well-matched at the initial shearing, and water still can be measured at the initial shearing and the shear contraction stage.These experimental data are then used to validate the performance of the proposed model, as described below.
At sufficiently low flow rates, the well-known cubic law is applicable to describe fluid flow through a single rock fracture since the inertial forces are negligible compared with viscous forces, as: where Q is the flow rate, w is the fracture width, e h is the hydraulic aperture, is the fluid viscosity, and ∇P is the pressure gradient, defined as the ratio of the pressure drops in the flow direction to the fracture length, that is ∇P = ΔP∕L.
It can be seen from Fig. 12 that the flow behaviour deviates from the linear relationship with the increase of the flow rate.This nonlinear relationship can be well described by Forchheimer's law as Eqs.( 9), ( 10), with the values of the correlation coefficient R 2 for all curves much close to 1, as listed in Fig. 12.
where A and B are the coefficients representing viscous and inertial effects, respectively, and represents the inertial resistance.
Furthermore, the Reynolds number ( Re ), which represents the nonlinearity in the fluid flow, is defined as the ratio of ( 8)  where is the fluid density, v is the average flow velocity, and D is the characteristic dimension of the flow system.The transmissivity ( T ) is used to describe the conductiv- ity of a fracture.T 0 represents the transmissivity when fluid flow in the linear regime.Thus, the normalized transmissivity ( T∕T 0 ) is commonly used to describe the nonlinear flow behaviour in fractures, defined as: The normalized transmissivity ( T∕T 0 ) can also be described as a function of Re as: where is the Forchheimer coefficient.The relationships between T∕T 0 and Re for two fractures under different nor- mal loads are shown in Fig. 13.It can be seen that T∕T 0 decreases with an increase in Re , which reveals the nonlin- earity of the flow.T∕T 0 = 0.9 indicates that the nonlinear (11) term ( BQ 2 ) contributes to 10% of the pressure drop, the Reynolds number at this point is referred to as the critical Reynolds number ( Re c ) (Yu et al. 2017), as shown by the dashed lines in Fig. 13.The flow in our tests shows obvious nonlinear behaviour.In addition, the T∕T 0 -Re curves gener- ally shift downwards as the shear displacement d increases, indicating that the nonlinearity of fluid flow is stronger at larger shear displacement.However, the investigation of the nonlinear flow equation is out of scope in this study.

Correlation with Experimental Results
The fitted quadratic polynomial expression in the form of Forchheimer's law consists of a linear term ( AQ ) and a non- linear term ( BQ 2 ) describing the viscous and inertial pres- sure drops, respectively, as shown in Fig. 14.Chen et al. (2015) indicated that the hydraulic aperture e h , representing the equivalent aperture of the parallel-plates model, can be back-calculated by substituting the slope of the linear regression line fitted to the linear portion of ∇P versus Q curves into the cubic law (Eq.( 8)), that is The measured e m and e h at different shear displacements were calculated with Eqs. ( 5) and ( 8).Here, e 0 is assumed to be equal to the initial hydraulic aperture back-calculated from the measured flow rate at the initial shearing by the cubic law.Based on the fracture surface scanning data and the above flow chart, the calculated e h was obtained from Eq. ( 6).The contact correction term should be (1 − 2.29c mob ) for fracture J1 and (1 − 2.22c mob ) for fracture J2.
Figure 15 shows the measured mechanical aperture (solid line) and the measured hydraulic aperture (dash line) versus shear displacement.It can be seen that the increase of the mechanical aperture is more significant than that of the hydraulic aperture after the peak shear stress (at around 1-2 mm shear displacement).Furthermore, comparisons of the predictions using the proposed model Eq. ( 6) (dash line) and the experimental results (hollow scatters) for two fractures under different normal stress levels are demonstrated in Fig. 15.The proposed model shows good agreement with the experimental results.It should be noted that each fracture specimen had initial apertures to varying degrees since the upper and lower fracture specimens always had a mismatch somewhat.Even though we tried to achieve a fully mated for the upper and lower fracture specimens before shearing in the laboratory tests, it is hard to ensure the fully mated state since the shear box is invisible.For a more straightforward comparison of the results under different normal stresses To further verify the proposed equation quantitatively, the estimations from the proposed model (Eq.( 6)) are compared with the results from other references.Even though extensive research on fluid flow through a single rough joint has been carried out over the past several decades and many equations for predicting the hydraulic aperture were proposed, few of them considered the complicated shear-flow coupled process.Olsson and Barton (2001)  It should be noted that Xiong's model and the proposed model both adopted the standard deviation e of the mean mechanical aperture to represent the variable aperture distribution.However, when determining e , Xiong's model considers the zero aperture as contact zones, whereas the proposed model only counts the void zones of the fractures, with the effect of contact zones being calculated with another term that incorporates the 3D representation of fracture contact and the consideration of the asperity degradation during shear.
Comparison between the calculated results with different equations and the measured hydraulic apertures are presented in Fig. 17.The solid line represents that the calculated e h from different equations is equal to the measured e h from experiments, indicating that the equation has per- fect capacity for predicting hydraulic conductivity.It can be obviously seen from Fig. 17a that Olsson and Barton's model significantly underestimates the hydraulic apertures.Figure 17b shows that Xiong's model displays roughly good predictive performance for experimental data, but there is a moderate overestimation.This discrepancy can be attributed to the omission of roughness degradation consideration.This also indicates that considering only the zero aperture in two dimensions is insufficient to describe the contact effects during shear.Therefore, the contact correction term that incorporates a 3D representation of fracture contact and accounts for asperity degradation during shear is necessary to improve the predictive capability of the model.The proposed equation shows better accuracy since the data points are all located near the solid line in Fig. 17c.The average estimation errors are further used to evaluate the prediction accuracy and can be obtained with Eq. ( 16): where e mea,i h is the measured hydraulic aperture, e cal,i h is the calculated hydraulic aperture, i represents the i th data point, and n represents the total number of the dataset we obtained.The average estimation error is 4.34% for the proposed model, while the average estimation errors are 61.23% for Olsson's equation and 10.16% for Xiong's equation.One can see that Olsson's equation displays a huge estimation error for experimental data.This may be caused by the dimensional inconsistency of the equation as Eq. ( 14).Xiong's equation shows good predictive performance as well.However, the proposed model in this study is still more accurate than Xiong's equation.
In summary, the above analysed and discussed the correlations between the calculated hydraulic apertures using the proposed model in this study and the measured hydraulic apertures from experiments, as well as the comparisons between the prediction performance of the proposed model with other published models.Results show that the proposed equation can provide a reliable prediction for the hydraulic aperture and hydraulic conductivity of a single rock joint.

Correlation with Experimental Data from Literature
To further verify the applicability of the proposed model, comparisons are made between the predictions from the proposed model and the experimental results performed by Wang et al. (2020b).They conducted a series of shear-flow tests on artificially created rough-walled fractures (Labelled as G3) with the surface morphology shown in Fig. 18.The size of the fracture surface is 200 × 100 × 100 mm, for which the JRC value of 7.36 was reported.The geometrical properties of G3 based on Grasselli's criterion are listed in Table 4. Therefore, the constant 1∕A 0 of the proposed model is given a value of 1.91 for G3.Plaster replicas were used to conduct the tests with the uniaxial compression strength c was 38.5 MPa, the tensile strength t was 2.5 MPa, and the internal friction angle of the joint surface was 60° (Wang et al. 2020a).
Figure 19 illustrates the comparison between predictions by the proposed model and experimental results by Wang et al. (2020b) in two cases on fracture surface G3.Although there are a few deviations at the larger displacement of case 1, the results show that the proposed model exhibits acceptable agreement with the experimental results reported in the literature, indicating that it has good applicability.

Conclusions
In this work, an improved model for estimating the hydromechanical behaviour of the rough-walled single rock joint is proposed.The evolutions of the variable aperture distribution and contact area ratio during shearing are considered to evaluate the hydraulic aperture, of which the variable aperture distribution changes can be calculated with a modified computational procedure, and the contact area ratio evolutions are calculated with the critical inclination angle changes based on Grasselli's criterion.The proposed model conforms to the dimensional consistency, while the model parameters are geometric parameters with a clear physical significance and can be readily obtained.
A total of 336 shear-flow testing data on artificial joint samples under constant normal loads are employed to verify the performance of the proposed equation for estimating hydraulic aperture.The experimental results of our tests are well captured by the proposed equation, suggesting that it can be used to effectively evaluate the hydraulic conductivity of a single rock fracture.Furthermore, comparisons between the proposed model and the existing models are also carried out based on the experimental results.The results show that the proposed model displays better accuracy and can provide a reasonable estimation of the shear-flow coupled behaviour of rock joints, thereby improving the reliability of stability and safety predictions for rock masses in engineering applications.
2c) apert denotes the standard deviation of the mean mechanical aperture

Fig. 1 a
Fig. 1 a The 3D optical scanner system, b the scanned rock joint surface J1, and c the scanned rock joint surface J2

Fig. 4 Fig. 5
Fig. 4 Aperture distributions at three specific shear displacements of two fracture surfaces measured by the modified algorithm

Fig. 7 Fig. 8
Fig. 7Relationships between two correction terms with mobilized aperture distribution and contact area ratio, respectively

Fig. 9
Fig. 9 Flowchart for the implementation of the proposed model

Fig. 10
Fig. 10 The coupled shear-flow experimental apparatus: a front view, b back view, c schematic diagram of the sealing system inside the shear box, d the seals inside the shear box, and e the all-round silicon rubber sealing after gas injection

Fig. 11
Fig. 11 Mechanical behaviours of fractures J1 and J2 in coupled shear-flow tests

Fig. 14
Fig. 14 Determination of the hydraulic aperture

Fig. 16 Fig. 17 c
Fig. 16 3D plot of the relationship between e 3 h ∕e 3 m and two correction terms (i.e., 1 − 1.0 e_mob e m and (1 − 1 A 0 c mob ) ) obtained from the proposed model and experimental results modified Barton's equation by incorporating the surface roughness JRC for estimating the hydraulic aperture considering the shear-flow coupled process.Xiong et al. (2011) improved Zimmerman's model which adopted the standard deviation e to represent the variable aperture distribution.Their equation takes account of the evolution of fracture aperture distribution during shear and its effects on the fluid flow behaviour.Thus, the equations proposed byOlsson and Barton, and  Xiong are adopted as comparisons with our proposed equation:Olsson and Barton's equation: JRC mob , s ≥ spXiong's equation: G3 under σ n0 = 1 MPa, k n = 1 MPa/mm

Table 1
Relationship between hydraulic aperture ( e h ) and mechanical aperture ( e m ) from previous works JRC mob , s ≥ sp JRC mob is the mobilized value of JRC , s is the shear displacement, and sp is the peak shear displacement

Table 2
Mechanical properties of fracture specimens

Table 3
Geometrical properties of two fracture surfaces based on Grasselli's criterion

Table 4
Geometrical properties of G3 based on Grasselli's criterion