# Numerical Simulation for Rock Fracture Viscoelastic Creep under Dry Conditions

Authors: Hao Kang; Herbert Einstein; Stephen Brown; John Germaine
Journal: Geofluids
 (2020)

---

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


---

## Body

## 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 Sections2.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 (H), root mean square roughness at a reference length scale (σ), and mismatch length scale (λmis) [28].The surface roughness of rock fractures typically follows a self-affine fractal distribution ([29, 31]). A self-affine surface can be described as
(1)zx~ε−Hzεx,where z is the surface elevation, ε is a constant for x direction scaling, and H is the Hurst exponent. The H value is between 0 and 1; and a larger H value corresponds to a smoother local surface profile.H 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 (A2) is defined as the power, and a plot of power versus the inverse of wavelength (q=2π/λ) is defined as the power spectrum (see Figure 2 for a schematic). For self-affine fracture surfaces, the power spectrum Cq can be related to the wavenumber (q):
(2)Cq~q−21+H.Figure 1
Schematic of surface profile breakdown. Each component is characterized by its wavelength (λ), amplitude (A), 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 frequencyq (1/λ), and the vertical axis is the power C (A2). The unit for q is (1/length).The upper frequency end,qmax, is defined as qmax=2π/λ1, where λ1 is the smallest representable or measurable size. The lower frequency end, qmin, is defined as qmin=2π/λL, where λL is the surface length.The second parameter is the mismatch lengthλmis. As shown in Figure 1, each waveform can be completely described by three parameters: wavelength (λ), amplitude (A), 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 λmis, 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 C and q, but independent phases. Figure 3 illustrates the concept of the mismatch length scale.Figure 3
Schematic of mismatch length scale,λmis. (a) Degree of correlation with respect to the mismatch length, λmis. (b) Example of wave components with different wavelengths. λ1, λ2, and λmis represent three different wavelengths. For the left wave components, λ1>λmis, and the two components match perfectly. For the right wave components, λ2<λmis, and the phases of the two components do not match (they are independent). φmis represents the mismatch in phase between the two components.
(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
(3)σ2=∫qminqmaxCqdq.After generating the surface profile, its root mean square value,σini, is calculated. Then, the surface height is scaled linearly as
(4)zfinal=ziniσiniσfinal,where zini is the initial surface height, zfinal is the final surface height, and σfinal is the designated root mean square roughness. Figure 4 illustrates the concepts of z (surface height) and h (aperture).Figure 4
Schematic ofz and h. zupper represents the height of the upper surface, zlower represents the height of the lower surface, and h 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 functionhx,y is defined as
(5)hx,y=h0x,y+uex,y−δ,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), h0x,y is the initial aperture, and uex,y 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:(6)px,y>0andhx,y=0,px,y=0andhx,y>0,where px,y is the contact pressure (in the normal direction) acting on location x,y. 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 displacementuex,y can be calculated from the Boussinesq and Cerrutti (1885) equation:
(7)uex,y=∫−∞+∞∫−∞+∞Kx,y,x′,y′px′,y′dx′dy′,where Kx,y,x′,y′ is the influence matrix for normal displacement at surface location Bx,y due to a normal unit stress applied at location Ax′,y′ and px′,y′ is the normal stress applied at location Ax′,y′. 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 locationAx′,y′, and uex,y represents the z direction elastic displacement at point Bx,y.The influence matrixKx,y,x′,y′ is defined as
(8)Kx,y,x′,y′=1−υ2πGx−x′2+y−y′2,where G 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 Figure6, the surface is discretized into a regular grid:
(9)xi=i∆x,i=1,2,⋯,N,yj=j∆y,j=1,2,⋯,M,where xi and yj are x and y coordinates, respectively; Δx and Δy are the grid spacings in the x and y directions, respectively; and N and M are the total number of grids in the x and y directions, respectively.Figure 6
Illustration of the grid discretization and discretized Boussinesq equation. Indicesi and k represent the coordinates in the x direction, with 1≤i≤N and 1≤k≤N. Indices j and l represent the coordinates in the y direction, with 1≤j≤M and 1≤l≤M. In Equation (13), the vertical pressure pk,l is applied on grid k,l, and uei,j represents the z direction elastic displacement on grid i,j. Δx and Δy are the grid sizes in the x direction and the y direction, respectively.The aperture field can be calculated as(10)hi,j=h0i,j+uei,j−δ.The boundary conditions become(11)pi,j>0andhi,j=0,(12)pi,j=0andhi,j>0.Love [33] discretized the Boussinesq and Cerrutti equation as
(13)uei,j=∑k=1N∑l=1Mpk,l×Ki,k,j,l,(14)Ki,k,j,l=1−ν2πGalnc+a2+c2d+a2+d2+blnd+b2+d2c+b2+c2+clna+a2+c2b+b2+c2+dlnb+b2+d2a+a2+d2,where
(15)a=xi−xk+∆x2,b=xi−xk−∆x2,c=yj−yl+∆y2,d=yj−yl−∆y2.Figure6 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, uei,j can be expressed as
(16)uei,j=FFT−1FFTKi,k,j,l×FFTpk,l,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 N2×M2 operations, while by using the FFT method, calculating Equation (13) requires N×M×lnN×M operations. Therefore, when N and M are large, using the FFT method can save computation time.In addition, force balance needs to be achieved:(17)F=∑k=1N∑l=1Mpk,l,where F is the applied external force on the two fracture surfaces and pk,l represents the pressure acting on grid k,l (see Equation (12) and Figure 5). The indices k and l 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(18)εt=∫0tJt−τdστdtdτ,σt=∫0tErt−τdετdtdτ,where Ert and Jt are the relaxation modulus function and creep compliance function, respectively. Ert (unit: pressure) describes the time-dependent stress response to a step-change in strain, and Jt (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 0≤τ≤t (this means that t is a constant and τ is the variable). Figure 7 explains the creep compliance and relaxation modulus.Figure 7
(a) Schematic for creep compliance under a constant stressσ0. J0+ represents the creep compliance at the moment when σ0 is applied, and J0+ is equivalent to the elastic modulus. (b) Schematic for relaxation modulus under a constant strain ε0. Er0+ represents the relaxation modulus at the moment when ε0 is applied, and Er0+ is equivalent to the elastic modulus.
(a)(b)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
(19)uex,y,t=∫0tJt−τ∂∂τ∫−∞+∞∫−∞+∞px′,y′,τ1−νπx−x′2+y−y′2dx′dy′dτ.In Equation (19), the shear modulus 1/2G (in Equation (8)) is replaced by the creep compliance Jt−τ.Rearranging Equation (19) yields
(20)uex,y,t=∫0t∫−∞+∞∫−∞+∞Jt−τ1−νπx−x′2+y−y′2∂px′,y′,τ∂τdx′dy′dτ.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 (∫0t∂px′,y′,τ/∂τdτ) is added in Equation (23). Therefore, the time duration also needs to be discretized into many small time steps. Here, the total time duration t is discretized into Nt time steps with a uniform interval Δt. It is assumed that Δt is sufficiently small that the pressure field at each time step is unchanged. In addition, according to the fundamental theorem of calculus, the term ∂px′,y′,τdτ/∂τ in Equation (20) can be expressed as px′,y′,τ+dτ−px′,y′,τ. Based on the sufficiently small time step assumption and the fundamental theorem of calculus, Equation (20) can be rewritten as a summation:
(21)uex,y,α∆t=∑α′=1α∫−∞+∞∫−∞+∞Jα∆t−α′∆t1−νπx−x′2+y−y′2×px′,y′,α′−px′,y′,α′−1dx′dy′,where α=1,2,⋯,Nt.Within each time duration, in a given surface, the pressure field does not change with time. Therefore,px′,y′,α′−px′,y′,α′−1 can be factored out from the integral. Equation (21) becomes
(22)uex,y,α∆t=∑α′=1αpx′,y′,α′−px′,y′,α′−1∫−∞+∞∫−∞+∞Jα∆t−α′∆t1−νπx−x′2+y−y′2dx′dy′.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:
(23)∫−∞+∞∫−∞+∞1−ν2πGx−x′2+y−y′2dx′dy′⟶Discretize∑k=1N∑l=1MKi,k,j,l.By multiplying both sides of Equation (23) by 2GJα−α′∆t, 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
(24)∫−∞+∞∫−∞+∞Jα∆t−α′∆t1−νπx−x′2+y−y′2dx′dy′⟶Discretize∑k=1N∑l=1M2GJα−α′∆tKi,k,j,l,where G is the shear modulus.Thus, Equation (22) can be rewritten as
(25)uei,j,α∆t=∑α′=1Nt∑k=1N∑l=1M2GJα−α′∆tKi,k,j,lpk,l,α′−pk,l,α′−1.Equation (25) can be decomposed into Equations (26) and (27):
(26)uei,j,α∆t=∑α′=1Ntueα′,(27)ueα′=∑k=1N∑l=1M2GJα−α′∆tKi,k,j,lpk,l,α′−pk,l,α′−1.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):
(28)ueα′=FFT−1FFT2GJα−α′∆tKi,k,j,l×FFTpk,l,α′−pk,l,α′−1,where ueα′, 2GJα−α′∆tKi,j,k,l, and pk,l,α′−pk,l,α′−1 in Equation (28) are equivalent to uei,j, Ki,k,j,l, and pk,l 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 historyFt, initial discretized aperture field hi,j, and the creep compliance function Jt for the viscoelastic material(2)
CalculateJα−α′Δt for each time step and the stiffness matrix Ki,k,j,l(3)
Use the method described in Section2.2; calculate the initial contact state and the pressure distribution (α=1). 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 pk,l,1,⋯,pk,l,α−1 are known, the pressure distribution field and displacement field at time step α can be calculated by the CG method(6)
Ifα reaches Nt (the total number of time steps), stop the program, and export the pressure and displacement field. Otherwise, go back to step 4Figure8 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. Figures9(a) and 9(b) show the schematic for the Maxwell model and SLS model, respectively.Figure 9
(a) Schematic of the Maxwell model. (b) Schematic of the SLS model.
(a)(b)The Maxwell model consists of a spring and a dashpot. The spring represents elasticity, with a shear modulus ofG; and the dashpot represents viscosity, with a viscosity of η (unit: time ∗ pressure). Under constant stress, the strain can be expressed as
(29)εt=σ01G+tη.Equation (29) indicates that under constant stress, the strain rate (creep rate) is constant.The creep compliance function can therefore be defined as(30)Jt=1G+tη,where t is the time. The relaxation time for the Maxwell model, T, is defined as
(31)T=ηG.The SLS model consists of two springs and one dashpot. The two springs represent elasticity, with an elastic modulus ofG1 and G2, respectively; and the dashpot represents viscosity, with a viscosity of η (unit: time ∗ pressure). Under constant stress σ0, the strain εt can be related to σ0 as
(32)σt=G1G2G1+G2εt+G1ηG1+G2dεtdt.The creep compliance function can be expressed as(33)Jt=1G1+1−e−tG2/ηG2,where t is the time. The relaxation time for the SLS model, T, is defined as
(34)T=ηG2.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
(35)pt,r=2πR1−υ∫0tGe−t−t′G/ηddt′a2t′−r21/2dt′,at=31−νRP41G+tη3,where p is the pressure, r is the distance from the center of the contacting region, t is the total time duration, at is the radius of the contacting region, R is the spherical diameter, G is the shear modulus, η is the viscosity, and P 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
(36)pt,r=2πR1−υ∫0tG1G1+G2G2+G1e−t−t′G1+G2/ηddt′a2t′−r21/2dt′,at=31−νRP41G1+1G21−e−tG2/η3,where p is the pressure, r is the distance from the center of the contacting region, t is the total time duration, at is the radius of the contacting region, R is the spherical diameter, G1 and G2 are the shear moduli, η is the viscosity, and P 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 P and R do not matter, because the pressure distribution p and the radius of the contacting region r are normalized.Figure 11
Comparison between analytical solutions ([39]; Johnson, 1985) and the numerical solutions for the Maxwell model. ph represents the maximum pressure when the time is zero (no creep, the solution is Hertzian contact solution), and rh 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. ph represents the maximum pressure when the time is zero (no creep, the solution is Hertzian contact solution), and rh represents the radius of the contacting region when the time is zero.The results in Figures11 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.


## 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 (H), root mean square roughness at a reference length scale (σ), and mismatch length scale (λmis) [28].The surface roughness of rock fractures typically follows a self-affine fractal distribution ([29, 31]). A self-affine surface can be described as
(1)zx~ε−Hzεx,where z is the surface elevation, ε is a constant for x direction scaling, and H is the Hurst exponent. The H value is between 0 and 1; and a larger H value corresponds to a smoother local surface profile.H 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 (A2) is defined as the power, and a plot of power versus the inverse of wavelength (q=2π/λ) is defined as the power spectrum (see Figure 2 for a schematic). For self-affine fracture surfaces, the power spectrum Cq can be related to the wavenumber (q):
(2)Cq~q−21+H.Figure 1
Schematic of surface profile breakdown. Each component is characterized by its wavelength (λ), amplitude (A), 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 frequencyq (1/λ), and the vertical axis is the power C (A2). The unit for q is (1/length).The upper frequency end,qmax, is defined as qmax=2π/λ1, where λ1 is the smallest representable or measurable size. The lower frequency end, qmin, is defined as qmin=2π/λL, where λL is the surface length.The second parameter is the mismatch lengthλmis. As shown in Figure 1, each waveform can be completely described by three parameters: wavelength (λ), amplitude (A), 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 λmis, 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 C and q, but independent phases. Figure 3 illustrates the concept of the mismatch length scale.Figure 3
Schematic of mismatch length scale,λmis. (a) Degree of correlation with respect to the mismatch length, λmis. (b) Example of wave components with different wavelengths. λ1, λ2, and λmis represent three different wavelengths. For the left wave components, λ1>λmis, and the two components match perfectly. For the right wave components, λ2<λmis, and the phases of the two components do not match (they are independent). φmis represents the mismatch in phase between the two components.
(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
(3)σ2=∫qminqmaxCqdq.After generating the surface profile, its root mean square value,σini, is calculated. Then, the surface height is scaled linearly as
(4)zfinal=ziniσiniσfinal,where zini is the initial surface height, zfinal is the final surface height, and σfinal is the designated root mean square roughness. Figure 4 illustrates the concepts of z (surface height) and h (aperture).Figure 4
Schematic ofz and h. zupper represents the height of the upper surface, zlower represents the height of the lower surface, and h 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 functionhx,y is defined as
(5)hx,y=h0x,y+uex,y−δ,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), h0x,y is the initial aperture, and uex,y 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:(6)px,y>0andhx,y=0,px,y=0andhx,y>0,where px,y is the contact pressure (in the normal direction) acting on location x,y. 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 displacementuex,y can be calculated from the Boussinesq and Cerrutti (1885) equation:
(7)uex,y=∫−∞+∞∫−∞+∞Kx,y,x′,y′px′,y′dx′dy′,where Kx,y,x′,y′ is the influence matrix for normal displacement at surface location Bx,y due to a normal unit stress applied at location Ax′,y′ and px′,y′ is the normal stress applied at location Ax′,y′. 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 locationAx′,y′, and uex,y represents the z direction elastic displacement at point Bx,y.The influence matrixKx,y,x′,y′ is defined as
(8)Kx,y,x′,y′=1−υ2πGx−x′2+y−y′2,where G 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 Figure6, the surface is discretized into a regular grid:
(9)xi=i∆x,i=1,2,⋯,N,yj=j∆y,j=1,2,⋯,M,where xi and yj are x and y coordinates, respectively; Δx and Δy are the grid spacings in the x and y directions, respectively; and N and M are the total number of grids in the x and y directions, respectively.Figure 6
Illustration of the grid discretization and discretized Boussinesq equation. Indicesi and k represent the coordinates in the x direction, with 1≤i≤N and 1≤k≤N. Indices j and l represent the coordinates in the y direction, with 1≤j≤M and 1≤l≤M. In Equation (13), the vertical pressure pk,l is applied on grid k,l, and uei,j represents the z direction elastic displacement on grid i,j. Δx and Δy are the grid sizes in the x direction and the y direction, respectively.The aperture field can be calculated as(10)hi,j=h0i,j+uei,j−δ.The boundary conditions become(11)pi,j>0andhi,j=0,(12)pi,j=0andhi,j>0.Love [33] discretized the Boussinesq and Cerrutti equation as
(13)uei,j=∑k=1N∑l=1Mpk,l×Ki,k,j,l,(14)Ki,k,j,l=1−ν2πGalnc+a2+c2d+a2+d2+blnd+b2+d2c+b2+c2+clna+a2+c2b+b2+c2+dlnb+b2+d2a+a2+d2,where
(15)a=xi−xk+∆x2,b=xi−xk−∆x2,c=yj−yl+∆y2,d=yj−yl−∆y2.Figure6 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, uei,j can be expressed as
(16)uei,j=FFT−1FFTKi,k,j,l×FFTpk,l,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 N2×M2 operations, while by using the FFT method, calculating Equation (13) requires N×M×lnN×M operations. Therefore, when N and M are large, using the FFT method can save computation time.In addition, force balance needs to be achieved:(17)F=∑k=1N∑l=1Mpk,l,where F is the applied external force on the two fracture surfaces and pk,l represents the pressure acting on grid k,l (see Equation (12) and Figure 5). The indices k and l 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(18)εt=∫0tJt−τdστdtdτ,σt=∫0tErt−τdετdtdτ,where Ert and Jt are the relaxation modulus function and creep compliance function, respectively. Ert (unit: pressure) describes the time-dependent stress response to a step-change in strain, and Jt (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 0≤τ≤t (this means that t is a constant and τ is the variable). Figure 7 explains the creep compliance and relaxation modulus.Figure 7
(a) Schematic for creep compliance under a constant stressσ0. J0+ represents the creep compliance at the moment when σ0 is applied, and J0+ is equivalent to the elastic modulus. (b) Schematic for relaxation modulus under a constant strain ε0. Er0+ represents the relaxation modulus at the moment when ε0 is applied, and Er0+ is equivalent to the elastic modulus.
(a)(b)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
(19)uex,y,t=∫0tJt−τ∂∂τ∫−∞+∞∫−∞+∞px′,y′,τ1−νπx−x′2+y−y′2dx′dy′dτ.In Equation (19), the shear modulus 1/2G (in Equation (8)) is replaced by the creep compliance Jt−τ.Rearranging Equation (19) yields
(20)uex,y,t=∫0t∫−∞+∞∫−∞+∞Jt−τ1−νπx−x′2+y−y′2∂px′,y′,τ∂τdx′dy′dτ.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 (∫0t∂px′,y′,τ/∂τdτ) is added in Equation (23). Therefore, the time duration also needs to be discretized into many small time steps. Here, the total time duration t is discretized into Nt time steps with a uniform interval Δt. It is assumed that Δt is sufficiently small that the pressure field at each time step is unchanged. In addition, according to the fundamental theorem of calculus, the term ∂px′,y′,τdτ/∂τ in Equation (20) can be expressed as px′,y′,τ+dτ−px′,y′,τ. Based on the sufficiently small time step assumption and the fundamental theorem of calculus, Equation (20) can be rewritten as a summation:
(21)uex,y,α∆t=∑α′=1α∫−∞+∞∫−∞+∞Jα∆t−α′∆t1−νπx−x′2+y−y′2×px′,y′,α′−px′,y′,α′−1dx′dy′,where α=1,2,⋯,Nt.Within each time duration, in a given surface, the pressure field does not change with time. Therefore,px′,y′,α′−px′,y′,α′−1 can be factored out from the integral. Equation (21) becomes
(22)uex,y,α∆t=∑α′=1αpx′,y′,α′−px′,y′,α′−1∫−∞+∞∫−∞+∞Jα∆t−α′∆t1−νπx−x′2+y−y′2dx′dy′.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:
(23)∫−∞+∞∫−∞+∞1−ν2πGx−x′2+y−y′2dx′dy′⟶Discretize∑k=1N∑l=1MKi,k,j,l.By multiplying both sides of Equation (23) by 2GJα−α′∆t, 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
(24)∫−∞+∞∫−∞+∞Jα∆t−α′∆t1−νπx−x′2+y−y′2dx′dy′⟶Discretize∑k=1N∑l=1M2GJα−α′∆tKi,k,j,l,where G is the shear modulus.Thus, Equation (22) can be rewritten as
(25)uei,j,α∆t=∑α′=1Nt∑k=1N∑l=1M2GJα−α′∆tKi,k,j,lpk,l,α′−pk,l,α′−1.Equation (25) can be decomposed into Equations (26) and (27):
(26)uei,j,α∆t=∑α′=1Ntueα′,(27)ueα′=∑k=1N∑l=1M2GJα−α′∆tKi,k,j,lpk,l,α′−pk,l,α′−1.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):
(28)ueα′=FFT−1FFT2GJα−α′∆tKi,k,j,l×FFTpk,l,α′−pk,l,α′−1,where ueα′, 2GJα−α′∆tKi,j,k,l, and pk,l,α′−pk,l,α′−1 in Equation (28) are equivalent to uei,j, Ki,k,j,l, and pk,l 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 historyFt, initial discretized aperture field hi,j, and the creep compliance function Jt for the viscoelastic material(2)
CalculateJα−α′Δt for each time step and the stiffness matrix Ki,k,j,l(3)
Use the method described in Section2.2; calculate the initial contact state and the pressure distribution (α=1). 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 pk,l,1,⋯,pk,l,α−1 are known, the pressure distribution field and displacement field at time step α can be calculated by the CG method(6)
Ifα reaches Nt (the total number of time steps), stop the program, and export the pressure and displacement field. Otherwise, go back to step 4Figure8 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. Figures9(a) and 9(b) show the schematic for the Maxwell model and SLS model, respectively.Figure 9
(a) Schematic of the Maxwell model. (b) Schematic of the SLS model.
(a)(b)The Maxwell model consists of a spring and a dashpot. The spring represents elasticity, with a shear modulus ofG; and the dashpot represents viscosity, with a viscosity of η (unit: time ∗ pressure). Under constant stress, the strain can be expressed as
(29)εt=σ01G+tη.Equation (29) indicates that under constant stress, the strain rate (creep rate) is constant.The creep compliance function can therefore be defined as(30)Jt=1G+tη,where t is the time. The relaxation time for the Maxwell model, T, is defined as
(31)T=ηG.The SLS model consists of two springs and one dashpot. The two springs represent elasticity, with an elastic modulus ofG1 and G2, respectively; and the dashpot represents viscosity, with a viscosity of η (unit: time ∗ pressure). Under constant stress σ0, the strain εt can be related to σ0 as
(32)σt=G1G2G1+G2εt+G1ηG1+G2dεtdt.The creep compliance function can be expressed as(33)Jt=1G1+1−e−tG2/ηG2,where t is the time. The relaxation time for the SLS model, T, is defined as
(34)T=ηG2.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
(35)pt,r=2πR1−υ∫0tGe−t−t′G/ηddt′a2t′−r21/2dt′,at=31−νRP41G+tη3,where p is the pressure, r is the distance from the center of the contacting region, t is the total time duration, at is the radius of the contacting region, R is the spherical diameter, G is the shear modulus, η is the viscosity, and P 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
(36)pt,r=2πR1−υ∫0tG1G1+G2G2+G1e−t−t′G1+G2/ηddt′a2t′−r21/2dt′,at=31−νRP41G1+1G21−e−tG2/η3,where p is the pressure, r is the distance from the center of the contacting region, t is the total time duration, at is the radius of the contacting region, R is the spherical diameter, G1 and G2 are the shear moduli, η is the viscosity, and P 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 P and R do not matter, because the pressure distribution p and the radius of the contacting region r are normalized.Figure 11
Comparison between analytical solutions ([39]; Johnson, 1985) and the numerical solutions for the Maxwell model. ph represents the maximum pressure when the time is zero (no creep, the solution is Hertzian contact solution), and rh 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. ph represents the maximum pressure when the time is zero (no creep, the solution is Hertzian contact solution), and rh represents the radius of the contacting region when the time is zero.The results in Figures11 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 Section2.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 (H), root mean square roughness (RMS), and mismatch length (λc) of each rock joint. The rock types include tuff, granodiorite, granite, basalt, chalk, and siltstone. Table 1 summarizes his analysis results. RMS and λc are normalized by the profile length.Table 1
Natural rock joint profiles analysis results (obtained from [28]).
Rock joint specimenProfile lengthL (mm)Hurst exponent (H)Normalized RMS (RMS/L)Normalizedλc (λc/L)Tuff13.00.59±0.330.007520.014-0.092Tuff13.00.55±0.310.005250.039-0.205Granodiorite13.00.91±0.400.017080.014-0.041Granodiorite13.00.68±0.360.008540.010-0.062Tuff52.00.70±0.250.014580.063-0.143Tuff52.00.83±0.270.013500.019-0.072Tuff52.00.73±0.260.018500.014-0.053Basalt13.00.85±0.420.017920.022-0.078Basalt13.00.77±0.390.010850.017-0.049Sandstone13.00.48±0.320.006540.069-0.317Chalk26.00.59±0.230.004000.007-0.015Tuff26.00.75±0.320.007270.011-0.029Tuff26.00.77±0.280.005650.018-0.092Siltstone13.00.52±0.320.003490.059-0.205Granite13.00.49±0.300.015460.100-0.461Granite13.00.51±0.310.014460.053-0.270Granite13.00.79±0.380.015690.055-1.000Rhyolitic dike13.00.66±0.370.007580.021-0.048Rhyolitic dike13.00.60±0.360.015690.020-0.076Rhyolitic dike13.00.64±0.360.006650.014-0.032Rhyolitic dike13.00.66±0.360.005140.139-1.000Rhyolitic dike13.00.65±0.340.006880.029-0.084Metasediment13.00.64±0.340.006540.050-0.125Based on the results in Table1, it can be concluded that the Hurst exponent is normally between 0.50 and 0.90; the normalized RMS (RMS/L) is normally between 0.005 and 0.015; and the normalized λc (λc/L) 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 H, RMS, and λc values. Table 2 summarizes the information of the seven pairs of the synthetic surfaces. For all the surfaces, the dimensions (profile lengths) in the x and y directions are 10 mm.Table 2
Summary of the seven synthetic surfaces.
Surface pair no.Profile lengthL (mm)HRMSλcRMS/LAbsolute value (μm)λc/LAbsolute value (μm)1100.60.005500.110002100.80.005500.110003101.00.005500.110004100.80.005500.220005100.80.005500.330006100.80.0101000.110007100.80.0151500.11000Table2 shows that in surface pairs 1, 2, and 3, the H value is varied; in surface pairs 2, 4, and 5, the λc 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. Figures13–15 show the effect of different H, λc, and RMS values on the aperture field, respectively. The following conclusions can be drawn:
(1)
From Figure13, as H decreases, the mean and standard deviation of the aperture increase(2)
From Figure14, as λc increases, the mean and standard deviation of the aperture increase(3)
From Figure15, 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)Figure 13
Effect of theH value on the aperture field. RMS=50μm and λc=1000μm. The scale of the color bar is between 0 and 130 μm. The number of grids in the x and y directions is 512. (a) Aperture field when H=0.6, which corresponds to surface pair 1 in Table 2; (b) aperture field when H=0.8, which corresponds to surface pair 2 in Table 2; (c) aperture field when H=1.0, which corresponds to surface pair 3 in Table 2. The H value is increased from (a) to (c).
(a)(b)(c)Figure 14
Effect of theλc value on the aperture field. RMS=50μm and H=0.8. The scale of the color bar is between 0 and 150 μm. The number of grids in the x and y directions is 512. (a) Aperture field when λc=1000μm, which corresponds to surface pair 2 in Table 2; (b) aperture field when λc=2000μm, which corresponds to surface pair 4 in Table 2; (c) aperture field when λc=3000μm, which corresponds to surface pair 5 in Table 2. The λc value is increased from (a) to (c).
(a)(b)(c)Figure 15
Effect of the RMS value on the aperture field.H=0.8 and λc=1000μm. The scale of the color bar is proportional to the RMS value. The number of grids in the x and y directions is 512. (a) Aperture field when RMS=50μm, which corresponds to surface pair 2 in Table 2; (b) aperture field when RMS=100μm, which corresponds to surface pair 6 in Table 2; (c) aperture field when RMS=150μm, which corresponds to surface pair 7 in Table 2. The RMS value is increased from (a) to (c).
(a)(b)(c)Table3 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.HRMS (μm)λc (μm)Mean aperture (μm)Standard deviation of aperture (μm)10.650100063.4114.2920.850100037.308.5731.050100021.895.1440.850200055.9415.0150.850300066.1020.1260.8100100074.5917.1570.81501000111.8925.72


### 3.2. Creep Simulation Results for the Maxwell Model

The creep deformations of the seven surface pairs are calculated using the Maxwell model. Table4 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 Section3 for the model explanation).
ParametersValueShear modulus,G (GPa)7.0Poisson’s ratio,υ0.25Viscosity,η (GPa∗s)2.0×107Relaxation time,τ=η/G (s)2.857×106The following parameters are also defined:(1)
Macroscopic stressσ1=F/A, where F is the total force applied on the fracture surface and A is the area of the fracture surface (10mm∗10mm)(2)
Contactratio%=No.ofcontactingcells/No.oftotalcells∗100. See the definition in Figure 6Figure16 summarizes the effect of different H, λc, and RMS values on the creep deformation. In the simulation, σ1 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 2τ.Figure 16
Creep deformation of the seven surface pairs (shown in Table3). The time t 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.
(a)
Mean aperture(b)
Contact ratioThe following conclusions can be drawn:(1)
Based on the curves of surface pairs 1, 2, and 3, whenH 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λc 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%Table5 summarizes the above conclusions.Table 5
Conclusions for the Maxwell model simulation.
ParametersAverage apertureContact ratioInitial valueDecrease rateInitial valueIncrease rateH↑↓↓↑↑λc↓↓↓↑↑RMS↓↓↓↑↑Note: “↓” means decrease and “↑”means increase.Figure17 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 H=1.0, RMS=50μm, λc=1000μm, and σ=10MPa, and the creep timeduration=2τ. 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 2τ. 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).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 is2τ). The number of grids in the x and y directions is 512.
(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. Table6 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 Section3 for the model explanation).
ParametersValueShear modulus,G1 (GPa)7.0Shear modulus,G2 (GPa)7.0Poisson’s ratio,υ0.25Viscosity,η (GPa∗s)2.0×107Relaxation time,τ=η2/G2 (s)2.857×106Figure18 summarizes the effect of different H, λc, and RMS values on the creep deformation. In the simulation, σ1 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 5τ. Compared with the Maxwell model, the time duration is extended from 2τ to 5τ to better illustrate the decreasing creep rate.Figure 18
Creep deformation of the seven surface pairs (shown in Table3). The time t 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.
(a)
Mean aperture(b)
Contact ratioThe following conclusions can be drawn:(1)
Based on the curves of surface pairs 1, 2, and 3, whenH 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λc 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 rateTable7 summarizes the above conclusions. The effects of H, RMS, and λc on the creep rate under the SLS model are similar to those under the Maxwell model.Table 7
Conclusions for the SLS model simulation.
ParametersAverage apertureContact ratioInitial valueDecrease rateInitial valueIncrease rateH↑↓↓↑↑λc↓↓↓↑↑RMS↓↓↓↑↑Note: “↓” means decrease and “↑” means increase.Figure19 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 H=1.0, RMS=50μm, λc=1000μm, and σ=10MPa, and the creep time duration=5τ. 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 5τ. 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).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 is5τ). The number of grids in the x and y directions is 512.
(a)(b)


## 3.1. Surface Generation Results

As described in Section2.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 (H), root mean square roughness (RMS), and mismatch length (λc) of each rock joint. The rock types include tuff, granodiorite, granite, basalt, chalk, and siltstone. Table 1 summarizes his analysis results. RMS and λc are normalized by the profile length.Table 1
Natural rock joint profiles analysis results (obtained from [28]).
Rock joint specimenProfile lengthL (mm)Hurst exponent (H)Normalized RMS (RMS/L)Normalizedλc (λc/L)Tuff13.00.59±0.330.007520.014-0.092Tuff13.00.55±0.310.005250.039-0.205Granodiorite13.00.91±0.400.017080.014-0.041Granodiorite13.00.68±0.360.008540.010-0.062Tuff52.00.70±0.250.014580.063-0.143Tuff52.00.83±0.270.013500.019-0.072Tuff52.00.73±0.260.018500.014-0.053Basalt13.00.85±0.420.017920.022-0.078Basalt13.00.77±0.390.010850.017-0.049Sandstone13.00.48±0.320.006540.069-0.317Chalk26.00.59±0.230.004000.007-0.015Tuff26.00.75±0.320.007270.011-0.029Tuff26.00.77±0.280.005650.018-0.092Siltstone13.00.52±0.320.003490.059-0.205Granite13.00.49±0.300.015460.100-0.461Granite13.00.51±0.310.014460.053-0.270Granite13.00.79±0.380.015690.055-1.000Rhyolitic dike13.00.66±0.370.007580.021-0.048Rhyolitic dike13.00.60±0.360.015690.020-0.076Rhyolitic dike13.00.64±0.360.006650.014-0.032Rhyolitic dike13.00.66±0.360.005140.139-1.000Rhyolitic dike13.00.65±0.340.006880.029-0.084Metasediment13.00.64±0.340.006540.050-0.125Based on the results in Table1, it can be concluded that the Hurst exponent is normally between 0.50 and 0.90; the normalized RMS (RMS/L) is normally between 0.005 and 0.015; and the normalized λc (λc/L) 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 H, RMS, and λc values. Table 2 summarizes the information of the seven pairs of the synthetic surfaces. For all the surfaces, the dimensions (profile lengths) in the x and y directions are 10 mm.Table 2
Summary of the seven synthetic surfaces.
Surface pair no.Profile lengthL (mm)HRMSλcRMS/LAbsolute value (μm)λc/LAbsolute value (μm)1100.60.005500.110002100.80.005500.110003101.00.005500.110004100.80.005500.220005100.80.005500.330006100.80.0101000.110007100.80.0151500.11000Table2 shows that in surface pairs 1, 2, and 3, the H value is varied; in surface pairs 2, 4, and 5, the λc 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. Figures13–15 show the effect of different H, λc, and RMS values on the aperture field, respectively. The following conclusions can be drawn:
(1)
From Figure13, as H decreases, the mean and standard deviation of the aperture increase(2)
From Figure14, as λc increases, the mean and standard deviation of the aperture increase(3)
From Figure15, 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)Figure 13
Effect of theH value on the aperture field. RMS=50μm and λc=1000μm. The scale of the color bar is between 0 and 130 μm. The number of grids in the x and y directions is 512. (a) Aperture field when H=0.6, which corresponds to surface pair 1 in Table 2; (b) aperture field when H=0.8, which corresponds to surface pair 2 in Table 2; (c) aperture field when H=1.0, which corresponds to surface pair 3 in Table 2. The H value is increased from (a) to (c).
(a)(b)(c)Figure 14
Effect of theλc value on the aperture field. RMS=50μm and H=0.8. The scale of the color bar is between 0 and 150 μm. The number of grids in the x and y directions is 512. (a) Aperture field when λc=1000μm, which corresponds to surface pair 2 in Table 2; (b) aperture field when λc=2000μm, which corresponds to surface pair 4 in Table 2; (c) aperture field when λc=3000μm, which corresponds to surface pair 5 in Table 2. The λc value is increased from (a) to (c).
(a)(b)(c)Figure 15
Effect of the RMS value on the aperture field.H=0.8 and λc=1000μm. The scale of the color bar is proportional to the RMS value. The number of grids in the x and y directions is 512. (a) Aperture field when RMS=50μm, which corresponds to surface pair 2 in Table 2; (b) aperture field when RMS=100μm, which corresponds to surface pair 6 in Table 2; (c) aperture field when RMS=150μm, which corresponds to surface pair 7 in Table 2. The RMS value is increased from (a) to (c).
(a)(b)(c)Table3 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.HRMS (μm)λc (μm)Mean aperture (μm)Standard deviation of aperture (μm)10.650100063.4114.2920.850100037.308.5731.050100021.895.1440.850200055.9415.0150.850300066.1020.1260.8100100074.5917.1570.81501000111.8925.72


## 3.2. Creep Simulation Results for the Maxwell Model

The creep deformations of the seven surface pairs are calculated using the Maxwell model. Table4 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 Section3 for the model explanation).
ParametersValueShear modulus,G (GPa)7.0Poisson’s ratio,υ0.25Viscosity,η (GPa∗s)2.0×107Relaxation time,τ=η/G (s)2.857×106The following parameters are also defined:(1)
Macroscopic stressσ1=F/A, where F is the total force applied on the fracture surface and A is the area of the fracture surface (10mm∗10mm)(2)
Contactratio%=No.ofcontactingcells/No.oftotalcells∗100. See the definition in Figure 6Figure16 summarizes the effect of different H, λc, and RMS values on the creep deformation. In the simulation, σ1 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 2τ.Figure 16
Creep deformation of the seven surface pairs (shown in Table3). The time t 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.
(a)
Mean aperture(b)
Contact ratioThe following conclusions can be drawn:(1)
Based on the curves of surface pairs 1, 2, and 3, whenH 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λc 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%Table5 summarizes the above conclusions.Table 5
Conclusions for the Maxwell model simulation.
ParametersAverage apertureContact ratioInitial valueDecrease rateInitial valueIncrease rateH↑↓↓↑↑λc↓↓↓↑↑RMS↓↓↓↑↑Note: “↓” means decrease and “↑”means increase.Figure17 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 H=1.0, RMS=50μm, λc=1000μm, and σ=10MPa, and the creep timeduration=2τ. 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 2τ. 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).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 is2τ). The number of grids in the x and y directions is 512.
(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. Table6 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 Section3 for the model explanation).
ParametersValueShear modulus,G1 (GPa)7.0Shear modulus,G2 (GPa)7.0Poisson’s ratio,υ0.25Viscosity,η (GPa∗s)2.0×107Relaxation time,τ=η2/G2 (s)2.857×106Figure18 summarizes the effect of different H, λc, and RMS values on the creep deformation. In the simulation, σ1 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 5τ. Compared with the Maxwell model, the time duration is extended from 2τ to 5τ to better illustrate the decreasing creep rate.Figure 18
Creep deformation of the seven surface pairs (shown in Table3). The time t 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.
(a)
Mean aperture(b)
Contact ratioThe following conclusions can be drawn:(1)
Based on the curves of surface pairs 1, 2, and 3, whenH 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λc 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 rateTable7 summarizes the above conclusions. The effects of H, RMS, and λc on the creep rate under the SLS model are similar to those under the Maxwell model.Table 7
Conclusions for the SLS model simulation.
ParametersAverage apertureContact ratioInitial valueDecrease rateInitial valueIncrease rateH↑↓↓↑↑λc↓↓↓↑↑RMS↓↓↓↑↑Note: “↓” means decrease and “↑” means increase.Figure19 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 H=1.0, RMS=50μm, λc=1000μm, and σ=10MPa, and the creep time duration=5τ. 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 5τ. 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).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 is5τ). The number of grids in the x and y directions is 512.
(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 differentH, λc, 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 2τ, the mean aperture, h, and the contact ratio, c, are normalized as
(37)hnor=helashelas−refh,(38)cnor=celascelas−refc,where hnor is the normalized mean aperture, helas is the mean aperture right after the elastic deformation (before creep), helas−ref is the mean aperture right after the elastic deformation for surface pair 2 (H=0.8, RMS=50μm, and λc=1000μm), cnor is the normalized contact ratio, celas is the contact ratio right after the elastic deformation (before creep), and celas−ref is the contact ratio right after the elastic deformation for surface pair 2.Figure20 summarizes the normalized mean aperture and contact ratio for different H, λc, and RMS values. σ1 is fixed at 10 MPa, and the time duration is from 0 to 2τ.Figure 20
Normalized mean aperture and contact ratio of the seven surface pairs (shown in Table3). The time t 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.
(a)
Normalized mean aperture(b)
Normalized contact ratioThe 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 differentH, λc, and RMS values can be normalized similarly as described in Section 4.1. For each surface pair, between time duration 0 and 5τ, the mean aperture, h, and the contact ratio, c, are normalized as
(39)hnor=helashelas−refh,(40)cnor=celascelas−refc,where hnor is the normalized mean aperture, helas is the mean aperture right after the elastic deformation (before creep), helas−ref is the mean aperture right after the elastic deformation for surface pair 2 (H=0.8, RMS=50μm, and λc=1000μm), cnor is the normalized contact ratio, celas is the contact ratio right after the elastic deformation (before creep), and celas−ref is the contact ratio right after the elastic deformation for surface pair 2.Figure21 summarizes the normalized mean aperture and contact ratio for different H, λc, and RMS values. σ1 is fixed at 10 MPa, and the time duration is from 0 to 5τ.Figure 21
Normalized mean aperture and contact ratio of the seven surface pairs (shown in Table3). The time t 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.
(a)
Normalized mean aperture(b)
Normalized contact ratioThe 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

Figure22 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 5τ. Other input parameters are H=1.0, RMS=50μm, λc=1000μm, and σ=10MPa. 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 5τ. 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,H=1.0, RMS=50μm, λc=1000μm, σ=10MPa, and the time duration is 5τ.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.


## 4.1. Creep Deformation Normalization: The Maxwell Model

The creep curves (mean aperture and contact ratio versus normalized time) under differentH, λc, 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 2τ, the mean aperture, h, and the contact ratio, c, are normalized as
(37)hnor=helashelas−refh,(38)cnor=celascelas−refc,where hnor is the normalized mean aperture, helas is the mean aperture right after the elastic deformation (before creep), helas−ref is the mean aperture right after the elastic deformation for surface pair 2 (H=0.8, RMS=50μm, and λc=1000μm), cnor is the normalized contact ratio, celas is the contact ratio right after the elastic deformation (before creep), and celas−ref is the contact ratio right after the elastic deformation for surface pair 2.Figure20 summarizes the normalized mean aperture and contact ratio for different H, λc, and RMS values. σ1 is fixed at 10 MPa, and the time duration is from 0 to 2τ.Figure 20
Normalized mean aperture and contact ratio of the seven surface pairs (shown in Table3). The time t 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.
(a)
Normalized mean aperture(b)
Normalized contact ratioThe 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 differentH, λc, and RMS values can be normalized similarly as described in Section 4.1. For each surface pair, between time duration 0 and 5τ, the mean aperture, h, and the contact ratio, c, are normalized as
(39)hnor=helashelas−refh,(40)cnor=celascelas−refc,where hnor is the normalized mean aperture, helas is the mean aperture right after the elastic deformation (before creep), helas−ref is the mean aperture right after the elastic deformation for surface pair 2 (H=0.8, RMS=50μm, and λc=1000μm), cnor is the normalized contact ratio, celas is the contact ratio right after the elastic deformation (before creep), and celas−ref is the contact ratio right after the elastic deformation for surface pair 2.Figure21 summarizes the normalized mean aperture and contact ratio for different H, λc, and RMS values. σ1 is fixed at 10 MPa, and the time duration is from 0 to 5τ.Figure 21
Normalized mean aperture and contact ratio of the seven surface pairs (shown in Table3). The time t 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.
(a)
Normalized mean aperture(b)
Normalized contact ratioThe 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

Figure22 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 5τ. Other input parameters are H=1.0, RMS=50μm, λc=1000μm, and σ=10MPa. 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 5τ. 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,H=1.0, RMS=50μm, λc=1000μm, σ=10MPa, and the time duration is 5τ.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 H (0.6, 0.8, and 1.0), root mean square roughness RMS (50 μm, 100 μm, and 150 μm), and mismatch length λc (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)
WhenH increases, the average aperture decreases, and the contact ratio increases faster with time(2)
Whenλc 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 (H between 0.6 and 1.0, RMS between 50 and 150 μm, and λc between 1000 and 3000 μm), macroscopic stresses (between 5 and 20 MPa), and the time durations (2τ for the Maxwell model and 5τ 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 regionWhile 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.



---

*Source: 8879890-2020-12-10.xml*