ThermalHydraulicMechanical (THM) Coupling Behaviour of Fractured Rock Masses
View this Special IssueResearch Article  Open Access
Hao Kang, Herbert Einstein, Stephen Brown, John Germaine, "Numerical Simulation for Rock Fracture Viscoelastic Creep under Dry Conditions", Geofluids, vol. 2020, Article ID 8879890, 30 pages, 2020. https://doi.org/10.1155/2020/8879890
Numerical Simulation for Rock Fracture Viscoelastic Creep under Dry Conditions
Abstract
In many rock engineering projects such as hydrocarbon extraction and geothermal energy utilization, the hydraulic and mechanical behavior of rock fractures significantly affects the safety and profitability of the project. In field conditions, the hydraulic and mechanical behavior of rock fractures changes with time (the rock fractures creep), and this creep is not negligible even under dry conditions. In addition, creep is strongly affected by the rock fracture surface geometry. However, there is not much literature systematically studying the effect of surface geometry on rock fracture creep under dry conditions. This paper presents the results of a numerical study considering the effect of surface geometry on rough fracture viscoelastic deformations. An inhouse numerical code has been developed to simulate the viscoelastic deformation of rough fractures. In addition, another numerical code has been developed to generate synthetic rough fracture surfaces by systematically changing the surface roughness parameters: the Hurst exponent, mismatch length, and root mean square roughness. The results indicate that by increasing the Hurst exponent or decreasing the mismatch length or decreasing the root mean square roughness, the fracture mean aperture decreases, and the contact ratio (number of contacting cells/total number of cells) increases faster with time.
1. Introduction
In many rock engineering projects such as hydrocarbon extraction, geothermal energy utilization, or nuclear waste storage, the hydraulic and mechanical behavior of rock fractures has a strong effect on the safety and profitability of the project. In field conditions, under constant tectonic stresses, rock fractures may creep. Creep, the deformation with time under constant stress, occurs in many types of rocks. The creep of intact rock has been studied extensively, and some authors proposed sophisticated creep models by conceptualizing the rock creep as a combination of springs and dashpots [1–4]. Also, creep rate (i.e., the change of deformation with time) may accelerate, stay constant, or decelerate. Such creep mechanisms often occur in rock fractures, essentially involving a timedependent change of aperture under constant load. Since creep, in time, affects hydraulic conductivity and thus affects the safety and profitability of rock engineering projects, there is significant interest in understanding the creep behavior of rock fractures.
Many papers have discussed numerical modeling and laboratory experiments on creep mechanisms of rock fractures. The creep of rock fractures containing fluids is a process of four interdependent mechanisms: mechanical compression, pressure solution, dissolution, and erosion. It appears that most experiments and numerical models focus on the mechanisms of pressure solution (at contacting asperities) and dissolution (on free surfaces). For dissolution, fracture permeability evolution under the influence of acidic solutions has been widely investigated [5–11]. A common conclusion is that for calcite rock, when the pH of the fluid in the fracture is less than 6, fracture aperture distribution becomes more heterogeneous, as a result of acidic dissolution. The effect of pressure solution on rock fracture creep has also been studied extensively [12–16]. A common finding is that pressure solution strongly affects the fracture mechanical and hydraulic aperture, especially for siliciclastic rock. In recent years, some numerical models, which simulate fracturescale pressure solution, have been proposed [17–23].
Based on the literature review on rock fractures in general, there is not much work systematically studying the creep mechanism caused by mechanical compression. Kang et al. [24] reported that for Musandam limestone fractures, the creep deformation under dry conditions is not negligible. Under dry conditions, mechanical compression is the only creep mechanism suggesting that it should be systematically studied. In addition, previous experiments suggest that the fracture surface geometry has a dominant effect on the hydraulic and mechanical behavior of rock fractures [24–27]. Therefore, the effect of surface geometry on fracture viscoelastic deformation should be systematically investigated.
Brown [28] proposed a simple mathematical model to geometrically describe rock fracture surfaces. In this model, a rock fracture surface can be completely described by three parameters: the Hurst exponent, the mismatch length, and the root mean square roughness. In this research, Brown’s [28] model will be used to generate synthetic fracture surfaces: the three abovementioned parameters will be varied systematically. In tribology, Chen et al. (2011) first numerically simulated the viscoelastic deformation of rough fracture surfaces. In this research, their method will be used to numerically determine the viscoelastic deformation for each surface geometry parameter combination.
The objective of this numerical study is to investigate the effect of surface geometry on rock fracture viscoelastic deformations. First, the method for surface geometry generation and viscoelastic deformation calculation will be explained in Sections 2.1–2.3. Then, the numerical results will be validated against analytical results for a given surface geometry in Section 2.4. Next, the viscoelastic deformation of different rough fracture surfaces (generated based on Brown’s model) will be compared and discussed in Sections 3 and 4. Finally, recommendations will be made in Section 5. The results of this numerical study will eventually provide a reference and basis for fracture viscoelastic deformation prediction.
2. Methodology
2.1. Methodology for Rough Fracture Surface Generation
Synthetic fracture surface pairs are generated based on Brown’s probabilistic model [28]. In this model, rough rock fracture surfaces can be approximated by a Gaussian height distribution and selfaffine organization over a large range of length scales ([28–30]). The rock fracture surface can be described by three parameters: Hurst exponent (), root mean square roughness at a reference length scale (), and mismatch length scale () [28].
The surface roughness of rock fractures typically follows a selfaffine fractal distribution ([29, 31]). A selfaffine surface can be described as where is the surface elevation, is a constant for direction scaling, and is the Hurst exponent. The value is between 0 and 1; and a larger value corresponds to a smoother local surface profile.
can be well quantified from the power spectrum of the surface. The power spectrum is calculated by decomposing the surface into a sum of sinusoidal waves, each with its amplitude, phase, and wavelength. Figure 1 shows the process of wave decomposition. In the wave decomposition process, the Fourier transform is used to decompose the random rough surface into a series of sinusoidal wave components. The squared amplitude of each wave component () is defined as the power, and a plot of power versus the inverse of wavelength () is defined as the power spectrum (see Figure 2 for a schematic). For selfaffine fracture surfaces, the power spectrum can be related to the wavenumber ():
The upper frequency end, , is defined as , where is the smallest representable or measurable size. The lower frequency end, , is defined as , where is the surface length.
The second parameter is the mismatch length . As shown in Figure 1, each waveform can be completely described by three parameters: wavelength (), amplitude (), and phase. In a natural rock joint, the two surfaces may not be perfectly matched; they may have relative shear displacements [28–30]. Rock fracture surface measurement results from Brown [28] and Matsuki et al. (1995) suggest that natural fractures are matched at long wavelengths, but the surfaces are not identical at short wavelengths. Based on these observations, Brown [28] proposed the parameter of critical wavelength , called the mismatch length scale. Above the mismatch length scale, the decomposed sinusoidal wavelengths are perfectly matched; they have the same phase, wavelength, and amplitude. Below the mismatch length scale, the decomposed sinusoidal wavelengths behave independently; they have identical relationship between and , but independent phases. Figure 3 illustrates the concept of the mismatch length scale.
(a)
(b)
The third parameter is the root mean square roughness (). This parameter is used to define the absolute scale of the surface height. can be expressed as
After generating the surface profile, its root mean square value, , is calculated. Then, the surface height is scaled linearly as where is the initial surface height, is the final surface height, and is the designated root mean square roughness. Figure 4 illustrates the concepts of (surface height) and (aperture).
In this research, rough fracture surface pairs will be generated based on Brown’s [28] model. The Hurst exponent, root mean square roughness, and mismatch length scale will be varied systematically to investigate the effect of surface roughness parameters on fracture viscoelastic deformations.
2.2. Methodology for Calculating Fracture Elastic Deformation
The method for calculating fracture elastic deformation is introduced in this section. An inhouse contact simulation code has been developed, and the method is similar to that formulated by Polonsky and Keer [32]. Here, only the key mathematical aspects will be presented; the details can be found elsewhere [32]. Polonsky and Keer’s method has been widely used in tribology and rock mechanics.
First, the aperture distribution function is defined as where is the rigidbody movement between two surfaces under applied stress in the direction normal to the fracture surface (compression is defined as positive), is the initial aperture, and is the elastic displacement field between the two contacting surfaces inside the contacting area.
For two rough surfaces contacting each other, the boundary conditions below need to be satisfied: where is the contact pressure (in the normal direction) acting on location . This indicates that the contact pressure is zero at noncontacting areas and is larger than zero at contacting areas. Compressive stress is defined as positive.
The vertical elastic displacement can be calculated from the Boussinesq and Cerrutti (1885) equation: where is the influence matrix for normal displacement at surface location due to a normal unit stress applied at location and is the normal stress applied at location . Figure 5 shows the schematic of the Boussinesq and Cerrutti equation.
The influence matrix is defined as where and are the shear modulus and Poisson’s ratio, respectively.
The elastic displacement of two rough surfaces with partial contact usually cannot be solved analytically, but can be solved numerically. For numerical simulation, the analytical equations need to be discretized into numerical equations. First, as shown in Figure 6, the surface is discretized into a regular grid: where and are and coordinates, respectively; and are the grid spacings in the and directions, respectively; and and are the total number of grids in the and directions, respectively.
The aperture field can be calculated as
The boundary conditions become
Love [33] discretized the Boussinesq and Cerrutti equation as where
Figure 6 illustrates the indices.
Andrews [34] first used the fast Fourier transform (FFT) method to solve the Boussinesq and Cerrutti equation in rock fractures. Later, Stanley and Kato [35] used this method in tribology. Fast Fourier transform is an algorithm that calculates the discrete Fourier transform of a sequence. In engineering, people often use Fourier analysis to convert a signal from its original domain (space or time) to a representation in the frequency domain. Discrete Fourier transform is obtained by decomposing the sequence of values into components of different frequencies, and using the fast Fourier transform can save the computation time. The detailed explanation of the FFT method can be found in Stanley and Kato [35, 36]. Using the FFT method, can be expressed as where FFT^{1} means inverse fast Fourier transform (IFFT). In mathematics, FFT changes the convolution (see Equation (13)) to a simple multiplication (see Equation (16)).
If the FFT method is not used, calculating Equation (10) requires operations, while by using the FFT method, calculating Equation (13) requires operations. Therefore, when and are large, using the FFT method can save computation time.
In addition, force balance needs to be achieved: where is the applied external force on the two fracture surfaces and represents the pressure acting on grid (see Equation (12) and Figure 5). The indices and are explained in Figure 5.
Equations (10), (11), (12), (16), and (17) are solved iteratively using the conjugate gradient (CG) algorithm proposed by Polonsky and Keer [32], where the details can be found.
2.3. Methodology for Calculating the Viscoelastic Deformation
Chen et al. (2011) [37] proposed a numerical method to solve the viscoelastic deformation of rough fracture surfaces. In this research, an inhouse viscoelastic contact simulation code has been developed, and the method is similar to that formulated by Chen et al. (2011). Here, the key mathematical concepts will be presented. Other details can be found in Chen et al.’s work (2011).
Before describing the methodology, it is essential to review the basic concepts of viscoelasticity. Viscoelastic responses can be modeled by a group of dashpots and springs, in which the springs represent elasticity and dashpots represent viscosity. In this research, we assume that the material is linearly viscoelastic. This means that the stress/strain response scales linearly with the stress/strain input and follows the law of linear superposition. For any arbitrary series of stress or strain inputs, the outputs can be expressed as where and are the relaxation modulus function and creep compliance function, respectively. (unit: pressure) describes the timedependent stress response to a stepchange in strain, and (unit: 1/pressure) describes the timedependent strain response to a stepchange in stress. represents the time variable in the integration, and its range is (this means that is a constant and is the variable). Figure 7 explains the creep compliance and relaxation modulus.
(a)
(b)
The analytical solution for the Boussinesq and Cerrutti equation can be modified by adding timedependent viscoelasticity. Adding the creep compliance function into the Boussinesq and Cerrutti equation (Equations (7) and (8)) yields
In Equation (19), the shear modulus (in Equation (8)) is replaced by the creep compliance .
Rearranging Equation (19) yields
Equation (20) cannot be solved analytically. However, it can be solved numerically, if it can be decomposed into a linear equation system. Compared with Equation (7), an additional integration over the time duration () is added in Equation (23). Therefore, the time duration also needs to be discretized into many small time steps. Here, the total time duration is discretized into time steps with a uniform interval . It is assumed that is sufficiently small that the pressure field at each time step is unchanged. In addition, according to the fundamental theorem of calculus, the term in Equation (20) can be expressed as . Based on the sufficiently small time step assumption and the fundamental theorem of calculus, Equation (20) can be rewritten as a summation: where .
Within each time duration, in a given surface, the pressure field does not change with time. Therefore, can be factored out from the integral. Equation (21) becomes
The integral part of Equation (22) can be discretized. From Section 2.2 (Equations (7), (8), and (13)), the following discretization relationship can be obtained:
By multiplying both sides of Equation (23) by , the lefthand side of Equation (26) becomes the integral part of Equation (22). Therefore, the integral part of Equation (22) can be discretized as where is the shear modulus.
Thus, Equation (22) can be rewritten as
Equation (25) can be decomposed into Equations (26) and (27):
It is worth noting that Equation (27) is analogous to Equation (13). Therefore, to save computation time, Equation (27) can also be solved using the FFT method, similar to Equation (16): where , , and in Equation (28) are equivalent to , , and in Equation (18), respectively.
Equations (26) and (28) will be used for numerical simulation. At each specific time step, Equations (10), (11), (12), (17), (26), and (28) will be solved iteratively using the CG method [32], and the pressure and displacement fields at each step can be obtained. Then, will be increased (a new time step will be added), and the new pressure and displacement fields will be solved.
Based on the above mathematical concepts, the numerical algorithm is developed. The main algorithm is summarized below: (1)Obtain the loading history , initial discretized aperture field , and the creep compliance function for the viscoelastic material(2)Calculate for each time step and the stiffness matrix (3)Use the method described in Section 2.2; calculate the initial contact state and the pressure distribution (). In this case, there is no creep effect, and the stress and displacement fields are calculated based on elasticity(4)Save the historical pressure distribution field(5)Now increase (apply a new time increment). Since the historical pressures are known, the pressure distribution field and displacement field at time step can be calculated by the CG method(6)If reaches (the total number of time steps), stop the program, and export the pressure and displacement field. Otherwise, go back to step 4
Figure 8 shows the flowchart of the main algorithm [38].
2.4. Validation of the Numerical Program
The results of the numerical program need to be validated against analytical solutions. In this research, the cases for validation will be the Maxwell and Standard Linear Solid (SLS) viscoelastic halfspaces indented by a spherical rigid body. Lee and Radok [39] developed the analytical solutions for the two abovementioned cases, and the pressure distribution along the contacting region will be compared.
It is essential to explain the concept of the Maxwell and Standard Linear Solid (SLS) constitutive models. Figures 9(a) and 9(b) show the schematic for the Maxwell model and SLS model, respectively.
(a)
(b)
The Maxwell model consists of a spring and a dashpot. The spring represents elasticity, with a shear modulus of ; and the dashpot represents viscosity, with a viscosity of (unit: time pressure). Under constant stress, the strain can be expressed as
Equation (29) indicates that under constant stress, the strain rate (creep rate) is constant.
The creep compliance function can therefore be defined as where is the time. The relaxation time for the Maxwell model, , is defined as
The SLS model consists of two springs and one dashpot. The two springs represent elasticity, with an elastic modulus of and , respectively; and the dashpot represents viscosity, with a viscosity of (unit: time pressure). Under constant stress , the strain can be related to as
The creep compliance function can be expressed as where is the time. The relaxation time for the SLS model, , is defined as
In the numerical code for the Maxwell model, Equation (30) is implemented in Equation (25), and the pressure and displacement fields are solved as described in Methodology; for the SLS model, Equation (33) is implemented in Equation (25), and the pressure and displacement fields are solved as described in Methodology.
Lee and Radok [39] stated that if a spherical rigid body is indenting into a viscoelastic flat surface (see Figure 10) and the flat surface satisfies the Maxwell model, the pressure distribution field can be expressed as where is the pressure, is the distance from the center of the contacting region, is the total time duration, is the radius of the contacting region, is the spherical diameter, is the shear modulus, is the viscosity, and is the total load. Figure 10 illustrates the physical parameters.
Lee and Radok [39] also suggested that if the flat surface satisfies the SLS model, the pressure distribution field can be expressed as where is the pressure, is the distance from the center of the contacting region, is the total time duration, is the radius of the contacting region, is the spherical diameter, and are the shear moduli, is the viscosity, and is the total load. The physical parameters are also illustrated in Figure 10.
Johnson (1985) [40] solved Lee and Radok’s analytical equations. Figures 11 and 12 compare the analytical and numerical solutions for the Maxwell and SLS models, respectively. The analytical solutions were calculated by Johnson (1985), and the numerical solutions were calculated by the code developed in this research. In the comparison, the absolute values of and do not matter, because the pressure distribution and the radius of the contacting region are normalized.
The results in Figures 11 and 12 suggest that the difference between the numerical and analytical solutions is less than 10%. Therefore, the numerical results match well with the analytical results.
3. Results
3.1. Surface Generation Results
As described in Section 2.1, synthetic rough fracture surfaces are generated based on Brown’s model [28]. Brown [28] analyzed the surface profiles of 23 natural rock joints and calculated the Hurst exponent (), root mean square roughness (RMS), and mismatch length () of each rock joint. The rock types include tuff, granodiorite, granite, basalt, chalk, and siltstone. Table 1 summarizes his analysis results. RMS and are normalized by the profile length.

Based on the results in Table 1, it can be concluded that the Hurst exponent is normally between 0.50 and 0.90; the normalized RMS () is normally between 0.005 and 0.015; and the normalized () is normally between 0.02 and 0.2. The synthetic surface pairs described in Methodology (Section 2.1) are generated based on the above values. Specifically, seven pairs of the synthetic surfaces are generated with three different , RMS, and values. Table 2 summarizes the information of the seven pairs of the synthetic surfaces. For all the surfaces, the dimensions (profile lengths) in the and directions are 10 mm.

Table 2 shows that in surface pairs 1, 2, and 3, the value is varied; in surface pairs 2, 4, and 5, the value is varied; in surface pairs 2, 6, and 7, the RMS value is varied.
For each pair of the surfaces, the aperture field can be calculated. Figures 13–15 show the effect of different , , and RMS values on the aperture field, respectively. The following conclusions can be drawn: (1)From Figure 13, as decreases, the mean and standard deviation of the aperture increase(2)From Figure 14, as increases, the mean and standard deviation of the aperture increase(3)From Figure 15, as RMS increases, the mean and standard deviation of the aperture increase. The magnitude of the aperture field scales linearly with the RMS value (this can be explained from Equation (4) in Section 2.1)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
Table 3 summarizes the mean and standard deviation of the aperture field for the seven surface pairs. Each of the calculated aperture field is implemented into the creep simulation code as the initial aperture field. The creep deformation of the aperture field is calculated as described in Sections 2.2 and 2.3.

3.2. Creep Simulation Results for the Maxwell Model
The creep deformations of the seven surface pairs are calculated using the Maxwell model. Table 4 summarizes the input parameters. The parameters are obtained from the testing results of the Vaca Muerta Shale published by Mighani et al. [41].

The following parameters are also defined: (1)Macroscopic stress , where is the total force applied on the fracture surface and is the area of the fracture surface ()(2). See the definition in Figure 6
Figure 16 summarizes the effect of different , , and RMS values on the creep deformation. In the simulation, is fixed at 10 MPa. The initial values of the average aperture and contact area correspond to the elastic deformation. The time duration is from 0 to .
(a) Mean aperture
(b) Contact ratio
The following conclusions can be drawn: (1)Based on the curves of surface pairs 1, 2, and 3, when increases, the average aperture decreases, and the contact ratio increases faster with time(2)Based on the curves of surface pairs 2, 4, and 5, when decreases, the average aperture decreases, and the contact ratio increases faster with time(3)Based on the curves of surface pairs 2, 6, and 7, when RMS decreases, the average aperture decreases, and the contact ratio increases faster with time(4)Under current surface parameters, macroscopic stresses, and the time durations, the contact ratio is generally less than 10%
Table 5 summarizes the above conclusions.
 
Note: “↓” means decrease and “↑”means increase. 
Figure 17 compares the contact area and local contact stress before and after the creep stage. The white areas correspond to noncontacting areas, and the colored areas correspond to contacting areas. The input parameters are , , , and , and the creep . The left figure corresponds to the stress distribution before the creep (the moment after the elastic deformation), and the right figure corresponds to the stress distribution after the time duration of . The scale of the color bar is 0 to 2000 MPa. After creep, the contact area increases, and the local contact stresses (the vertical stress in the contacting cells, see the definition of cell in Figure 6) decrease (the color becomes darker).
(a)
(b)
3.3. Creep Simulation Results for the Standard Linear Solid (SLS) Model
The creep deformations of the seven surface pairs are also calculated by the SLS model. Table 6 summarizes the input parameters. The parameters are obtained from the testing results of the Vaca Muerta Shale published by Mighani et al. [41].

Figure 18 summarizes the effect of different , , and RMS values on the creep deformation. In the simulation, is fixed at 10 MPa. The initial values of the average aperture and contact area correspond to the elastic deformation. The time duration is from 0 to . Compared with the Maxwell model, the time duration is extended from to to better illustrate the decreasing creep rate.
(a) Mean aperture
(b) Contact ratio
The following conclusions can be drawn: (1)Based on the curves of surface pairs 1, 2, and 3, when increases, the average aperture decreases, and the contact ratio increases faster with time(2)Based on the curves of surface pairs 2, 4, and 5, when decreases, the average aperture decreases, and the contact ratio increases faster with time(3)Based on the curves of surface pairs 2, 6, and 7, when RMS decreases, the average aperture decreases, and the contact ratio increases faster with time(4)Under current surface parameters, macroscopic stresses, and time durations, the contact ratio is generally less than 10%(5)Under current surface parameters, macroscopic stresses, and time durations, the creep rate (aperture decrease rate) decreases significantly with time. The main reason is that the SLS model assumes exponentially decreasing creep rate
Table 7 summarizes the above conclusions. The effects of , RMS, and on the creep rate under the SLS model are similar to those under the Maxwell model.
 
Note: “↓” means decrease and “↑” means increase. 
Figure 19 compares the contact area and local contact stress before and after the creep stage. The white areas correspond to noncontacting areas, and the colored areas correspond to contacting areas. The input parameters are , , , and , and the creep time . The left figure corresponds to the stress distribution before the creep (the moment after the elastic deformation), and the right figure corresponds to the stress distribution after the time duration of . The scale of the color bar is 0 to 2000 MPa. After creep, the contact area increases, and the local contact stresses (the vertical stress in the contacting cells, see the definition of cell in Figure 6) decrease (the color becomes darker).
(a)
(b)
4. Brief Discussion
4.1. Creep Deformation Normalization: The Maxwell Model
The creep curves (mean aperture and contact ratio versus normalized time) under different , , and RMS values can be normalized. For each surface pair (see the definition in Sections 2.1 and 3.1; each surface pair generates one aperture field), between time duration 0 and , the mean aperture, , and the contact ratio, , are normalized as where is the normalized mean aperture, is the mean aperture right after the elastic deformation (before creep), is the mean aperture right after the elastic deformation for surface pair 2 (, , and ), is the normalized contact ratio, is the contact ratio right after the elastic deformation (before creep), and is the contact ratio right after the elastic deformation for surface pair 2.
Figure 20 summarizes the normalized mean aperture and contact ratio for different , , and RMS values. is fixed at 10 MPa, and the time duration is from 0 to .
(a) Normalized mean aperture
(b) Normalized contact ratio
The results show that by normalizing the values as described in Equations (37) and (38), the curves for different surface parameters fall into a narrow region.
4.2. Creep Deformation Normalization: The SLS Model
The creep curves (mean aperture and contact ratio versus normalized time) under different , , and RMS values can be normalized similarly as described in Section 4.1. For each surface pair, between time duration 0 and , the mean aperture, , and the contact ratio, , are normalized as where is the normalized mean aperture, is the mean aperture right after the elastic deformation (before creep), is the mean aperture right after the elastic deformation for surface pair 2 (, , and ), is the normalized contact ratio, is the contact ratio right after the elastic deformation (before creep), and is the contact ratio right after the elastic deformation for surface pair 2.
Figure 21 summarizes the normalized mean aperture and contact ratio for different , , and RMS values. is fixed at 10 MPa, and the time duration is from 0 to .
(a) Normalized mean aperture
(b) Normalized contact ratio
The results show that by normalizing the values as described in Equations (39) and (40), the curves for different surface parameters fall into a narrow region.
4.3. Contact Pressure Evolution during the Creep Stage: One Example
Figure 22 compares the histogram of local contact stress (only for the contacting cells, see the definition of cell in Figure 6) before and after the creep stage. The horizontal axis represents the local contact stress in each contacting cell, and the vertical axis represents the normalized frequency (the number of contacting cells within the stress range/the total number of contacting cells). The SLS model is used as the viscoelastic model, and the time duration is . Other input parameters are , , , and . The blue histogram corresponds to the local stress histogram before the creep (the moment right after the elastic deformation), and the orange histogram corresponds to the local stress histogram after the time duration of . This plot implies that after creep, the local stresses in the contact region decrease (the histogram shifts leftwards).
From the histograms, it is worth noting that the local contact stress can be larger than 200 MPa. However, for both the Maxwell and SLS models, the input elastic modulus and viscosity are obtained from triaxial tests with confining and differential pressures less than 50 MPa. Therefore, if possible, creep data from triaxial tests with confining and differential pressures larger than 200 MPa should be used to obtain more accurate simulation results.
4.4. Limitations for the Numerical Method
In this simulation, the contacting asperities deform viscoelastically. This implies that there is no upper limit for the compressive stress in the contacting cells. For some synthetic surface pairs, the compressive stress in a few contacting cells reaches 1.5 GPa. In reality, under such high compressive stresses, the cells may deform plastically, and the contact ratio will further increase. Therefore, if plastic deformation is ignored, the compressive stress in the contacting cells may be overestimated, and the contact ratio may be underestimated. In addition, asperity breakage or cracking is ignored. In reality, under high compressive stress, the asperities may break, which will affect the contact ratio and contact stress distribution.
5. Conclusions and Recommendations
In this research, an inhouse numerical code has been developed to generate synthetic rock fracture surfaces based on Brown’s [28] model, and another code has been developed to simulate viscoelastic deformations of rough fractures. Seven synthetic fracture surface pairs are generated using different values of the Hurst exponent (0.6, 0.8, and 1.0), root mean square roughness RMS (50 μm, 100 μm, and 150 μm), and mismatch length (1000 μm, 2000 μm, and 3000 μm). The viscoelastic deformations of the seven surface pairs are simulated based on the Maxwell and Standard Linear Solid (SLS) models. For both models, the following conclusions can be drawn: (1)When increases, the average aperture decreases, and the contact ratio increases faster with time(2)When decreases, the average aperture decreases, and the contact ratio increases faster with time(3)When RMS decreases, the average aperture decreases, and the contact ratio increases faster with time(4)For the surface parameters ( between 0.6 and 1.0, RMS between 50 and 150 μm, and between 1000 and 3000 μm), macroscopic stresses (between 5 and 20 MPa), and the time durations ( for the Maxwell model and for the SLS model) used in the simulation, the contact ratio is generally less than 10%(5)The creep curves for different surface parameters can be normalized so the curves fall into a very narrow region
While the results are very informative, future work would be beneficial. Currently, the numerical code only simulates the deformation using the Maxwell and SLS models. Other viscoelastic models, such as the power law model and the Burgers model, could be implemented. In addition, the code only considers the viscoelastic deformation, while in reality, the contacting asperities may deform plastically. The code could be modified to simulate the viscoelastoplastic deformation of fracture surfaces. Last but not least, the code calculates the creep deformation under dry conditions. Pressure solution equations can be added in addition to the elastic deformation code to simulate the fracture creep deformation under wet conditions.
Data Availability
The data used to support the results of this research are available from the corresponding author upon request.
Conflicts of Interest
None of the authors have any conflicts of interest.
Acknowledgments
This work was partially funded by the Abu Dhabi National Oil Company (ADNOC) (project ID: 015RCM2015). The authors gratefully acknowledge Yves Bernabe, Brian Evans, and John Williams for their useful suggestions.
References
 L. F. Fan, G. W. Ma, and L. N. Y. Wong, “Effective viscoelastic behaviour of rock mass with doublescale discontinuities,” Geophysical Journal International, vol. 191, no. 1, pp. 147–154, 2012. View at: Publisher Site  Google Scholar
 H. Sone and M. D. Zoback, “Timedependent deformation of shale gas reservoir rocks and its longterm effect on the in situ state of stress,” International Journal of Rock Mechanics and Mining Sciences, vol. 69, pp. 120–132, 2014. View at: Publisher Site  Google Scholar
 Y. Zhao, Y. Wang, W. Wang et al., “Modeling of nonlinear rheological behavior of hard rock using triaxial rheological experiment,” International Journal of Rock Mechanics and Mining Sciences, vol. 2017, no. 93, pp. 66–75, 2017. View at: Publisher Site  Google Scholar
 Y. Zhao, L. Zhang, W. Wang, W. Wan, and W. Ma, “Separation of elastoviscoplastic strains of rock and a nonlinear creep model,” International Journal of Geomechanics, vol. 18, no. 1, article 04017129, 2018. View at: Publisher Site  Google Scholar
 H. Deng, J. P. Fitts, D. Crandall, D. McIntyre, and C. A. Peters, “Alterations of fractures in carbonate rocks by CO_{2}acidified brines,” Environmental Science and Technology, vol. 49, no. 16, pp. 10226–10234, 2015. View at: Publisher Site  Google Scholar
 R. L. Detwiler, “Experimental observations of deformation caused by mineral dissolution in variableaperture fractures,” Journal of Geophysical Research, vol. 113, no. B8, p. B08202, 2008. View at: Publisher Site  Google Scholar
 W. B. Durham, W. L. Bourcier, and E. A. Burton, “Direct observation of reactive flow in a single fracture,” Water Resources Research, vol. 37, no. 1, pp. 1–12, 2001. View at: Publisher Site  Google Scholar
 J. E. Elkhoury, P. Ameli, and R. L. Detwiler, “Dissolution and deformation in fractured carbonates caused by flow of CO_{2}rich brine under reservoir conditions,” International Journal of Greenhouse Gas Control, vol. 16, pp. S203–S215, 2013. View at: Publisher Site  Google Scholar
 M. GarciaRios, L. Luquot, J. M. Soler, and J. Cama, “The role of mineral heterogeneity on the hydrogeochemical response of two fractured reservoir rocks in contact with dissolved CO_{2},” Applied Geochemistry, vol. 84, pp. 202–217, 2017. View at: Publisher Site  Google Scholar
 C. Noiriel, P. Gouze, and B. Made, “3D analysis of geometry and flow changes in a limestone fracture during dissolution,” Journal of Hydrology, vol. 486, pp. 211–223, 2013. View at: Publisher Site  Google Scholar
 A. Polak, D. Elsworth, J. Liu, and A. S. Grader, “Spontaneous switching of permeability changes in a limestone fracture with net dissolution,” Water Resources Research, vol. 40, no. 3, 2004. View at: Publisher Site  Google Scholar
 N. M. Beeler and S. H. Hickman, “Stressinduced, timedependent fracture closure at hydrothermal conditions,” Journal of Geophysical Research, vol. 109, no. B2, article B02211, 2004. View at: Publisher Site  Google Scholar
 T. Ishibashi, T. P. McGuire, N. Watanabe, N. Tsuchiya, and D. Elsworth, “Permeability evolution in carbonate fractures: competing roles of confining stress and fluid PH,” Water Resources Research, vol. 49, no. 5, pp. 2828–2842, 2013. View at: Publisher Site  Google Scholar
 T. P. McGuire, D. Elsworth, and Z. Karcz, “Experimental measurements of stress and chemical controls on the evolution of fracture permeability,” Transport in Porous Media, vol. 98, no. 1, pp. 15–34, 2013. View at: Publisher Site  Google Scholar
 H. Yasuhara, A. Polak, Y. Mitani, A. Grader, P. Halleck, and D. Elsworth, “Evolution of fracture permeability through fluid rock reaction under hydrothermal conditions,” Earth and Planetary Science Letters, vol. 244, no. 1–2, pp. 186–200, 2006. View at: Publisher Site  Google Scholar
 H. Yasuhara, N. Kinoshita, H. Ohfuji, M. Takahashi, K. Ito, and K. Kishida, “Longterm observation of permeability in sedimentary rocks under hightemperature and stress conditions and its interpretation mediated by microstructural investigations,” Water Resources Research, vol. 51, no. 7, 2015. View at: Publisher Site  Google Scholar
 Y. Bernabé and B. Evans, “Numerical modelling of pressure solution deformation at axisymmetric asperities under normal load,” Geological Society, London, Special Publications, vol. 284, no. 1, pp. 185–205, 2007. View at: Publisher Site  Google Scholar
 P. S. Lang, A. Paluszny, and R. W. Zimmerman, “Hydraulic sealing due to pressure solution contact zone growth in siliciclastic rock fractures,” Journal of Geophysical Research: Solid Earth, vol. 120, no. 6, pp. 4080–4101, 2015. View at: Publisher Site  Google Scholar
 B. Li, Z. Zhao, Y. Jiang, and L. Jing, “Contact mechanism of a rock fracture subjected to normal loading and its impact on fast closure behavior during initial stage of fluid flow experiment,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 39, no. 13, pp. 1431–1449, 2015. View at: Publisher Site  Google Scholar
 I. Neretnieks, “Stressmediated closing of fractures: impact of matrix diffusion,” Journal of Geophysical Research: Solid Earth, vol. 119, no. 5, pp. 4149–4163, 2014. View at: Publisher Site  Google Scholar
 H. Yasuhara and D. Elsworth, “Evolution of permeability in a natural fracture: significant role of pressure solution,” Journal of Geophysical Research, vol. 109, p. B03204, 2004. View at: Google Scholar
 Z. Zhao, L. Liu, I. Neretnieks, and L. Jing, “Solute transport in a single fracture with timedependent aperture due to chemically medicated changes,” International Journal of Rock Mechanics and Mining Sciences, vol. 66, pp. 69–75, 2014. View at: Publisher Site  Google Scholar
 L. Zou, B. Li, Y. Mo, and V. Cvetkovic, “A highresolution contact analysis of roughwalled crystalline rock fractures subject to normal stress,” Rock Mechanics and Rock Engineering, vol. 53, no. 5, pp. 2141–2155, 2020. View at: Publisher Site  Google Scholar
 H. Kang, J. T. Germain, and H. H. Einstein, Creep Behavior of Musandam Limestone Fractures under Confining Pressures Less Than 1 MPa, Rock Mechanics/Geomechanics Symposium, U.S, 53rd edition, 2019.
 M. Iwano and H. H. Einstein, “Laboratory experiments on geometric and hydromechanical characteristics of three fractures in granodiorite,” in Proceedings of the Eighth International Congress of the ISRM, Japan, 1995. View at: Google Scholar
 S. R. Brown, A. Caprihan, and R. Hardy, “Experimental observation of fluid flow channels in a single fracture,” Journal of Geophysical Research, vol. 103, no. B3, pp. 5125–5132, 1998. View at: Publisher Site  Google Scholar
 N. Watanabe, N. Hirano, and N. Tsuchiya, “Determination of aperture structure and fluid flow in a rock fracture by highresolution numerical modeling on the basis of a flowthrough experiment under confining pressure,” Water Resources Research, vol. 44, 2008. View at: Publisher Site  Google Scholar
 S. R. Brown, “Simple mathematical model of a rough fracture,” Journal of Geophysical Research, vol. 100, no. B4, pp. 5941–5952, 1995. View at: Publisher Site  Google Scholar
 S. R. Brown and C. H. Scholz, “Broad bandwidth study of the topography of natural rock surfaces,” Journal of Geophysical Research, vol. 90, no. B14, pp. 12575–12,582, 1985. View at: Publisher Site  Google Scholar
 S. R. Brown and C. H. Scholz, “Closure of random elastic surfaces in contact,” Journal of Geophysical Research, vol. 90, no. B7, pp. 5531–5545, 1985. View at: Publisher Site  Google Scholar
 B. N. J. Persson, O. Albohr, U. Tartaglino, A. I. Volokitin, and E. Tosatti, “On the nature of surface roughness with application to contact mechanics, sealing, rubber friction and adhesion,” Journal of Physics: Condensed Matter, vol. 17, no. 1, pp. R1–R62, 2005. View at: Publisher Site  Google Scholar
 I. A. Polonsky and L. M. Keer, “A numerical method for solving rough contact problems based on the multilevel multisummation and conjugate gradient techniques,” Wear, vol. 231, no. 2, pp. 206–219, 1999. View at: Publisher Site  Google Scholar
 A. E. H. Love, “Stress produced in a semiinfinite solid by pressure on part of the boundary,” Philosophical Transactions of Royal Society, vol. A228, no. 337, pp. 54–59, 1929. View at: Google Scholar
 D. J. Andrews, “On modeling closure of rough surfaces in contact,” Eos Trans, AGU, vol. 69, pp. 14261427, 1988. View at: Google Scholar
 H. M. Stanley and T. Kato, “An FFTbased method for rough surface contact,” Journal of Tribology, vol. 119, no. 3, pp. 481–485, 1997. View at: Publisher Site  Google Scholar
 S. B. Liu, Q. Wang, and G. Liu, “A versatile method of discrete convolution and FFT (DCFFT) for contact analyses,” Wear, vol. 243, no. 12, pp. 101–111, 2000. View at: Publisher Site  Google Scholar
 W. Wayne Chen, Q. Jane Wang, Z. Huan, and X. Luo, “Semianalytical viscoelastic contact modeling of polymerbased materials,” Journal of Tribology, vol. 133, no. 4, article 041404, 2011. View at: Publisher Site  Google Scholar
 S. Spinu and D. Cerlinca, “Modelling of rough contact between linear viscoelastic materials,” Modelling and Simulation in Engineering, vol. 2017, Article ID 2521903, 11 pages, 2017. View at: Publisher Site  Google Scholar
 E. H. Lee and J. R. M. Radok, “The contact problem for viscoelastic bodies,” Journal of Applied Mechanics, vol. 27, no. 3, pp. 438–444, 1960. View at: Publisher Site  Google Scholar
 K. L. Johnson, Contact Mechanics, Cambridge University Press, London, 2012.
 S. Mighani, Y. Bernabe, A. Boulenouar, U. Mok, and B. Evans, “Creep deformation in Vaca Muerta Shale from nanoindentation to triaxial experiments,” Journal of Geophysical Research: Solid Earth, vol. 124, no. 8, pp. 7842–7868, 2019. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2020 Hao Kang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.