Geofluids
Volume 2020 (2020), Article ID 8879890, 30 pages
https://doi.org/10.1155/2020/8879890

Research Article

Numerical Simulation for Rock Fracture Viscoelastic Creep under Dry Conditions

Hao Kang,¹ Herbert Einstein,¹ Stephen Brown,¹ and John Germaine²

¹Massachusetts Institute of Technology, Cambridge, MA, USA
²Tufts University, Medford, MA, USA

Correspondence should be addressed to Hao Kang; kangh@mit.edu

Received 1 October 2020; Revised 17 November 2020; Accepted 23 November 2020; Published 10 December 2020

Academic Editor: Yanlin Zhao

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.

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 in-house 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 time-dependent 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 fracture-scale 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 above-mentioned 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 self-affine 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 self-affine fractal distribution ([29, 31]). A self-affine 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 self-affine fracture surfaces, the power spectrum can be related to the wavenumber ():

Figure 1: Schematic of surface profile breakdown. Each component is characterized by its wavelength (), amplitude (), and phase (relative position of the first peak of the sinusoid to that of all others).

Figure 2: Schematic of the power spectrum. The horizontal axis is the frequency (), and the vertical axis is the power (). The unit for is (1/length).

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)

Figure 3: Schematic of mismatch length scale, . (a) Degree of correlation with respect to the mismatch length, . (b) Example of wave components with different wavelengths. , , and represent three different wavelengths. For the left wave components, , and the two components match perfectly. For the right wave components, , and the phases of the two components do not match (they are independent). represents the mismatch in phase between the two components.

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).

Figure 4: Schematic of and . represents the height of the upper surface, represents the height of the lower surface, and represents the aperture (the gap between the two surfaces).

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 in-house 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 rigid-body 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.

Figure 5: Illustration of the Boussinesq solution. A normal stress field p is applied surrounding the location , and represents the direction elastic displacement at point .

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.

Figure 6: Illustration of the grid discretization and discretized Boussinesq equation. Indices and represent the coordinates in the direction, with and . Indices and represent the coordinates in the direction, with and . In Equation (13), the vertical pressure is applied on grid , and represents the direction elastic displacement on grid . and are the grid sizes in the direction and the direction, 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⁻¹ 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 in-house 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 time-dependent stress response to a step-change in strain, and (unit: 1/pressure) describes the time-dependent strain response to a step-change 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)

Figure 7: (a) Schematic for creep compliance under a constant stress . represents the creep compliance at the moment when is applied, and is equivalent to the elastic modulus. (b) Schematic for relaxation modulus under a constant strain . represents the relaxation modulus at the moment when is applied, and is equivalent to the elastic modulus.

The analytical solution for the Boussinesq and Cerrutti equation can be modified by adding time-dependent 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 left-hand 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].

Figure 8: Flowchart of the main algorithm.

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 half-spaces indented by a spherical rigid body. Lee and Radok [39] developed the analytical solutions for the two above-mentioned 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)

Figure 9: (a) Schematic of the Maxwell model. (b) Schematic of the SLS model.

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.

Figure 10: Illustration of the physical parameters in Equations (30) and (31).

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.

Figure 11: Comparison between analytical solutions ([39]; Johnson, 1985) and the numerical solutions for the Maxwell model. represents the maximum pressure when the time is zero (no creep, the solution is Hertzian contact solution), and represents the radius of the contacting region when the time is zero.

Figure 12: Comparison between analytical solutions ([39]; Johnson, 1985) and the numerical solutions for the SLS model. represents the maximum pressure when the time is zero (no creep, the solution is Hertzian contact solution), and represents the radius of the contacting region when the time is zero.

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.

Table 1: Natural rock joint profiles analysis results (obtained from [28]).

Rock joint specimen
Profile length (mm)
Hurst exponent ()
Normalized RMS ()
Normalized ()

Tuff
13.0

0.00752
0.014-0.092

Tuff
13.0

0.00525
0.039-0.205

Granodiorite
13.0

0.01708
0.014-0.041

Granodiorite
13.0

0.00854
0.010-0.062

Tuff
52.0

0.01458
0.063-0.143

Tuff
52.0

0.01350
0.019-0.072

Tuff
52.0

0.01850
0.014-0.053

Basalt
13.0

0.01792
0.022-0.078

Basalt
13.0

0.01085
0.017-0.049

Sandstone
13.0

0.00654
0.069-0.317

Chalk
26.0

0.00400
0.007-0.015

Tuff
26.0

0.00727
0.011-0.029

Tuff
26.0

0.00565
0.018-0.092

Siltstone
13.0

0.00349
0.059-0.205

Granite
13.0

0.01546
0.100-0.461

Granite
13.0

0.01446
0.053-0.270

Granite
13.0

0.01569
0.055-1.000

Rhyolitic dike
13.0

0.00758
0.021-0.048

Rhyolitic dike
13.0

0.01569
0.020-0.076

Rhyolitic dike
13.0

0.00665
0.014-0.032

Rhyolitic dike
13.0

0.00514
0.139-1.000

Rhyolitic dike
13.0

0.00688
0.029-0.084

Metasediment
13.0

0.00654
0.050-0.125

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: Summary of the seven synthetic surfaces.

Surface pair no.
Profile length (mm)

RMS

Absolute value (μm)

Absolute value (μm)

1
10
0.6
0.005
50
0.1
1000

2
10
0.8
0.005
50
0.1
1000

3
10
1.0
0.005
50
0.1
1000

4
10
0.8
0.005
50
0.2
2000

5
10
0.8
0.005
50
0.3
3000

6
10
0.8
0.010
100
0.1
1000

7
10
0.8
0.015
150
0.1
1000

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)

Figure 13: Effect of the value on the aperture field. and . The scale of the color bar is between 0 and 130 μm. The number of grids in the and directions is 512. (a) Aperture field when , which corresponds to surface pair 1 in Table 2; (b) aperture field when , which corresponds to surface pair 2 in Table 2; (c) aperture field when , which corresponds to surface pair 3 in Table 2. The value is increased from (a) to (c).

(a)

(b)

(c)

Figure 14: Effect of the value on the aperture field. and . The scale of the color bar is between 0 and 150 μm. The number of grids in the and directions is 512. (a) Aperture field when , which corresponds to surface pair 2 in Table 2; (b) aperture field when , which corresponds to surface pair 4 in Table 2; (c) aperture field when , which corresponds to surface pair 5 in Table 2. The value is increased from (a) to (c).

(a)

(b)

(c)

Figure 15: Effect of the RMS value on the aperture field. and . The scale of the color bar is proportional to the RMS value. The number of grids in the and directions is 512. (a) Aperture field when , which corresponds to surface pair 2 in Table 2; (b) aperture field when , which corresponds to surface pair 6 in Table 2; (c) aperture field when , which corresponds to surface pair 7 in Table 2. The RMS value is increased from (a) to (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.

Table 3: Mean and standard deviation of the aperture field for the seven surface pairs.

Surface pair no.

RMS (μm)
(μm)
Mean aperture (μm)
Standard deviation of aperture (μm)

1
0.6
50
1000
63.41
14.29

2
0.8
50
1000
37.30
8.57

3
1.0
50
1000
21.89
5.14

4
0.8
50
2000
55.94
15.01

5
0.8
50
3000
66.10
20.12

6
0.8
100
1000
74.59
17.15

7
0.8
150
1000
111.89
25.72

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].

Table 4: Input parameters for the Maxwell model (see Section 3 for the model explanation).

Parameters
Value

Shear modulus, (GPa)
7.0

Poisson’s ratio,
0.25

Viscosity, ()

Relaxation time, (s)

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

Figure 16: Creep deformation of the seven surface pairs (shown in Table 3). The time is normalized by . (a) Mean aperture changing with time for the seven surface pairs; (b) contact ratio changing with time for the seven surface pairs.

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.

Table 5: Conclusions for the Maxwell model simulation.

Parameters
Average aperture
Contact ratio

Initial value
Decrease rate
Initial value
Increase rate

↑
↓
↓
↑
↑

↓
↓
↓
↑
↑

RMS↓
↓
↓
↑
↑

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)

Figure 17: (a) Contact area and local contact stress before creep (the moment after the elastic deformation); (b) contact area and local contact stress after creep (the time duration is ). The number of grids in the and directions is 512.

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].

Table 6: Input parameters for the SLS model (see Section 3 for the model explanation).

Parameters
Value

Shear modulus, (GPa)
7.0

Shear modulus, (GPa)
7.0

Poisson’s ratio,
0.25

Viscosity, ()

Relaxation time, (s)

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

Figure 18: Creep deformation of the seven surface pairs (shown in Table 3). The time is normalized by . (a) Mean aperture changing with time for the seven surface pairs; (b) contact ratio changing with time for the seven surface pairs.

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.

Table 7: Conclusions for the SLS model simulation.

Parameters
Average aperture
Contact ratio

Initial value
Decrease rate
Initial value
Increase rate

↑
↓
↓
↑
↑

↓
↓
↓
↑
↑

RMS↓
↓
↓
↑
↑

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)

Figure 19: (a) Contact area and local contact stress before creep (the moment after the elastic deformation); (b) contact area and local contact stress after creep (the time duration is ). The number of grids in the and directions is 512.

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

Figure 20: Normalized mean aperture and contact ratio of the seven surface pairs (shown in Table 3). The time is normalized by . (a) Normalized mean aperture changing with time for the seven surface pairs; (b) normalized contact ratio changing with time for the seven surface pairs. The formulae for normalized mean aperture and normalized contact ratio are shown in Equations (37) and (38), respectively.

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

Figure 21: Normalized mean aperture and contact ratio of the seven surface pairs (shown in Table 3). The time is normalized by . (a) Normalized mean aperture changing with time for the seven surface pairs; (b) normalized contact ratio changing with time for the seven surface pairs. The formulae for normalized mean aperture and normalized contact ratio are shown in Equations (39) and (40), respectively.

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).

Figure 22: Histogram of local contact stress before and after the creep stage. The SLS model is used as the input, , , , , and the time duration is .

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 in-house 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 visco-elasto-plastic 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: 015-RCM-2015). 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 double-scale discontinuities,” Geophysical Journal International, vol. 191, no. 1, pp. 147–154, 2012.

H. Sone and M. D. Zoback, “Time-dependent deformation of shale gas reservoir rocks and its long-term effect on the in situ state of stress,” International Journal of Rock Mechanics and Mining Sciences, vol. 69, pp. 120–132, 2014.

Y. Zhao, Y. Wang, W. Wang et al., “Modeling of non-linear 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.

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.

H. Deng, J. P. Fitts, D. Crandall, D. McIntyre, and C. A. Peters, “Alterations of fractures in carbonate rocks by CO₂-acidified brines,” Environmental Science and Technology, vol. 49, no. 16, pp. 10226–10234, 2015.

R. L. Detwiler, “Experimental observations of deformation caused by mineral dissolution in variable-aperture fractures,” Journal of Geophysical Research, vol. 113, no. B8, p. B08202, 2008.

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.

J. E. Elkhoury, P. Ameli, and R. L. Detwiler, “Dissolution and deformation in fractured carbonates caused by flow of CO₂-rich brine under reservoir conditions,” International Journal of Greenhouse Gas Control, vol. 16, pp. S203–S215, 2013.

M. Garcia-Rios, 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₂,” Applied Geochemistry, vol. 84, pp. 202–217, 2017.

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.

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.

N. M. Beeler and S. H. Hickman, “Stress-induced, time-dependent fracture closure at hydrothermal conditions,” Journal of Geophysical Research, vol. 109, no. B2, article B02211, 2004.

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.

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.

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.

H. Yasuhara, N. Kinoshita, H. Ohfuji, M. Takahashi, K. Ito, and K. Kishida, “Long-term observation of permeability in sedimentary rocks under high-temperature and stress conditions and its interpretation mediated by microstructural investigations,” Water Resources Research, vol. 51, no. 7, 2015.

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.

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.

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.

I. Neretnieks, “Stress-mediated closing of fractures: impact of matrix diffusion,” Journal of Geophysical Research: Solid Earth, vol. 119, no. 5, pp. 4149–4163, 2014.

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.

Z. Zhao, L. Liu, I. Neretnieks, and L. Jing, “Solute transport in a single fracture with time-dependent aperture due to chemically medicated changes,” International Journal of Rock Mechanics and Mining Sciences, vol. 66, pp. 69–75, 2014.

L. Zou, B. Li, Y. Mo, and V. Cvetkovic, “A high-resolution contact analysis of rough-walled crystalline rock fractures subject to normal stress,” Rock Mechanics and Rock Engineering, vol. 53, no. 5, pp. 2141–2155, 2020.

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.

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.

N. Watanabe, N. Hirano, and N. Tsuchiya, “Determination of aperture structure and fluid flow in a rock fracture by high-resolution numerical modeling on the basis of a flow-through experiment under confining pressure,” Water Resources Research, vol. 44, 2008.

S. R. Brown, “Simple mathematical model of a rough fracture,” Journal of Geophysical Research, vol. 100, no. B4, pp. 5941–5952, 1995.

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.

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.

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.

I. A. Polonsky and L. M. Keer, “A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques,” Wear, vol. 231, no. 2, pp. 206–219, 1999.

A. E. H. Love, “Stress produced in a semi-infinite solid by pressure on part of the boundary,” Philosophical Transactions of Royal Society, vol. A228, no. 337, pp. 54–59, 1929.

D. J. Andrews, “On modeling closure of rough surfaces in contact,” Eos Trans, AGU, vol. 69, pp. 1426-1427, 1988.

H. M. Stanley and T. Kato, “An FFT-based method for rough surface contact,” Journal of Tribology, vol. 119, no. 3, pp. 481–485, 1997.

S. B. Liu, Q. Wang, and G. Liu, “A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses,” Wear, vol. 243, no. 1-2, pp. 101–111, 2000.

W. Wayne Chen, Q. Jane Wang, Z. Huan, and X. Luo, “Semi-analytical viscoelastic contact modeling of polymer-based materials,” Journal of Tribology, vol. 133, no. 4, article 041404, 2011.

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.

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.

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.