Study on wave propagation across a single rough fracture by the modified thin-layer interface model

Journal of Applied Geophysics 110 (2014) 106–114

Contents lists available at ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

Study on wave propagation across a single rough fracture by the modified thin-layer interface model J.C. Li a,⁎, H.B. Li a, J. Zhao b a b

State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China Department of Civil Engineering, Monash University, Building 60, Clayton, VIC 3800, Australia

a r t i c l e

i n f o

Article history: Received 29 July 2014 Accepted 8 September 2014 Available online 18 September 2014 Keywords: Wave propagation Fracture Thin-layer interface model Fracture roughness Fracture weathering

a b s t r a c t Wave propagation across a fracture is dependent on the specific stiffness of the fracture, which is implicitly related to the fracture surface characteristics, such as roughness, matching and wall strength. In this study, a rough fracture with different void spaces and contact areas is modeled by two smooth surfaces separated by square column asperities with different heights, that is, the fracture is modeled as a thin-layer interface. The material properties of the asperities also vary to represent different weathering grades of fracture surfaces. The wave propagation equation is established by analyzing the interaction between a stress wave and the rough fracture. By back analysis, the relation between the normal stress and the closure of the fracture is derived, where the effects of the contact area, the fracture thickness, the frequency of the incidence and the weathering grade of the asperities on wave propagation are studied. The theoretical solutions of the thin-layer interface model (TLIM) are compared with those based on the zero-thickness interface model (ZTIM) which is commonly used currently. The study concludes that the TLIM is able to represent the behavior of fractures with different types of surface roughness and weathering. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Fractures or joints in a rock mass govern the physical and mechanical properties of the rock mass, including wave propagation in the rock mass. It is essential to determine the physical and mechanical behaviors of fractures, especially the seismic characteristics, which are associated with the problems in geology, geophysics, mining and earthquake engineering. The mechanical properties of fractures are dependent on the geometry of the void space, the contact area and the filling material. During various geological processes, such as weathering, shearing, loading and thermal cycles, fracture surfaces may be altered and the asperities on fracture surfaces may be crushed. The various geological processes influence the matching and roughness of fracture surfaces (Zhao, 1997). Fracture appearance changes as weathering proceeds (Ehlen, 2002). Barton and Choubey (1977) presented a joint roughness coefficient (JRC) to describe the roughness of joint surfaces. Bandis et al. (1983) carried out a laboratory investigation for natural unfilled joints to analyze the influence of weathering on the joint normal stiffness, and developed a hyperbolic relation between normal stress and closure of the joints. Brown and Scholz (1986) studied the effect of rough and mismatched surfaces on joint closure by comparing experimental and theoretical results. Zhao (1997) proposed a joint matching coefficient ⁎ Corresponding author. E-mail address: [email protected] (J.C. Li).

http://dx.doi.org/10.1016/j.jappgeo.2014.09.005 0926-9851/© 2014 Elsevier B.V. All rights reserved.

(JMC) based on the percentage of joint surfaces in contact. By coupling JMC and JRC, matched and mismatched joints can be represented. Based on static equilibrium and deformation compatibility, the Hertzian contact theory (Timoshenko and Goodier, 1951) for two curved smooth surfaces under a normal load was adopted by many authors to establish different joint contact theories. For example, Greenwood and Williamson (1966) derived an analytical half-joint model for the contact of a rough surface and a flat elastic surface. Greenwood and Tripp (1971) presented an analytical full-joint model for the contact of two rough surfaces. Brown and Scholz (1985) extended the model by Greenwood and Williamson (1966) and presented a model with general form for two random nominally flat elastic surfaces. From the Hertzian contact theory, Misra and Marangos (2011) derived the normal stiffness of asperities as a function of the normal deformation, which was then utilized to analyze the reflection and transmission of plane waves across fractures. In addition, Hudson et al. (1997) proposed an approximately circular contact model for faults with imperfect facial contact to analyze the mean transmission properties. Besides the spatial geometry of the contact area, the deformation of a joint is also dependent on the roughness of joint surfaces and the filling material. By using a similar numerical approach by Hopkins (2000), Pyrak-Nolte and Morris (2000) modeled a fracture with two half-spaces separated by cylindrical asperities to investigate the relation between the stiffness and fluid flow through the fracture. Based on the modified SHPB test equipments, studying for the dynamic normal properties of filled joints and longitudinal- (P-) wave propagation

J.C. Li et al. / Journal of Applied Geophysics 110 (2014) 106–114

across the filled fractures was conducted by Li and Ma (2009), Ma et al. (2011) and Wu et al. (2013). A fracture was usually modeled as an imperfectly bonded interface or a displacement discontinuous boundary to investigate wave propagation. For example, a quasi-static model was proposed by Baik and Thompson (1984) to analyze ultrasonic scattering. In their study, the imperfect interface was modeled as two half-spaces connected by distributed springs. Based on the displacement discontinuity method (DDM) (Miller, 1977; Schoenberg, 1980), many researchers analytically studied P-wave propagation across unfilled fractures (e.g. Fan et al., 2012; J. Zhao et al., 2006; Li et al., 2010, 2012; Pyrak-Nolte et al., 1990; Schoenberg, 1980; X.B. Zhao et al., 2006) and filled fractures (Fan and Wong, 2013; Wu et al., 2014). In the DDM, the stresses on the two surfaces of a fracture are assumed to be continuous while the displacements are discontinuous. To simplify the problem, the aperture of each fracture was considered to be zero in the foregoing analytical studies for wave propagation. The fracture in the existing analytical methods for wave propagation was modeled as a zero-thickness interface. In addition, some numerical methods (Jiao et al., 2005; Wu and Fan, 2014; Zhao et al., 2008) were adopted to simulate wave propagation across a jointed rock mass, where the joints were also modeled as zero-thickness and elastic interfaces. Different from the assumption of the zero-thickness interface model, Desai et al. (1984) presented a thin-layer interface concept for a joint or an interface between two solids. For the thin-layer interface concept, the joint is assumed as an equivalent solid or continuum medium with a finite and small thickness. Li et al. (2013, 2014) developed the concept to model rock joints as thin-layer interfaces. The thin-layer interface model for filled joints was approved to be affective to analyze wave propagation across a filled rock joint and to predict the mechanical property of the joint. This study addresses the effects of surface roughness on wave propagation across unfilled fractures and on the relationship between the closure and normal stress. The thin-layer interface model is improved to analyze the fracture treated as square column bars sandwiched between two plane surfaces. Except the longest columns, the cross sectional area of the columns (asperities) is assumed to follow the normal distribution. The materials of the columns are assigned to represent four weathering grades. The medium on both sides of the fracture is intact rocks with the same properties. Based on the equation of motion for one dimensional problem, the interaction between a P-wave and a welded interface of two contacting columns with different cross sections is analyzed. The wave propagation equation across the interface is established. The normal stress and closure for the fracture are calculated herein. Comparison is then carried out between the modified thin-layer interface model (TLIM) and the zero-thickness interface model (ZTIM). The causes of the discrepancy between two interface models and limitations of the present approach are to be discussed. The theoretical solution will be compared with experimental measurements in a forthcoming paper.

107

Left Right surface surface Column asperities

vI

vT

vR

H0 hj

Rock

Fracture

Rock

Fig. 1. Schematic illustration for a fracture with surface roughness.

not considered during wave propagation, the columns can be assumed to be uncorrelated. For simplification, the scattered waves on the two sides of the fracture are ignored and only 1D wave propagation problem is considered in the paper. When a P-wave impinges on the left side of the fracture, the wave propagates through the columns in contact with the other side. Wave reflection and transmission occur at the two contact surfaces, as shown in Fig. 1. For one unit cross sectional area of the rocks before and after the two surfaces, the initial contact area between the columns and planes is denoted by A0, which equals to the sum of the cross sectional area of the asperities with height H0. In this study, along the direction of wave propagation, the thickness of the fracture (i.e. the length of the longest columns) is divided into N sub-layers, and each sub-layer has the same thickness denoted as Δh. The height of the columns is denoted as Hj, which is equal to (N − j)Δh and j = [0, N]. H0 and HN are the heights of the longest and shortest columns, respectively. The columns and the adjacent rocks (blocks with plane surfaces) are assumed to deform elastically during wave propagation.

2.2. Basic equations for stresses and particle velocities In order to analyze the effect of fracture geometry on wave propagation, wave propagation across two connecting columns with thickness Δh is analyzed, as shown in Fig. 2. Based on the 1D wave propagation theory, there are two running waves propagating in two opposite directions in a continuous medium, which are termed as the right- and the left-running waves, and denoted

2. Theoretical formulations 2.1. Problem description A fracture with rough surfaces is schematically illustrated in Fig. 1. To simulate the fracture geometry in 3D problems, the rough fracture is represented by two plane surfaces separated by columns of the same cross section area (in 2D problems, they can be stripes) of different heights. The square columns denote the asperity of the fracture surfaces. The separation distance of the two plane surfaces is represented by H0, which also denotes the thickness of the fracture. The longest columns (asperities) initially in contact with both plane surfaces therefore have a height of H0. If the radial deformation of any column is

Fig. 2. Schematic illustration of wave propagation across welded elastic columns and blocks with different cross sections.

108

J.C. Li et al. / Journal of Applied Geophysics 110 (2014) 106–114

as RR and LR, respectively. When an RR or LR impinges on a discontinuous interface between two blocks, both reflection and transmission take − place (Bedford and Drumheller, 1994; Kolsky, 1953). v− r,j and vl,j in Fig. 2 + and v are the RR represent the RR and LR before the jth interface, v+ r,j l,j and LR after the jth interface, where j = [0, N]. The particle + velocities, v− j and vj , before and after the interface can be respectively expressed as −

−

−

þ

þ

þ

v j ðt Þ ¼ vr; j ðt Þ−vl; j ðt Þ and v j ðt Þ ¼ vr; j ðt Þ−vl; j ðt Þ:

ð1Þ

If the two columns besides the jth interface are welded, the force and the velocity at the discontinuous interface should be continuous, i.e. −

þ

−

þ

F j ðt Þ ¼ F j ðt Þ; v j ðt Þ ¼ v j ðt Þ:

ð2Þ

From Eq. (2), there is −

−

þ

þ

S j σ j ðt Þ ¼ S j σ j ðt Þ

ð3Þ

− − + + + − + where σ− j = σr,j + σl,j , σj = σr,j + σl,j , and the symbols Sj and Sj denote the cross sectional areas of the columns before and after the jth interface, respectively. According to the conservation of momentum on the wave front, the stress σ on the wave front is

σ ðt Þ ¼ ρcvðt Þ

−

þ

þ

vr; j ðt Þ−vl; j ðt Þ ¼ vr; j ðt Þ−vl; j ðt Þ:

i h i − − − − − − − þ þ þ þ þ þ þ ρ j c j vr; j ðt Þ þ ρ j c j vl; j ðt Þ S j ¼ ρ j c j vr; j ðt Þ þ ρ j c j vl; j ðt Þ S j

ð6Þ

− where ρ− j and cj are the density and wave velocity of the column + material before the jth interface in Fig. 2, ρ+ j and cj are the density and wave velocity of the material after the interface. At any time ti, Eqs. (5) and (6) can be expressed as the matrix



þ

vr; j ðt Þ − vl; j ðt Þ

(

 ¼A

−1

−

vr; j ðt Þ B þ vl; j ðt Þ

A¼

1

ð7Þ

þ þ þ 4 ρj cj Sj − − ρ− j cj Sj

1

3

2

1

5 4 −1 ; B ¼ 1

1

þ þ þ ρj cj Sj − − − − ρj cj Sj

3 5:

ð8Þ

Eq. (7) indicates wave propagation across an interface between two columns with different cross-sections. For wave propagation across the jth interface, the LR and RR before and after the interface can be obtained, respectively, according to the time-shifting function, that is −

þ

−

vR ðt Þ ¼ vl;0 ðt Þ:

ð10Þ

At the right surface of the fracture shown in Fig. 1, v+ l,N(t) is zero and the transmitted wave is þ

vT ðt Þ ¼ vr;N ðt Þ:

þ

−

vr; j ðt i Þ ¼ vr; j−1 ðt i −Δt Þ and vl; j ðt i Þ ¼ vl; jþ1 ðt i −Δt Þ

ð9Þ

where the subscripts “j − 1” and “j + 1” denote the (j − 1)th and (j + 1)th interfaces, respectively, and j = [1, N − 1]; Δt is the shifting time for the waves between two adjacent interfaces, i.e. Δt = Δh/c− j and Δt = Δh/c+ j for the RR and LR, respectively, in Eq. (9); Δh denotes the thickness of each column or sub-layer of the fracture shown in Fig. 1. Eq. (9) shows that the wave arriving at one interface can be obtained from the wave caused from the adjacent interfaces.

ð11Þ

At the left surface of the fracture, the strain of the rock εL(t) is εL ðt Þ ¼

vI ðt Þ vR ðt Þ þ cr cr

ð12Þ

and the particle velocity for the rock vL(t) is vL ðt Þ ¼ vI ðt Þ−vR ðt Þ:

ð13Þ

The strain and the particle velocity of the rock, εR(t) and vR(t), at the right surface of the fracture can respectively be given by

)

where A and B are the matrix parameters and 2

For wave propagation across the fracture shown in Fig. 1, if v− r,0(t) represents an incident wave vI(t), the reflected wave from the left surface of the fracture is

ð5Þ

Similarly, the continuity of the force on the interface is derived from Eqs. (3) and (4), h

2.3. Normal stress and closure of a fracture

ð4Þ

where ρ is the density of one medium, c is the P-wave velocity in the medium, and v(t) is the particle velocity. The velocity continuity at the interface can be derived from Eqs. (1) and (2), that is −

Eqs. (7) and (9) express a recursive computational method in time domain. The time interval dt = Δh/cr, and cr and ca are defined as the P-wave propagation velocities in rocks and the asperities, respectively. With the initial conditions vkr,j(t0) and vkl,j(t0) (k = − and +, and j = [0, N]) for all + interfaces and the boundary conditions v− r,0(ti) and vl,N(ti) at any time ti for the left and right surfaces of the fracture, Eqs. (7) and (9) with the matrix parameters A and B shown in Eq. (8) are applied to determine the particle velocities at the adjacent interfaces for times ti − Δt and ti + Δt. Through the recursive Eqs. (7) and (9), the transmitted wave v+ r,N(ti) after the right surface and the reflected wave v− l,0(ti) before the left surface can be obtained. The above analysis is available for incident wave propagating from the left surface. When the incident wave propagates from the right surface, the wave reflection at the free surface of the columns for non-zero aperture asperities should be considered with Eqs. (7) and (9) to analyze P-wave propagation across a fracture.

vT ðt Þ cr

ð14Þ

vR ðt Þ ¼ vT ðt Þ:

ð15Þ

ε R ðt Þ ¼

From the Hooke's law, the average normal stress on the fracture at time ti is σ ðt i Þ ¼

Er E ½ε þ εR  ¼ r ½vI ðt i Þ þ vR ðt i Þ þ vT ðt i Þ: 2 L 2cr

ð16Þ



When the strain rate of the fracture ε ðt Þ is obtained from 

ε ðt Þ ¼

1 1 ½v −vR  ¼ ½v ðt Þ−vR ðt Þ−vT ðt Þ; H0 L H0 I

ð17Þ

we can calculate the normal closure of the fracture from t0 to ti, Z ΔLðt i Þ ¼ εH0 ¼ H 0

ti 0

Z 

ε dt ¼

ti 0

½vI ðt Þ−vR ðt Þ−vT ðt Þdt

where H0 is the initial thickness of the fracture.

ð18Þ

J.C. Li et al. / Journal of Applied Geophysics 110 (2014) 106–114

0.04

processes, and the apertures are often several centimeters wide in highly weathered granitic rock. Based on a series of field investigation and laboratory tests, Zhao et al. (1994) studied different weathering patterns for granites. Some material parameters for the granites reported by Zhao et al. (1994) are adopted for the column asperities in Fig. 1, e.g., the Young's modulus Ea is set at 70, 35, 17.5 and 8.75 GPa for the weathering grades I, II, III and IV, and the corresponding density ρa is 2650, 2600, 2550 and 2450 kg/m3, respectively. The P-wave propagation velocities for the asperities of different weathering grades can be from the ffiffiffiffiffiffiffiffiffiffiffiffiffi pcalculated Young's modulus and the density, that is, ca ¼ Ea =ρa . The rocks on both sides of the fracture are assumed to be the granite with weathering grade I.

μ =0.1 μ =0.15 μ =0.2 μ =0.5

0.02

ω =50%

0.01

0.00 0

20

40

60

80

3.2. Geometry distribution

100

Asperity cross section area Aj (%) Fig. 3. Distribution of the cumulated cross sectional area Aj for non-zero apertures in one fracture.

3. Properties of fracture asperities in the analysis 3.1. Material property Fractures underwent various alteration processes, particularly weathering that reduces the fracture asperity strength, fracture surface roughness and matching (Zhao, 1997). By investigating the effects of weathering on fractures in granitic rocks, Ehlen (2002) observed that fracture surfaces in slightly weathered rock are stained, the apertures of major fractures change from narrow to wide during weathering

As illustrated in Fig. 1, some of the asperities (columns) between two plane surfaces are initially in contact with rocks (zero apertures) and some are not (non-zero apertures), whose heights are Hj (j = [0, N]). Here we only consider the rocks with one unit area. When j ranges from 1 to N, Hj also represents the distance from the jth interface between two adjacent sub-layers to the right surface. In Eq. (8), when j = 0, the cumulative cross sectional area S− 0 before the left surface shown in Fig. 1 is equal to one and the corresponding cumulative cross sectional area S+ 0 after the left surface is A0. When j = N, the cumulative cross sectional area before and after the right surface N−1

N

m¼1

m¼1

þ are S− N ¼ A0 þ ∑ Am and SN ¼ A0 þ ∑ Am ¼ 1, respectively. When

j = [1, N − 1], the corresponding cumulative cross sectional areas j−1

þ before and after the jth interface are S− j ¼ A0 þ ∑ Am and S j ¼ A0 þ m¼1

20

15

10

A0=40% A0=30%

5

A0=20% A0=10%

0 0.000

Normal stress (MPa)

Normal stress (MPa)

20

15

10

0.004

0.006

0.008

A0=30%

0.002

0.004

0.006

0.008

0.010

Closure (mm)

Closure (mm)

(a) Weathering grade I

(b) Weathering grade II 20

20

Normal stress (MPa)

A0=20% A0=10%

0.000

0.010

A0=40%

5

0

0.002

15

A0=40%

10

A0=30% A0=20%

5

A0=10% 0 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035

Closure (mm)

(c) Weathering grade III

Normal stress (MPa)

Probability density

0.03

109

15

10

A0=40% A0=30%

5

A0=20% A0=10%

0 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035

Closure (mm)

(d) Weathering grade IV

Fig. 4. The normal stress–closure curve for fractures with varying contact areas (H0 = 2 mm).

J.C. Li et al. / Journal of Applied Geophysics 110 (2014) 106–114

Weathering grade I Weathering grade II Weathering grade III Weathering grade IV

15000

1.0

Transmission coefficient

Normal stiffness (GPa/m)

110

10000

5000

0

10

20

30

40

0.9 Weathering grade I Weathering grade II Weathering grade III Weathering grade IV

0.8

TLIM ZTIM (Pyrak-Nolte et al., 1990)

0.7

0

1

2

3

4

5

6

7

8

Joint thickness (mm)

A0 (%)

Fig. 7. Relation between transmission coefficient and fracture thickness (A0 = 10%, f = 2000 Hz).

Fig. 5. Stiffness of fractures (H0 = 2 mm).

j

∑ Am , respectively. Since the size distribution for the non-zero

m¼1

aperture asperities can be normal or log-normal distribution (Pyrak-Nolte and Morris, 2000), Aj (j = [1, N]) for non-zero aperture

asperities is assumed to satisfy the normal distribution, as shown in Fig. 3, where ω and μ denote the location and the squared scale parameters of the normal distribution function, and ρA denotes the probability density. The effects of asperity weathering grade, initial contact area and fracture thickness on the wave propagation and the fracture mechanical properties will be analyzed in the following section.

1.00

Transmission coefficient

4. Effect of constant percentage of contact area

0.98

0.96

Weathering grade I Weathering grade II Weathering grade III Weathering grade IV

0.94

TLIM ZTIM (Pyrak-Nolte et al., 1990)

0.92

0.90

0

2000

4000

6000

8000

10000

Frequency (Hz)

(a) A0 =40%

Transmission coefficient

1.0 0.9

4.1. Normal stress and closure of fractures

0.8

Weathering grade I Weathering grade II Weathering grade III Weathering grade IV

0.7 0.6 0.5

At the initial condition, the asperities with height H0 are contacted with the rocks at the two surfaces of the fracture, as shown in Fig. 1. When an incident wave vI(t) impinges on the fracture, transmitted and reflected waves are created from the left surface and propagate in two opposite directions. The created transmitted wave propagates in the asperities with height H0 as a newly incident wave until it arrives at the right surface of the fracture, where newly transmitted and reflected waves are produced. During the wave propagation process, multiple transmitted and reflected waves are repeatedly given rise from the two fracture surfaces. The distance of the two fracture surfaces changes during wave propagation. The ratio of the cumulative contact area of the asperities to the adjacent rock with one unit area is defined to be the percentage of contact area. In this section, the cumulative contact area is assumed to be constant and equal to that at the initial condition, that is, the percentage of contact area is always A0 during wave propagation across the fracture. The incident wave is a harmonic wave with the amplitude I = 1 m/s, and the frequencies 100 and 1000 Hz, respectively.

TLIM ZTIM (Pyrak-Nolte et al., 1990)

0

2000

4000

6000

8000

10000

Frequency (Hz)

(b) A0 =10% Fig. 6. Relation between transmission coefficient and frequency of incident wave (H0 = 2 mm).

When the incident harmonic wave vI(t) impinges on the left surface of the fracture shown in Fig. 1, for a given weathering grade and a given percentage of contact area, the reflected and transmitted waves caused from the left and right surfaces of the fracture can be calculated from the wave propagation Eqs. (7) to (9). From Eqs. (16) and (18), the normal stress σ(t) and the corresponding closure ΔL(t) during the wave propagation process can be obtained. The stiffness kn is defined as the ratio of the normal stress to the closure, i.e. kn = σ(t)/ΔL(t). Fig. 4 shows the relation between the normal stress σ(t) and the closure ΔL(t) of fractures with thickness 2 mm, where four weathering grades of the asperities are considered and the frequencies of the incident harmonic wave are 100 and 1000 Hz, respectively. It can be observed from Fig. 4 that ΔL(t) linearly increases with the normal stress σ(t) on the fracture. From Fig. 4, kn is calculated and found to be constant when the weathering grade and the constant percentage of

J.C. Li et al. / Journal of Applied Geophysics 110 (2014) 106–114

Incident wave Transmitted wave Reflected wave

0.8 0.6 0.4 0.2

1.0

Particle velocity (m/s)

Particle velocity (m/s)

1.0

111

Incident wave Transmitted wave Reflected wave

0.8 0.6 0.4 0.2 0.0

0.0 0.00

0.05

0.10

0.15

-0.2 0.00

0.20

0.05

(a) Weathering grade I

0.4 0.2 0.0

1.0

Particle velocity (m/s)

Particle velocity (m/s)

0.6

0.20

(b) Weathering grade II

Incident wave Transmitted wave Reflected wave

0.8

0.15

Time (ms)

Time (ms)

1.0

0.10

Incident wave Transmitted wave Reflected wave

0.8 0.6 0.4 0.2 0.0 -0.2

-0.2 0.00

0.05

0.10

0.15

0.20

-0.4 0.00

0.05

(c) Weathering grade III

0.10

0.15

0.20

Time (ms)

Time (s)

(d) Weathering grade IV

Fig. 8. Incident, transmitted and reflected waves (f = 5000 Hz, H0 = 2 mm, A0 = 10%).

the contact area are given. In other words, kn is not related to the incident wave and the wave propagation process. Fig. 5 shows the relation between the stiffness kn and the constant percentage of the contact area for four weathering grades. It can be observed that kn for the fracture with asperities of weathering grade I and constant percentage of contact area A0 10% is the same with that of the fracture with weathering grade II and A0 20% or the fracture with weathering grade III and A0 40%, which is around 4500 GPa/m. 4.2. Transmission coefficient The transmission coefficient Tp − p is defined as the ratio of the magnitudes between the transmitted and incident waves. Pyrak-Nolte et al. (1990) derived the solution for the transmission coefficient based on the displacement discontinuity method for fractures with the zero-thickness interface model (ZTIM). Comparison is then carried out for Tp − p from the present approach and the existing method (Pyrak-Nolte et al., 1990), where the fracture stiffness kn shown in Fig. 5 is adopted for the existing method. Fig. 6 shows the variations of the transmission coefficient Tp − p with the frequency of the incident wave when the fracture thickness H0 is 2 mm. Fig. 6(a) and (b) plots the results for the constant percentages of contact area A0 = 40% and 10%, respectively. For A0 = 10% and f = 2000 Hz, the relation between the transmission coefficient and the fracture thickness is shown in Fig. 7. It can be seen from Figs. 6 and 7 that for the two interface models the tendencies of Tp − p with the variation of the fracture thickness or the frequency are very similar, that is, Tp − p decreases as the incident frequency or fracture thickness increases. Meanwhile, some discrepancy is observed in Tp − p between the two interface models, and the

discrepancy increases with larger thickness of the fracture or higher frequency of the incidence. For a thinner fracture or a lower frequency of the incidence, such as f less than 1000 Hz and H0 about 1 mm, Tp − p from the two methods are very close, while the discrepancy becomes much obvious for a wider fracture and a higher frequency. 5. Effect of increased percentage of contact area With increasing normal stress on the fracture, the aperture of the fracture will be closed, which results in the augment of the contact area between the two fracture surfaces. The incident wave in this section is assumed to be a half-cycle sinusoidal wave, i.e. vI(t) = I sin(2πft), where t ranges from 0 to 1/2f. 5.1. Calculation process Step 1: At the initial time t0, for a percentage of initial contact area A0, the reflected and transmitted waves before the left and after the right surfaces of the fracture can be calculated from Eqs. (7) to (9) when the incident wave is given. Step 2: The normal stress on the fracture and the closure of the two surfaces of the fracture can be obtained from Eqs. (16) and (18), respectively. Step 3: Compare the fracture closure ΔL with the aperture for the columns with height (N − j)Δh, where j = 1 for the first layer. If ΔL ≥ jΔh, the tops of the columns have been in contact with the left surface. Substitute the percentage of initial contact j

area by the new contact area A0 þ ∑ Am . m¼1

112

J.C. Li et al. / Journal of Applied Geophysics 110 (2014) 106–114

25 Constante contact area Increased contact area Curve fitting for

Normal stress (MPa)

Normal stress (MPa)

20

25

15

10

5

20

15

10

Constant contact area Increased contact area Curve fitting for

5

0 0.00

0.01

0.02

0 0.00

0.03

0.01

0.02

Closure (mm)

Closure (mm)

(a) H0 =2 mm

(a) H0=2 mm

0.03

20

Normal stress (MPa)

Normal stress (MPa)

15

10

Constant contact area Increased contact area Curve fitting for

5

0 0.00

0.02

0.04

0.06

0.08

0.10

0.12

0.14

15

10

0 0.00

Closure (mm)

Constant contact area Increased contact area Curve fitting for

5

0.02

0.04

0.06

0.08

Closure (mm)

(b) H0 =8 mm

(b) H0=8 mm

Fig. 9. The relation between the normal stress and closure of a fracture (weathering grade IV, A0 = 10%).

Step 4: Go back to Step 1 and continue the next time step to calculate the reflected and transmitted waves repeatedly, where the perj

centage of initial contact area A0 is substituted by A0 þ ∑ Am . m¼1

Step 5: Establish the relationship between the normal stress and the closure. 5.2. Normal stress and closure of fractures

Fig. 8 shows the reflected and transmitted waves for wave propagation across the asperities with four weathering grades, where the fracture thickness H0 is 2 mm, the initial percentage of contact area A0 is 10% and the wave frequency f is 5000 Hz. It can be observed from the four figures that the incident wave is almost transmitted across the fracture for the weathering grade I. When the weathering grade changes from I to IV, the incident wave attenuates more and the transmitted wave form gradually becomes wider. For the fracture with weathering grade IV, the relation between the normal stress and the closure of the fracture is shown in Figs. 9 and 10 for the percentage of initial contact area A0 10% and 20%, respectively, where the two random parameters μ and ω are 0.5 and 50%, respectively.

Fig. 10. The relation between pressure and closure of a fracture (weathering grade IV, A0 = 20%).

The relation between the normal stress and closure of a fracture with constant and increased contact areas can be found in Figs. 9 and 10, where the results for a constant contact area of a fracture have been studied in Section 4.1. If the contact area varies with the change of the normal load, the closure increases nonlinearly with increasing normal stress, while for a constant contact area the relation between the normal stress and the closure is linear. In another word, the increased contact area causes nonlinear normal stress–closure curves. If the effect of increased contact area is considered, the initial slope of the normal stress–closure curve is almost the same with the stiffness of the fracture with constant contact area. Meanwhile, for a given normal stress on a fracture, the closure of the fracture under the effect of increased contact area is smaller than that of the fracture with constant contact area. In another word, the two surfaces of the fracture are closer if the contact area is constant than that if the contact area increases with increasing normal stress. The nonlinear stress–closure relation of the fracture shown in Figs. 9 and 10 can be fitted with hyperbolic or logarithmic curves (Bandis et al., 1983; Brown and Scholz, 1986). In the present study, the relation between the normal stress and the closure under the effect of increased contact area is fitted with hyperbolic curves, which were used to describe the mechanical properties of nonlinear fractures by X.B. Zhao et al. (2006). The stiffness of the fracture calculated in Fig. 5 is adopted as the initial stiffness kni in the hyperbolic curve, and the maximum

J.C. Li et al. / Journal of Applied Geophysics 110 (2014) 106–114

kni (GPa/m) Lmax/H0

H0 = 8 mm A0 = 20%,

A0 = 10%,

A0 = 20%,

556 0.072

1113 0.07

140 0.08

280 0.12

closure L max is calculated and shown in Table 1. The corresponding fitting curves are also illustrated in Figs. 9 and 10 for fractures with two initial contact areas, respectively. Figs. 9 and 10 show that the increase of the normal stress on the fracture results in gradual closure of the fracture. When an incident wave impinges on the fracture, the tops of the asperities with height H0 contact first with the rocks. Then the aperture between the asperities with height hj and the left surface of the fracture is reduced with increasing normal load. If the closure exceeds some of the aperture, more asperities will be in contact with the left surface of the fracture, resulting in more contact area. 5.3. Transmission coefficient The variation of transmission coefficient with the frequency of the incidence is shown in Fig. 11, where two random parameters, μ 0.5 and ω 50%, are adopted for the distribution of the asperities with weathering grade IV. Similar tendency can be observed in Figs. 6 and 11 that the transmission coefficient calculated from the two interface models for fractures drops with increasing frequency of the incidence. For a given frequency, the transmission coefficient for wave propagation across the fracture with H0 = 2 mm is greater than that for the fracture with H0 = 8 mm. This is caused by the mechanical property of the fractures, that is, the slope of the normal stress versus closure curve for the fracture with H0 = 8 mm is smaller than that for the fracture with H0 = 2 mm, which has been shown in Figs. 9 and 10. Comparison also shows that the transmission coefficient calculated from the present approach is smaller than that from the DDM coupled with the method of characteristic (X.B. Zhao et al., 2006), in which the fracture is modeled as a zero-thickness interface. The transmitted wave is the superposition of two parts: one is the result from wave direct propagation across the fracture, and the second part is the result from the wave reflections between the two surfaces of the columns. These two parts arrive at different times, which may cause the difference between the transmission coefficients based on the two interface models. Moreover, this difference becomes more obvious when the fracture thickness or the frequency of the incident wave is bigger. When the initial contact area is 10% and the fracture thickness is 2 mm, for a given frequency of the incident wave, the transmission coefficients Tp − p from the two interface models shown in Fig. 11(a) are greater than those in Fig. 6. For example, in Fig. 11 Tp − p is 0.861 and 0.804 from the TLIM and ZTIM, respectively, while Tp − p in Fig. 6 is 0.793 and 0.6 from the two interface models, respectively. It can be concluded that the theoretical results of Tp − p considering the effect of increased contact area of the fracture are greater than those only including the effect of initial contact area of the fracture.

H0= 2 mm

1.0

Transmission coefficient

H0 = 2 mm A0 = 10%,

Fig. 11 shows that the theoretical results based on TLIM and ZTIM are almost the same when the fracture thickness is 2 mm and the incident frequency is less than 3000 Hz, or when the fracture thickness is 8 mm and the incident frequency is less than 1000 Hz. Otherwise, the discrepancy of the transmission coefficient Tp − p between two methods appears and increases with increasing fracture thickness. Hence, if the effect of a fracture thickness on wave propagation can't be omitted, the TLIM should be adopted. Fig. 5 shows that the fractures with percentages of contact area 10%, 20% and 40% for weathering grades I, II and III, respectively, have the same stiffness if the contact area remains constant and is equal to the initial value. If the force on the two surfaces of the fracture is balanced, the normal stress σR on the fracture surface is equal to the normal stress σJ on the asperities multiplied by the initial percentage of contact area, that is, σR = σJ ⋅ A0. As the closure of the fracture equals to the relative normal displacement of the asperities with height H0 , i.e. ΔL = σ J H 0/E J , there is ΔL = σR H 0 / (EJ A0). If the contact area is assumed to be constant during loading and equal to the initial value, the stiffness of the fracture kn can be predicted as k n = σ r /ΔL = E J A0/H 0, indicating that k n is directly proportional to EJ A0. Therefore, for a given thickness of a fracture, the stiffness for different contact areas and weathering grades are identical if the values of EJ A0 are the same, which can be observed in Fig. 5.

0.8

H0= 8 mm

0.6

0.4

TLIM ZTIM (Zhao XB et al., 2006) TLIM ZTIM (Zhao XB et al., 2006)

0.2

0.0

0

4000

6000

8000

10000

(a) A0 =10% 1.0

H0= 2 mm

0.8

H0= 8 mm 0.6

0.4

TLIM ZTIM (Zhao XB et al., 2006) TLIM ZTIM (Zhao XB et al., 2006)

0.2

6. Discussions 0.0

During wave propagation across the asperities of a fracture, multiple reflections between two surfaces of the fracture occur. For a given incidence, a wider layer gives rise to more time delay for wave propagation. Hence, for the fracture with H0 = 8 mm, the discrepancy for Tp − p between two interface models is bigger than those for the fracture with H0 = 1 and 2 mm shown in Figs. 7 and 11, respectively.

2000

Frequency (Hz)

Transmission coefficient

Table 1 Parameters of the hyperbolic fitting curves for the normal stress–closure relation (for fractures with weathering grade IV).

113

0

2000

4000

6000

8000

10000

Frequency (Hz)

(b) A0 =20% Fig. 11. Transmission coefficients based on TLIM and ZTIM (weathering grade IV).

114

J.C. Li et al. / Journal of Applied Geophysics 110 (2014) 106–114

7. Conclusions The present paper intends to find a new method to study wave propagation across rough fractures if the mechanical behavior of the fracture is unavailable in advance. In this paper, a rough fracture is modeled by two smooth surfaces separated by square column asperities with different heights, that is, the fracture is modeled as a thin-layer interface. Based on the thin-layer interface model (TLIM), the effect of the fracture thickness is demonstrated by the wave propagation equation Eq. (9). Compared to the zero-thickness interface model (ZTIM), the TLIM can be applied to analyze P-wave propagation across a joint with rough surface without knowing its dynamic property in advance. The material properties of the asperities are assigned to represent four weathering grades. The present study shows that the transmission coefficient is related to both the fracture thickness and the frequency of incidence, and decreases with increasing values of these two parameters. The transmission coefficients are also related to the material properties of the asperities. For a given fracture thickness, the transmission coefficient is the highest for wave propagation across fresh fractures, and gradually decreases when the weathering grade ranges from II to IV. By back analysis, the relation between the normal stress and the closure of the fracture can be estimated from the present approach. The theoretical results show that the relation between the normal stress and the closure is linear if the contact area remains constant, while the closure of the fracture nonlinearly increases with increasing normal stress on the fracture if the contact area varies during wave propagation. Nature fractures present different surface profiles on each side and different degrees of matching. When the joint roughness coefficient (JRC) and joint matching coefficient (JMC) are considered to describe the geometry and contact of the fractures, the approach presented in the paper can then be applied in reality. Further exploration is underway to extend the current approach for more complicated geometry of fractures and effect of incidence angle, wave conversions, effective shear stiffness and correlation of the asperities on wave propagation. Acknowledgments The authors acknowledge the support from the Chinese National Science Research Fund (41272348, 51025935) and the Major State Basic Research Project of China (2010CB732001). References Baik, J.M., Thompson, R.B., 1984. Ultrasonic scattering from imperfect interfaces: a quasi-static model. J. Nondestruct. Eval. 4 (3–4), 177–196. Bandis, S.C., Lumsden, A.C., Barton, N.R., 1983. Fundamental of rock joint deformation. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 20, 249–268. Barton, N.R., Choubey, V., 1977. The shear strength of rock joints in theory and practice. Rock Mech. 10, 1–54. Bedford, A., Drumheller, D.S., 1994. Introduction to Elastic Wave Propagation. John Wiley, New York. Brown, S.R., Scholz, C.H., 1985. Closure of random elastic surfaces in contact. J. Geophys. Res. 90 (B7), 5531–5545.

Brown, S., Scholz, C., 1986. Closure of rock joints. J. Geophys. Res. 91 (B5), 4939–4948. Desai, C.S., Zaman, M.M., Lightner, J.G., Siriwardance, H.J., 1984. Thin-layer element for interfaces and joints. Int. J. Numer. Anal. Methods Geomech. 8 (1), 19–43. Ehlen, J., 2002. Some effects of weathering on joints in granitic rocks. Catena 49, 91–109. Fan, L.F., Wong, L.N.Y., 2013. Stress wave transmission across a filled joint with different loading/unloading constitutive relationships. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 60, 227–234. Fan, L.F., Ma, G.W., Wong, L.N.Y., 2012. Effective viscoelastic behaviour of rock mass with double-scale discontinuities. Geophys. J. Int. 191 (1), 147–154. Greenwood, J., Tripp, J.H., 1971. The contact of two nominally flat rough surfaces. Proc. Inst. Mech. Eng. 85, 625–633. Greenwood, J.A., Williamson, J.B.P., 1966. Contact of nominally flat surfaces. Proc. R. Soc. Lond. A29, 300–319. Hopkins, D.L., 2000. The implication of joint deformation in analyzing the properties and behavior of fractured rock masses, underground excavations, and faults. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 37, 175–202. Hudson, J.A., Liu, E., Crampin, S., 1997. The mean transmission properties of a fault with imperfect facial contact. Geophys. J. Int. 129 (3), 720–726. Jiao, Y.Y., Fan, S.C., Zhao, J., 2005. Numerical investigation of joint effect on shock wave propagation in jointed rock masses. J. Test. Eval. 33 (3), 197–203. Kolsky, H., 1953. Stress Waves in Solids. Clarendon Press, Oxford. Li, J.C., Ma, G.W., 2009. Experimental study of stress wave propagation across a filled rock joint. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 46 (3), 471–478. Li, J.C., Ma, G.W., Zhao, J., 2010. An equivalent viscoelastic model for rock mass with parallel joints. J. Geophys. Res. 115 (B03305). Li, J.C., Li, H.B., Ma, G.W., Zhao, J., 2012. A time-domain recursive method to analyze transient wave propagation across rock joints. Geophys. J. Int. 188 (2), 631–644. Li, J.C., Wu, W., Li, H.B., Zhu, J.B., Zhao, J., 2013. A thin-layer interface model for wave propagation through filled rock joints. J. Appl. Geophys. 91, 31–38. Li, J.C., Li, H.B., Jiao, Y.Y., Liu, Y.Q., Xia, X., Yu, C., 2014. Analysis for oblique wave propagation across filled joints based on thin-layer interface model. J. Appl. Geophys. 102, 39–46. Ma, G.W., Li, J.C., Zhao, J., 2011. Three-phase medium model for filled rock joint and interaction with stress waves. Int. J. Numer. Anal. Methods Geomech. 35, 97–110. Miller, R.K., 1977. An approximate method of analysis of the transmission of elastic waves through a frictional boundary. J. Appl. Math. 44 (4), 652–656. Misra, A., Marangos, O., 2011. Rock-joint micromechanics: relationship of roughness to closure and wave propagation. Int. J. Geomech. 11, 431–439 (Special Issue: Material and Computer Modeling). Pyrak-Nolte, L.J., Morris, J.P., 2000. Single fractures under normal stress: the relation between fracture specific stiffness and fluid flow. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 20, 245–262. Pyrak-Nolte, L.J., Myer, L.R., Cook, N.G.W., 1990. Anisotropy in seismic velocities and amplitudes from multiple parallel fractures. J. Geophys. Res. 95 (B7), 11345–11358. Schoenberg, M., 1980. Elastic wave behavior across linear slip interfaces. J. Acoust. Soc. Am. 68 (5), 1516–1521. Timoshenko, S., Goodier, J.N., 1951. Theory of Elasticity. McGraw-Hill, New York. Wu, Z.J., Fan, L.F., 2014. The numerical manifold method for elastic wave propagation in rock with time-dependent absorbing boundary conditions. Eng. Anal. Bound. Elem. 46, 41–50. Wu, W., Li, J.C., Zhao, J., 2013. Seismic response of adjacent filled parallel rock fractures with dissimilar properties. J. Appl. Geophys. 96, 33–37. Wu, W., Li, J.C., Zhao, J., 2014. Role of filling materials in a P-wave interaction with a rock fracture. Eng. Geol. 172, 77–84. Zhao, J., 1997. Joint surface matching and shear strength part A: joint matching coefficient. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 34, 173–178. Zhao, J., Broms, B.B., Zhou, Y., Choa, V., 1994. A study of the weathering of the Bukit Timah granite, part B: field and laboratory investigation. Bull. Int. Assoc. Eng. Geol. 50 (1), 105–111. Zhao, J., Zhao, X.B., Cai, J.G., 2006a. A further study of P-wave attenuation across parallel fractures with linear deformational behavior. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 43, 776–788. Zhao, X.B., Zhao, J., Cai, J.G., 2006b. P-wave transmission across fractures with nonlinear deformational behavior. Int. J. Numer. Anal. Methods Geomech. 30 (11), 1097–1112. Zhao, X.B., Zhao, J., Cai, J.G., Hefny, A.M., 2008. UDEC modeling on wave propagation across fractured rock masses. Comput. Geotech. 35 (1), 97–104.