Numerical Simulation for Rock Fracture Viscoelastic Creep under Dry Conditions

In many rock engineering projects such as hydrocarbon extraction and geothermal energy utilization, the hydraulic and mechanical behavior of rock fractures signi ﬁ cantly a ﬀ ects the safety and pro ﬁ tability of the project. In ﬁ eld 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 a ﬀ ected by the rock fracture surface geometry. However, there is not much literature systematically studying the e ﬀ ect of surface geometry on rock fracture creep under dry conditions. This paper presents the results of a numerical study considering the e ﬀ ect 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.


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][2][3][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][6][7][8][9][10][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][13][14][15][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][18][19][20][21][22][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][25][26][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. 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][29][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].

Methodology
The surface roughness of rock fractures typically follows a self-affine fractal distribution ( [29,31]). A self-affine surface can be described as 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 (A 2 ) 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 CðqÞ can be related to the wavenumber (q): Decompose via Fourier Transform

Random Function
Fourier Components 2A 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).
q max q min Figure 2: Schematic of the power spectrum. The horizontal axis is the frequency q (1/λ), and the vertical axis is the power C (A 2 ). The unit for q is (1/length).
The upper frequency end, q max , is defined as q max = 2π/λ 1 , where λ 1 is the smallest representable or measurable size. The lower frequency end, q min , is defined as q min = 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][29][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.  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.  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, σ ini , is calculated. Then, the surface height is scaled linearly as where z ini is the initial surface height, z final 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). 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.

Methodology for Calculating Fracture Elastic
Deformation. The method for calculating fracture elastic (1,1) x y-direction x-direction y Figure 6: Illustration of the grid discretization and discretized Boussinesq equation. Indices i 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 p k,l is applied on grid ðk, lÞ, and ðu e Þ i,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.
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 h ðx, yÞ is defined as where δ is the rigid-body movement between two surfaces under applied stress in the direction normal to the fracture  Figure 7: (a) Schematic for creep compliance under a constant stress σ 0 . Jð0 + Þ represents the creep compliance at the moment when σ 0 is applied, and Jð0 + Þ is equivalent to the elastic modulus. (b) Schematic for relaxation modulus under a constant strain ε 0 . E r ð0 + Þ represents the relaxation modulus at the moment when ε 0 is applied, and E r ð0 + Þ is equivalent to the elastic modulus.
Start. Obtain load history F(t), initial aperture h i, j , and creep compliance function J(t) Calculate J(( -') t) for each time step and the stiffness matrix K i, k, j, l .
Calculate initial contact state (elastic deformation); obtain pressure and aperture field.
Increase and the displacement field.
. Calculate the pressure p ( ) Stop. Export the results.

No
Is smaller than N t ? surface (compression is defined as positive), h 0 ðx, yÞ is the initial aperture, and u e ðx, 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: where pðx, 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 displacement u e ðx, yÞ can be calculated from the Boussinesq and Cerrutti (1885) equation: where Kðx, y, x′, y′Þ is the influence matrix for normal displacement at surface location Bðx, yÞ due to a normal unit stress applied at location Aðx′, y′Þ and pðx′, y′Þ is the normal stress applied at location Aðx ′ , y ′ Þ. Figure 5 shows the schematic of the Boussinesq and Cerrutti equation. The influence matrix Kðx, y, x ′ , y ′ Þ is defined as 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 Figure 6, the surface is discretized into a regular grid: where x i and y j are x and y coordinates, respectively; Δx and Δy are the grid spacings in the x and y directions, respec-tively; and N and M are the total number of grids in the x and y directions, respectively. The aperture field can be calculated as The boundary conditions become p i,j > 0 and h i,j = 0, ð11Þ Love [33] discretized the Boussinesq and Cerrutti equation as  (30) and (31).
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 calcu-lates 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, ðu e Þ i,j can be expressed as where FFT -1 means inverse fast Fourier transform (IFFT). In mathematics, FFT changes the convolution (see Equation (13)) to a simple multiplication (see Equation (16)). If the FFT method is not used, calculating Equation (10) requires N 2 × M 2 operations, while by using the FFT method, calculating Equation (13) requires N × M × ln ðN × 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: where F is the applied external force on the two fracture surfaces and p k,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.  [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 E r ðtÞ and JðtÞ are the relaxation modulus function and creep compliance function, respectively. E r ðtÞ (unit: pressure) describes the time-dependent stress response to a step-change in strain, and JðtÞ (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 vari-able). Figure 7 explains the creep compliance and relaxation 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 1/2G (in Equation (8)) is replaced by the creep compliance Jðt − τÞ.
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 ( Ð t 0 ð∂pðx ′ , 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 N t 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 ∂pðx ′ , y ′ , τÞdτ/∂τ in Equation (20) can be expressed as pðx ′ , y ′ , τ + dτÞ − p ðx ′ , y ′ , τÞ. Based on the sufficiently small time step Table 1: Natural rock joint profiles analysis results (obtained from [28]).   (20) can be rewritten as a summation:

Rock joint specimen
where α = 1, 2, ⋯, N t . Within each time duration, in a given surface, the pressure field does not change with time. Therefore, ½pðx ′ , y ′ , α ′ Þ − pðx ′ , y ′ , α ′ − 1Þ 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: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 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 where G 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 ðu e Þ α′ , 2GJððα − α ′ ÞΔtÞK i,j,k,l , and ðp k,l,α′ − p k,l,α′−1 Þ in Equation (28) are equivalent to ðu e Þ i,j , K i,k,j,l , and p k,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 history FðtÞ, initial discretized aperture field h i,j , and the creep compliance function JðtÞ for the viscoelastic material

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. The Maxwell model consists of a spring and a dashpot. The spring represents elasticity, with a shear modulus of G; 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 t is the time. The relaxation time for the Maxwell model, T, is defined as The SLS model consists of two springs and one dashpot. The two springs represent elasticity, with an elastic modulus of G 1 and G 2 , 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 The creep compliance function can be expressed as  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). where t is the time. The relaxation time for the SLS model, T, 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 p is the pressure, r is the distance from the center of the contacting region, t is the total time duration, aðtÞ 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. Lee and Radok [39] also suggested that if the flat surface satisfies the SLS model, the pressure distribution field can be expressed as where p is the pressure, r is the distance from the center of the contacting region, t is the total time duration, aðtÞ is the radius of the contacting region, R is the spherical diameter, G 1 and G 2 are the shear moduli, η is the viscosity, and P is the total load. The physical parameters are also illustrated in Figure 10.  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). 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.
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.

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 (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 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).    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 (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 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. Figures 13-15 show the effect of different H, λ c , and RMS values on the aperture field, respectively. The following conclusions can be drawn: (1) From Figure 13, as H decreases, the mean and standard deviation of the aperture increase (2) From Figure 14, as λ c increases, the mean and standard deviation of the aperture increase   Table 3). 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.
Note: "↓" means decrease and "↑"means increase. RMS value (this can be explained from Equation (4) in Section 2.1) 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.

Creep Simulation Results for the Maxwell Model.
The creep deformations of the seven surface pairs are calculated using the Maxwell model. Table 4 summarizes the input parameters. The parameters are obtained from the testing results of the Vaca Muerta Shale published by Mighani et al. [41].
The following parameters are also defined: (1) Macroscopic stress σ 1 = F/A, where F is the total force applied on the fracture surface and A is the area of the fracture surface (10 mm * 10 mm)  The following conclusions can be drawn: (1) Based on the curves of surface pairs 1, 2, and 3, when H 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% Table 5 summarizes the above conclusions. 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 H = 1:0, RMS = 50 μm, λ c = 1000 μm, and σ = 10 MPa, and the creep time duration = 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).

Creep Simulation Results for the Standard Linear Solid
(SLS) Model. The creep deformations of the seven surface pairs are also calculated by the SLS model. Table 6 summarizes the input parameters. The parameters are obtained from the testing results of the Vaca Muerta Shale published by Mighani et al. [41]. Figure 18 summarizes the effect of different 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.
The following conclusions can be drawn: (1) Based on the curves of surface pairs 1, 2, and 3, when H 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 rate Table 7 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. 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 H = 1:0, RMS = 50 μm, λ c = 1000 μm, and σ = 10 MPa, and the creep time duration = 5τ. The left figure corresponds to where h nor is the normalized mean aperture, h elas is the mean aperture right after the elastic deformation (before creep), h elas−ref is the mean aperture right after the elastic   Table 3). 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.
Note: "↓" means decrease and "↑" means increase.    Figure 20 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τ.
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.

Creep Deformation Normalization:
The SLS Model. The creep curves (mean aperture and contact ratio versus normalized time) under different H, λ 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 where h nor is the normalized mean aperture, h elas is the mean aperture right after the elastic deformation (before creep), h elas−ref is the mean aperture right after the elastic deformation for surface pair 2 (H = 0:8, RMS = 50 μm, and λ c = 1000 μm), c nor is the normalized contact ratio, c elas is the contact ratio right after the elastic deformation (before creep), and c elas−ref 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 H, λ c , and RMS values. σ 1 is fixed at 10 MPa, and the time duration is from 0 to 5τ.   Table 3). 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.  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.

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 5τ. Other input parameters are H = 1:0, RMS = 50 μm, λ c = 1000 μm, and σ = 10 MPa. 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).
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.   Table 3). 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. 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.

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