# Effect of Surface Morphology on the Dissipation During Shear and Slip Along a Rock-Rock Interface that Contains a Visco-elastic Core

Luiza Angheluta, Thibault Candela, Joachim Mathiesen, Francois Renard

► To cite this version:

Luiza Angheluta, Thibault Candela, Joachim Mathiesen, Francois Renard. Effect of Surface Morphology on the Dissipation During Shear and Slip Along a Rock-Rock Interface that Contains a Visco-elastic Core. Pure and Applied Geophysics, 2011, 168 (12), pp.2335-2344. 10.1007/s00024-011-0272-8. insu-00679237

HAL Id: insu-00679237

https://insu.hal.science/insu-00679237

Submitted on 21 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effect of Surface Morphology on the Dissipation During Shear and Slip Along a Rock–Rock Interface that Contains a Visco-elastic Core

LUIZA ANGHELUTA,<fnref n="1" /> THIBAULT CANDELA,<fnref n="2" /> JOACHIM MATHIESEN,<fnref n="1" />,<fnref n="3" /> and FRANÇOIS RENARD<fnref n="1" />,<fnref n="2" />

*Abstract*—High resolution topography measurements of the Vuache–Sillingy fault (Alps, France) reveal a characteristic roughness of the fault zone. We investigate the effect of roughness on the rheology of a planar shear configuration by using a model system consisting of a visco-elastic layer embedded into a rigid solid. The model is discussed in the context of several geological cases: a damage fault zone, a fault smeared with a clay layer, and a shear zone with strain weakening. Using both analytical approaches and finite element simulations, we calculate to linear order the relation between wall roughness and the viscous dissipation in the fault zone as well as the average shear rate.

## 1. Introduction

In the Earth’s crust, many systems involve deformation along pre-existing interfaces. The fact that these interfaces, loaded by normal or shear stress or a combination of both, locally alter the stress transmission may lead to several important geological phenomena that occur on all scales. At the grain scale, stress is concentrated along grain boundaries where localized deformation may occur either by dissolution-precipitation processes or by surface diffusion (ANGHELUTA *et al.*, 2008; RUTTER, 1976). At the outcrop scale, localized dissolution along existing planes leads to the formation of pressure solution seams or stylolites (ARTHAUD and MATTAUER, 1969; RAILSBACK and ANDREWS, 1995). At the lithospheric scale, shear displacement is also widely observed,

along faults in the upper-crust, or in shear zones at greater depth – for a review, see for example (CHESTER *et al.*, 2004). In a recent work, we have shown that such mechanical interfaces can become morphologically unstable and roughen with time (ANGHELUTA *et al.*, 2010). All these observations involve deformation past a mechanical discontinuity.

In the present study, we characterize the interaction between the morphology of the interface and the corresponding shear or sliding resistance in a set-up where a visco-elastic layer is embedded between two rigid bodies. For this, we use analytical approaches and finite element simulations in two dimensions to estimate how surface roughness affects the effective shear flow properties. We consider here two elastic solids, separated by a fluid layer of finite thickness *H*. The contact layer is approximated by a Maxwell visco-elastic rheology representing a core of a fault zone where dissipative processes occur (CHERY *et al.*, 2004). The model could be applied to various geological cases: faults that contain a damage zone, faults smeared with a clay layer, shear zones with strain weakening, or stylolites with a clay interface. Note also that a similar approach has been applied in studies of basal flow of glaciers along rough surfaces (KAMB, 1970) and the formation of residual stresses due to slip on wavy faults (CHESTER and CHESTER, 2000; SAUCIER *et al.*, 1992). Here we consider a full visco-elastic rheology and extend on previous analyses by performing a numerical modelling in geometries sampled in the field.

The paper is divided into three sections. In Sect. 2, we present geological observations of rough fault surfaces and how high resolution topography measurements can help characterizing the self-affine property of the slip surface. In Sect. 3, we consider

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹ Physics of Geological Processes, University of Oslo, P.O. Box 1048, Blindern 0316 Oslo, Norway. E-mail: mathies@nbi.dk</footnote>

<footnote n="2">² University Joseph Fourier, GrenobleI/CNRS/OSUG/ISTerre, Grenoble, France.</footnote>

<footnote n="3">³ Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, 2100 Copenhagen, Denmark.</footnote>
relations between this roughness and the rheology of the fault by introducing a simple model of a fault gouge. The effect of small perturbations of the fault surface on the fault dissipation is studied within the linear analysis. The analysis can be extended into the finite amplitude regime using numerical simulations. We compute numerically the stress distribution near a rough surface and the average flow rate in the Newtonian limit and compare it with the linear theoretical analysis. Finally, Sect. 4 offers concluding remarks.

## 2. Rough Geological Interfaces

### 2.1. Observations of Rough Faults

In the upper crust, many processes happen at interfaces which generally have mechanical properties different than the bulk. Often these interfaces are not flat and have developed complex rough morphologies. Active fault surfaces are known to show corrugation at all scales (POWER et al., 1987), see Fig. 1. Their complex geometry results from the interplay of abrasion processes, crack multi-branching and damage during rupture propagation (POWER et al., 1988; SAGY et al., 2007), and healing processes during the interseismic periods (RENARD et al., 2000). With the recent development of accurate Light Detection and Ranging (LiDAR) laser devices, it is possible to measure accurately the topography of such surfaces at all scales (Fig. 2). These measurements show a non-flat topography from the micrometer scale to the scale of several tens of meters (CANDELA et al., 2009). Moreover, detailed surface rupture mapping on the kilometer scale also reveals non planar geometry (KLINGER., 2010). In several cases, the geometry of slip surfaces has been observed to scale with different power-law exponents in the direction of slip and perpendicular to it (POWER et al., 1987; RENARD et al., 2006).

In such faults, a weaker layer or damage zone [from meters to tens of meters wide (SHIPTON et al., 2006)] is sandwiched between less damaged rock bodies (CAINE et al., 1996; CHESTER and LOGAN, 1986). This zone may creep slowly, either permanently (i.e. the creeping segment of the San Andreas Fault, north of Parkfield, CA, USA) or for a short period after a major earthquake, where afterslip is often measured [for a

review on afterslip processes see, for example, (PRITCHARD and SIMONS, 2006)]. These observations indicate that the motion along the fault is able to overcome the morphological roughness asperities, usually without any emission of seismic waves. Therefore, some time-dependent dissipative deformation mechanisms are at work during aseismic slip.

At the outcrop scale, it is also common to observe, in faulted sedimentary basins, that clay layers smear a fault interface (EGHOLM and CLAUSEN, 2008) and lubricate it (see Fig. 1d). In this case, the clay layer, with a visco-plastic behavior, is often strongly dragged, sheared and internally deformed during slip, as shown experimentally by CUISIAT and SKURVEIT (2010). In general, the clay layer must deform along the interface such to overcome the possible roughness of fault interfaces.

### 2.2. Roughness of the Vuache–Sillingy Strike-Slip Fault

The roughness of a slip surface of the Vuache–Sillingy has been measured for spatial wavelengths covering more than 6 decades. This fault is considered as a model system and the results could be extrapolated to other faults which show similar scaling properties (POWER et al., 1987; SAGY et al., 2007). This strike-slip fault, with a small normal component, is located near Annecy in the French Alps (RENARD et al., 2006) and exposes well-protected slip surfaces in carbonate rocks (Fig. 2a). The fault surface has been measured at different scales using three high resolution devices. At the outcrop scale, the morphology of the slip surface was measured using a LiDAR device [see (RENARD et al., 2006) for more details]. In the present study, we complement these outcrop data with laboratory scale measurements. We have measured the topography of several hand samples using a laser profilometer, with a height resolution down to 1 µm and spatial increments of 30 µm; and a white light interferometer, with a height resolution down to 1 nanometer and spatial increments of 0.5 µm. The results of each topography measurement is a Digital Elevation Model (DEM) of the slip surface (Fig. 2b).

The slip surface shows corrugations at all scales and a slight anisotropy is observed, due to slip along a
<figure><img src="image_3.png" /><figcaption>Figure 1</figcaption></figure>

Several fault zones in various geological environments show morphological corrugations at all scales. a Strike-slip plane in silt rocks of the Corona Heights Fault (San Francisco, CA, USA) showing corrugations at the meter scale. b Zoom on the slip surface, showing striations and asperities at the centimeter scale (pencil for scale). c Strike-slip plane of the La Cléry Fault (Vercors, France), in limestones, show large wavelengths topographical variations (hammer for scale). d Red clay layer drag along the Taxarchion Fault (Corinth, Greece) smear this normal fault plane

well-defined direction. We have extracted several hundred profiles in this direction of slip to analyze the statistical properties of the fault roughness (Fig. 2c). A Fourier power spectrum method was used, which has been shown to be robust and reliable, to characterize scaling properties of surfaces (CANDELA et al., 2009). For all 2D DEM data, the topography of the slip surfaces shows a self-affine geometry, demonstrated by a linear relationship when plotting on log-log axes the Fourier power spectrum versus the spatial wavenumber (Fig. 2d). The different slip surfaces analyzed cover approximately 6 decades of length scales and show a scaling relationship with a Hurst exponent $H_0$ close to 0.6 for profiles along the slip direction. These data confirm and extend

previous studies of scaling properties of fault surfaces (POWER et al., 1987; RENARD et al., 2006).

### 3. Visco-Elastic Shear Flow Between Rough Walls

#### 3.1. Setup

We now consider a simple model of a fault consisting of two undeformable plates separated by a fault gouge or fluid layer of thickness $H$ as sketched in Fig. 3. Either one or both of the plates are assumed to have a rough surface and the middle layer is assumed to satisfy a Maxwell visco-elastic rheology. Using this system, we analyze how the rough walls
<figure><img src="image_4.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_7.png" /></figure>
<caption>Figure 2</caption>

Surface roughness results of the Vuache-Sillingy strike-slip fault (French Alps). a The fault surface consists of many discrete slip surfaces at all scales separating lenses of variably deformed fault rock. The back rectangle corresponds to the surface shown on b. b Examples of Digital Elevation Model (DEM) at the outcrop scale (LiDAR) and at the laboratory scale (laser profilometer, white light interferometer). c Representative 1D self-affine profiles of the slip surface extracted and detrended from the DEM (b) along the direction of slip. Magnified portion of the profiles at the LiDAR scale (up) has a statistically similar appearance to the entire profile when using a self-affine transformation with a Hurst exponent equal to 0.6. d Fourier power spectrum calculated for the fault surface along slip, covering six orders of magnitude of spatial wavelengths. Power-law fit (thick gray line) with a prescribed roughness exponent (H = 0.6), connecting the field and laboratory data, is shown on plot for eye guidance. The inset displays an example of the height elevation Z (y-axis) versus wavelength (x-axis) of a rough profile. Contours (black dotted line) representing constant elevation Z to wavelength ratio, reflecting a self-similar behavior, are provided to allow easier interpretation of the spectra. Black arrows (at the bending of spectra) indicate the level of noise of the LiDAR and the lower limit for the fit performed at the WLI scale

change the mean shear stress when the top plate moves with a fixed velocity. In other words, we study an effective rheology of a system where the roughness of the outer plates couple to the flow properties of the full system. In general, assuming that the average separation distance *H* is fixed, we find that, for a small amplitude roughness of the outer surfaces, the viscous dissipation increases while the flow rate decreases with the amplitude.

The governing equations are given by the mass and momentum conservation laws supplemented by the rheological equation of state. Moreover, the

<figure><img src="image_8.png" /><figcaption>Figure 3</figcaption></figure>

Model system of a visco-elastic layer trapped between two rigid planes. The lower plane is assumed to have a rough surface and a simple shear is applied to the system by moving the upper plane at a constant velocity. This system is utilized in an analysis of the effective rheology of a fault system

visco-elastic layer is assumed to be incompressible with constant density, which implies that

$$ \partial_i v_i = 0, \qquad (1) $$

and the momentum conservation in the limit where fluid inertia can be neglected is satisfied by the relation

$$ \partial_i T_{ij} = 0, \qquad (2) $$

where $T_{ij}$ is the fluid stress tensor. The stress components can be decomposed into a homogeneous isotropic part and deviatoric components $T_{ij} = -p\delta_{ij} + \tau_{ij}$, where the deviatoric components $\tau_{ij}$ can be related to the strain rate by the Oldroyd-type equations of state (OLDROYD, 1950; SHANKAR and KUMAR, 2004) for a Maxwell visco-elastic fluid

$$ \frac{\mu}{G} (\partial_t \tau_{ij} + v_k \partial_k \tau_{ij} - \partial_k v_i \tau_{kj} - \partial_k v_j \tau_{ki}) + \tau_{ij} = \mu (\partial_i v_j + \partial_j v_i), \qquad (3) $$

where *G* is the elastic shear modulus and *μ* is the Newtonian viscosity. These equations are brought in a non-dimensional form by rescaling the spacial coordinates in units of the layer thickness *H*, the velocity in units of the upper boundary velocity *V*, time and stresses are given in units of *H/V*, *μV/H*, respectively. The rough bottom surface is positioned in dimensionless units at *z*<sub>b</sub>(*x*). With these rescalings, Eq. 3 is equivalent to

$$ W (\partial_t \tau_{ij} + v_k \partial_k \tau_{ij} - \partial_k v_i \tau_{kj} - \partial_k v_j \tau_{ki}) + \tau_{ij} = \partial_i v_j + \partial_j v_i, \qquad (4) $$

where the variables now are dimensionless and the $W = \frac{\mu V}{HG}$ is the Weissenberg number representing the ratio between the stress relaxation time and convec- tive timescale. Newtonian rheology is obtained in the limit of instant stress relaxation, i.e. $W = 0$. To model the relative slip between plates, we assume that the top plate moves at a constant velocity, while the bottom plate is fixed, namely

$$ v_x(x, 1) = 1, \quad v_z(x, 1) = 0 \qquad (5) $$

$$ v_x(x, z_b) = 0, \quad v_z(x, z_b) = 0. \qquad (6) $$

In the case where the confining surfaces are flat, i.e. $z_b = 0$, the steady state solution of the flow is obtained by matching the boundary conditions, and thus $v_x^{(0)} = z$, $v_z^{(0)} = 0$. Here we have introduced an
upper index which refers to the solutions around a flat interface. If we now insert the fluid velocity into the equations of state, Eq. 4, we determine the stress components, i.e. $\tau_{xz}^{(0)} = 1$, $\tau_{xx}^{(0)} = 2W$ and $\tau_{zz}^{(0)} = 0$. Notice the jump in the normal stress components $\tau_{xx}^{(0)} - \tau_{zz}^{(0)}$, which vanishes in the Newtonian limit.

## 3.2. Small Perturbations to a Flat Fault Interface

While the flow field readily follows when the interfaces are flat, the calculation for a system with rough interfaces becomes more involved. We now derive expressions for the steady state dynamics and the corresponding stress state in the limit where the interfaces are roughened by small amplitude perturbations. We restrict ourselves to consider perturbations of the lower interface only, which is assumed to have a height profile in the z-direction given by the expression $z_b = \sum_n A_n e^{inqx} + c.c.$, where $A_n$ is the amplitude (in units of H), $n$ is the Fourier mode and $q = 2\pi/\lambda$ is the characteristic dimensionless wavenumber. The linear regime is set by the condition $nqA_n \ll 1$ for all $n$. In particular, we consider a sinusoidal perturbation, for which $A_1 = \text{const}$ and $A_n = 0$ for $n > 1$.

For small morphological amplitudes $|A_n| \ll 1$, all the relevant fields (velocity and stress) are expanded around the unperturbed state as:

$$F(x, z) = F^{(0)}(x, z) + \sum_{n=1}^{\infty} (F^{(n)}(z)e^{inqx} + c.c.), \quad (7)$$

where the function $F(x, z)$ is introduced as a generic expression for the variable under consideration. The disturbance field is decomposed into Fourier modes denoted by $F^{(n)}(z)$, which are determined from the linearized governing equations. After substituting the Fourier modes for stress and velocity into Eqs. 1, 2 and 4, the subsequent equations can be reduced to a single equation in $v_z^{(n)}(z)$ which has a general solution given by (GORODTSOV and LEONOV, 1967; SHANKAR and KUMAR, 2004)

$$v_z^{(n)}(z) = B_1 z e^{nqz} + B_2 z e^{-nqz} + B_3 e^{nq(-iW + \sqrt{1+W^2})z} \\ + B_4 e^{-nq(iW + \sqrt{1+W^2})z}. \quad (8)$$

The other component $v_x^{(n)}(z)$ follows directly from Eq. 1, namely

$$v_x^{(n)}(z) = \frac{i}{qn} \frac{dv_z^{(n)}(z)}{dz}. \quad (9)$$

The coefficients $B_k$, with $k = 1, \dots, 4$, are obtained by inserting the expressions for $v_x^{(n)}(z)$ and $v_z^{(n)}(z)$ into the boundary conditions from Eqs. 5–6. To the first order in the surface amplitude, the velocity modes satisfy the following boundary conditions

$$v_x^{(n)}(0) + A_n = 0, \quad v_z^{(n)}(0) = 0 \quad (10)$$

$$v_x^{(n)}(1) = 0, \quad v_z^{(n)}(1) = 0, \quad (11)$$

where the nontrivial equation for $v_x^{(n)}$ follows from an expansion at the perturbed interface,

$$v_x(x, z_b) = v_x^{(0)}(x, 0) + \sum_n \left[ \partial_z v_x^{(0)}(x, 0) A_n + v_x^{(n)}(0) + \mathcal{O}(A_n^2) \right] e^{inqx} \\ + c.c., \quad (12)$$

and using the planar Couette solution, $\partial_z v_x^{(0)}(x, z) = 1$. Solving the system from Eqs. 10–11, we determine the coefficients $B_k$ as a function of amplitudes $A_n$, Weissenberg number $W$ and characteristic wave-number $q$.

## 3.3. Effective Flow and Energy Dissipation

In this section, we study the relation between wall roughness, mean flow rate and bulk energy dissipation. It is determined by the strain rate field given as the gradient of the velocity field, $e_{ij} = \partial_j v_i + \partial_i v_j$, where $\{i, j\} := \{x, y\}$. The shear rate, for instance, is reconstructed from the velocity Fourier modes and assumes the form

$$e_{xz}(x, z) = 1 + 2 \sum_n \Re \left\{ [\partial_z v_x^{(n)}(z) + inq v_z^{(n)}(z)] e^{inqx} \right\}. \quad (13)$$

Similar expressions apply for the other strain rate components $e_{zz}$ and $e_{xx}$. In general, the above expression depends on the detailed shape of the rough surface; however, we shall here consider an interface described by a single-mode profile $z_b = 2A_1 \cos(qx)$. Then, the mean flow rate is
<figure><img src="image_9.png" /><figcaption>Figure 4
Viscous shear flow along a lower rough surface. Left panel shows the flow profile in the case of a single harmonic perturbation. Right panel corresponds to a shear flow past a random surface with a Hurst exponent $H = 0.6$. The color scale represents the magnitude of the dimensionless shear $\tau_{xz}$ and the velocity field is shown in the vicinity of the rough surface</figcaption></figure>

obtained by the area average of the shear strain rate and given by

$$ \langle e_{xz} \rangle = 1 + \frac{q}{2\pi} \int_{0}^{2\pi/q} dx \int_{z_b}^{1} dz e_{xz}(x, z) \\ \approx 1 + f(q, W) A_1^2, \quad (14) $$

where $f(q, W)$ appears from the lowest order Taylor expansion in the roughness amplitude $A_1$. Similarly the mean energy dissipation rate in the bulk, $|\dot{E}|$, is computed as the area integral of $|\sum_{i,j} \partial_i v_j|^2$, thus

<figure><img src="image_10.png" /><figcaption>Figure 5
Mean shear stress as a function of the roughness amplitude in the viscous limit $W = 0$ and for a harmonic perturbation with $q = 1$. Data from numerical simulations are shown by red dots and the black curve correspond to the theoretical value obtained from the linear analysis, Eq. 20. It is observed that the linear analysis breaks down for a ratio between roughness amplitudes and layer thickness around 0.1</figcaption></figure>

$$ |\dot{E}| = \frac{1}{2} \left[ 1 + \frac{q}{2\pi} \int_{0}^{2\pi/q} dx \int_{0}^{1} dz (2e_{xz}^2 + e_{xx}^2 + e_{zz}^2) \right], \quad (15) $$

which to the lowest order in $A_1$ becomes

$$ |\dot{E}| \approx \frac{1}{2} [1 + g(q, W) A_1^2]. \quad (16) $$

The bulk energy dissipation depends both on the fluid rheology and the surface roughness. Although the function $g(q, W)$ follows directly from the linear expansion, it cannot be represented in a simple and short form.

### 3.4. Newtonian Limit

Relatively simple expressions can be obtained in the Newtonian limit, i.e. when $W = 0$. From Eq. 4, we observe that, in this case, the strain rate is the same as the deviatoric stress $\tau_{ij}$. In particular, the shear strain rate equals to

$$ e_{xz}(x, z) = 1 - \frac{4qA_1 \cos(qx)}{\cosh(2q) - 1 - 2q^2} [q(2-z) \cosh(qz) \\ + (2q^2z - 2q^2 - 1) \sinh(qz) \\ + qz \cosh(qz - 2q) + \sinh(qz - 2q)]. \quad (17) $$

Straightforwardly, we can determine the shear drag by evaluating the above expression on top and bottom surfaces, with the final result given by

$$ \tau_{xz}(x, 1) = 1 - 8qA_1 \frac{q \cosh(q) - \sinh(q)}{\cosh(2q) - 1 - 2q^2} \cos(qx) \quad (18) $$

$$ \tau_{xz}(x, z_b) = 1 + 4qA_1 \frac{\sinh(q) - 2q}{\cosh(2q) - 1 - 2q^2} \cos(qx). \quad (19) $$

The wall drag varies linearly with the roughness amplitude and frequency and, in a one mode approximation, alternates between regions of maximum and minimum resistance. These regions are located oppositely on the two plates. Namely, the maximum drag is set at protrusion of the rough bottom surface, while the same point on the top flat surface corresponds to a minimum drag.
<figure><img src="image_11.png" /><figcaption>Figure 6
Left panel $g(q,0)$ showing the increase in the mean viscous dissipation with the amplitude and wavenumber. Right panel The coefficient $f(W, q)$ from Eq. 14 as a function of wavenumber $W$ and $q$</figcaption></figure>

Also, inserting Eq. 17 into Eq. 14, we find that the effective shear flow depends on the amplitude to the lowest order,

$$ \langle e_{xz} \rangle = 1 - 4q \frac{\sinh(2q) - 2q}{\cosh(2q) - 1 - 2q^2} A_1^2. \quad (20) $$

We notice that in Eq. 20, the mean flow rate decreases monotonically with amplitude $A_1$ and wavenumber $q$ and at some point will start having negative values corresponding to a change in the flow direction. This means that the linear prediction will break down at finite amplitudes and one would expect that the mean shear will reach a minimum saturation value. To verify this, we resort to numerical simulations to study the steady state flow properties around a surface with a larger roughness amplitude.

A decrease in the mean flow rate is associated with an increase of the mean viscous dissipation, when $g(q, 0) > 0$ as seen in Fig. 6.

### 3.5. Numerical Approach

In order to simplify the analysis, we shall here assume that the embedded layer in our model configuration effectively behaves like a viscous fluid. That is, we consider the case where $W = 0$ for which the constitutive equations, Eqs. 1–4, reduce to those of an incompressible Stokes flow. We discretized the equations using a Galerkin finite element scheme on an unstructured triangular mesh. In general, the incompressibility condition is difficult to tackle

numerically, since it leads to a singular matrix in the pressure equation. One way to resolve this problem is to use a mixed finite element formulation with quadratic velocity shape functions and a discontinuous linear interpolation for the pressure degrees of freedom. Our numerical implementation follows the algorithm proposed in (DABROWSKI et al., 2008).

A snapshot of the shear stress field for a finite amplitude roughness is illustrated in Fig. 4. We observe that the presence of a rough interface locally disturbs the flow profile of the viscous layer. This in turn, will lead to changes in the overall shear resistance of the system and thereby increase the energy dissipation in the fault.

A comparison between the theoretical prediction given in Eq. 20 and the numerically computed mean shear stress is shown in Fig. 5. In general, the increase in shear resistance becomes less pronounced as the roughness amplitude is increased. It might be speculated that in a real system the resistance might even start to decrease for large roughness amplitudes, since vortices may form and remain trapped in valleys of the surface morphology, e.g. (SKJETNE et al., 1999).

### 3.6. Effect of Fluid Elasticity

For a finite $W$ number, we notice that, compared to the Newtonian limit, there is an extra contribution to the strain rate due to stress relaxation, as shown in
Eq. 4. Thus, measuring the shear stress is not the same as measuring the strain rate. That being said, the mean flow rate also decreases with the amplitude for $W > 0$, as shown in Fig. 6, albeit at a slower rate compared to the Newtonian limit. The fact that the shear flow rate at a given roughness is larger for a visco-elastic fluid than for a purely viscous one can be related to a positive contribution of the stress relaxation as suggested by Fig. 6.

## 4. Concluding Remarks

Recent measurements of fault slip surfaces have revealed a morphological roughness that spans a wide range of scales. For most active faults, it is likely that the roughness may develop on a time scale different than that introduced by the shear rate. The exact mechanisms leading to rough interfaces are largely unknown and may be related to both damage accumulation and recovery. On short time scales, rapid rupture can damage wall rocks and produce abrasive wear, while on a longer time scale, the branching of fractures or healing processes may corrugate the slip surface.

Here, we have considered the implications of this roughness on the dynamics of a fault by introducing a simple model consisting of a visco-elastic layer sandwiched between two rigid plates. The main result of our analysis is that the mean shear flow rate decreases with increasing roughness amplitude, while at the same time there is an increase in the mean bulk viscous dissipation. In the Newtonian limit, the maximum shear flow rate is attained when the interface is flat and gradually decreases as the roughness amplitude is increased. Numerical simulations at finite amplitude suggest that the mean flow may approach a minimum value which is independent of the amplitude. At a finite stress relaxation, mean flow rate is relatively higher than in the viscous limit, suggesting that the elastic modes enhance the total deformation rate.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

## REFERENCES

L. ANGHELUTA, E. JETTESTUEN, J. MATHIESEN, F. RENARD, and JAMTVEIT. *B. Stress-driven phase transformation and the roughening of solid-solid interfaces*. Physical Review Letters, 100:096106, 2008.

L. ANGHELUTA, J. MATHIESEN, C. MISBAH, and F. RENARD. *Morphological instabilities of stressed and reactive geological interfaces*. Journal of Geophysical Research, 115: doi: 10.1029/2009JB006880, 2010.

F. ARTHAUD and M. MATTAUER. *Exemples de stylolites d'origine tectonique dans le Languedoc, leur relation avec la tectonique cassante*. Bulletin de la Société Géologique de France, 11:738-744, 1969.

J.S. CAINE, J.P. EVANS, and C.B. FORSTER. *Fault zone architecture and permeability structure*. Journal of Structural Geology, 24:1025-1028, 1996.

T. CANDELA, F. RENARD, M. BOUCHON, A. BROUSTE, D. MARSAN, J. SCHMITTSBUHL, and C VOISIN. *Characterization of fault roughness at various scales: Implications of three-dimensional high resolution topography measurements*. Pure and Applied Geophysics, 166: doi:10.1007/s00024-009-0521-2, 2009.

J. CHERY, M. D. ZOBACK, and S. HICKMAN. *A mechanical model of the san andreas fault and safod pilot hole stress measurements*. Geophysical Research Letters, 31:L15S13, doi:10.1029/2004GL019521, 2004.

F.M. CHESTER, J.S. CHESTER, D.L. KIRSCHNER, S.E. SCHULZ, and J.P. EVANS. *Rheology and deformation of the lithosphere at continental margins*, chapter Structure of large-displacement, strike-slip fault zones in the brittle continental crust, page 223. Columbia University Press, New York, 2004.

F.M. CHESTER and J.M. LOGAN. *Implications for mechanical-proPERTIES of brittle faults from observations of the punchbowl fault zone, California*. Pure and Applied Geophysics, 124:79-106, 1986.

F.M. CHESTER and J.S. CHESTER. *Stress and deformation along wavy frictional faults*. Journal of Geophysical Research, 105(B10): 23421, 2000.

F. CUISIAT and S. SKURTVEIT. An experimental investigation of the development and permeability of clay smears along faults in uncemented sediments. *Journal of Structural Geology*, pages 1-14, doi:10.1016/j.jsg.2009.12.005, 2010.

M. DABROWSKI, M. KROTKIEWSKI, and D.W. SCHMID. *Milamin: Matlab-based finite element method solver for large problems*. Geochem. Geophys. Geosyst., 9:Q04030, doi:10.1029/2007GC001719, 2008.

D.L. EGHOLM, O.R. CLAUSEN, M. SANDIFORD, KIRSTENSEN M. B., and J.A. KORSTGARD. *The mechanics of clay smearing along faults*. Geology, 36:7877908, 2008.

V.A. GORODTSOV and A.I. LEONOV. *On a linear instability of a plane parallel couette flow of viscoelastic fluid*. Journal of Applied Mathematics and Mechanics, 31:310-319, 1967.

B. KAMB. *Sliding motion of glaciers: theory and observation*. Reviews of Geophysics, 8(4):673-728, 1970.

Y. KLINGER. *Relation between continental strike-slip earthquake segmentation and thickness of the crust*. Journal of Geophysical Research, 114, 2010.
J.G. OLDROYD. *On the formulation of rheological equations of state.*
Proceedings of the Royal Society of London Series A - Mathe-
matical and Physical Sciences, 200:523-541, 1950.

W.-L. POWER, T.-E. TULLIS, S.-R. BROWN, G.-N. BOITNOTT, and C.-H. SCHOLZ. *Roughness of natural fault surfaces.* Geophysical Research Letters, 14:29-32, 1987.

W.-L. POWER, T.-E. TULLIS, and J.-D. WEEKS. *Roughness and wear
during brittle faulting.* Journal of Geophysical Research,
93:15268-15278, 1988.

M.E. PRITCHARD and M. SIMONS. *An aseismic slip pulse in northern
chile and along-strike variations in seismogenic behavior.*
Journal of Geophysical Research, 111:B08405, doi:10.1029/2006
JB004258, 2006.

L.B. RAILSBACK and L.M. ANDREWS. *Tectonic stylolites in the
undeformed Cumberland Plateau of southern Tennessee.* Journal
of Structural Geology, 17:911-915, 1995.

F. RENARD, J.-P. GRATIER, and B. JAMTVEIT. *Kinetics of crack-seal-
ing, intergranular pressure solution, and compaction around
active faults.* Journal of Structural Geology, 22:1395-1407, 2000.

F. RENARD, C. VOISIN, D. MARSAN, and J. SCHMITTBUHL. *High res-
olution 3d laser scanner measurements of strike-slip faults*

quantify its morphological anisotropy at all scales. Geophysical
Research Letters, 33:L04305, 2006.

E.H. RUTTER. *The kinetics of rock deformation by pressure solution.*
Philosophical Transactions of the Royal Society of London,
283:203–219, 1976.

A. SAGY, E.E. BRODSKY, and G.J. AXEN. *Evolution of fault-surface roughness with slip.* Geology, 35:283-286, 2007.

F. SAUCIER, E. HUMPHREYS, and R. WELDON II. *Stress near geo-
metrically complex strike-slip faults: application to the San
Andreas Fault at Cajon Pass, southern California.* Journal of
Geophysical Research, 97(B4):5081–5094, 1992.

V. SHANKAR and S. KUMAR. *Instability of viscoelastic plane couette flow past a deformable wall.* Journal of Non-Newtonian fluid mechanics, 116:371-393, 2004.

Z. SHIPTON, A.M. SODEN, and KIRKPATRICK J. D. *Earthquakes: Radiated energy and the physics of faulting*, chapter How thick is a fault ? Fault displacement-thickness scaling revisited, page 193. American Geophysical Union, 2006.

E. SKJETNE, A. HANSEN, and J.S. GUDMUNDSSON. *High-velocity flow in a rough fracture.* Journal of Fluid Mechanics, 383:1-28, 1999.