Rocks and Water

• About
• Bio
• Research
• Publications
• Global Block Model
• Archive

Evolution of fault stress through the earthquake cycle: a Mohr circle perspective

Richard Styron

Mon 06 January 2014

In my recent work on topographic stresses and their interaction with faulting, I have been thinking a bit about what stresses look like on a fault during the earthquake cycle. I wrote this IPython notebook to explore this topic a little bit.

A note about reading my posts that are IPython notebooks: the formatting is better if the menubar on the right side of this blog is minimized, or the browser window is pretty wide. Otherwise, code is often wrapped strangely, and the images extend way off into the margins.

First let's import some modules.

%pylab inline
%config InlineBackend.figure_format = 'svg'
#config InlineBackend.figure_format = 'retina'

Populating the interactive namespace from numpy and matplotlib

import numpy as np
import matplotlib.animation as animation
from tempfile import NamedTemporaryFile
from IPython.display import HTML

The animation stuff here is based on Jake Vanderplas's example. Thanks, Jake!

VIDEO_TAG = """<video controls>
<source src="data:video/x-m4v;base64,{0}" type="video/mp4">
Your browser does not support the video tag.
</video>"""

def anim_to_html(anim):
    if not hasattr(anim, '_encoded_video'):
        with NamedTemporaryFile(suffix='.mp4') as f:
            anim.save(f.name, fps=20, dpi=100,
                      extra_args=['-vcodec', 'libx264',
                                  '-pix_fmt', 'yuv420p'])
            video = open(f.name, "rb").read()
        anim._encoded_video = video.encode("base64")
        
    return VIDEO_TAG.format(anim._encoded_video)


def display_animation(anim):
    plt.close(anim._fig)
    return HTML(anim_to_html( anim) )

animation.Animation._repr_html_ = anim_to_html

Mohr circles of the earthquake cycle

One fixed stress component

My first notions of this were that one of the principal stresses stays constant throughout the earthquake cycle, while the other principal stress varies$^1$. Once $\tau = \tau_f$, i.e. the shear stress $\tau$ reaches the failure shear stress $\tau_f$ for that particular normal stress (this is the Mohr-Coulomb failure envelope), there is an earthquake.

$^1$I am not sure this is correct, and Eric Hetland disagrees with me on this; there are some particular implications for normal stress drop during earthquakes that may be problematic. I have also implemented a model (at the end) where the normal stress stays constant. I also am not sure this is correct...

I'm also not going to provide a much of a description of Mohr circles here. Suffice it to say that they are a graphical and quantitative representation of the state of stress at a point, and the normal and shear stresses on a plane of any orientation are represented as points on the Mohr circle itself. The venerable Steven Dutch has a good page here for an introduction.

Thrust fault

Fault dips 30$^\circ$. Maximum principal stress $s_1$ is horizontal, and minimum principal stress $s_2$ is vertical.

During the course of the interseismic period, $s_1$ increases with time to failure; $s_2$ is constant in this model.

s2 = 20
s1 = np.array([20, 30, 40, 50, 60, 70, 80, 90, 100])

s_ = (s1 + s2) / 2.
r_ = (s1 - s2) / 2.

s_f = s_ + (r_ * np.cos(2 * np.pi / 3) )
t_f = r_ * np.sin(2 * np.pi / 3)

stress relationships

$\bar{s} = (s_1 + s_2) / 2 $ mean stress

$ r = (s_1 - s_2) / 2 $ radius of Mohr circle

$ \sigma_f = \bar{s} + r \ \cos 2 \theta $ normal stress at failure

$ \tau_f = r \ \sin 2 \theta $ shear stress at failure

fig = figure()

for i, s in enumerate(s_):
    fig.gca().add_artist( Circle( (s, 0), r_[i], fill=False, 
                                 color='b', alpha=0.4) )
    
    plot([s, s_f[i]], [0, t_f[i]], 'c--', alpha=0.4)
    

# 
plot(s_f, t_f, 'co-', label='stress on fault')

# plot failure envelope, coeff. fric. = 0.6
plot(np.array([0,120]), (np.tan( np.pi/6.)* np.array([0,120]) + 11.8), 
     'r', label=r'static/ failure friction $\mu_s$')

plot(np.array([0,120]), (np.tan( np.pi/8.)* np.array([0,120]) + 11.8),
     'orange', label='kinematic/ arrest friction $\mu_k$')

axvline(color='k')
axhline(color='k')
axis('equal')
xlabel(r'$\sigma$', fontsize=18)
ylabel(r'$\tau$', fontsize=18)
legend(loc='upper left')

show()

Mohr circles showing the stress at different points in the earthquake cycle. Cyan circles indicate $\tau$ and $\sigma$ on the fault plane; the solid cyan line connecting them represents the stress path the fault can take. The dashed lines go to the mean stress for each Mohr circle.

OK, so what does this show? In this plot, we see that as $s_1$ increases from its initial equality with $s_2$, all stresses not aligned with $s_2$ (which is roughly vertical for shallow thrust faulting) increase. So the horizontal $s_1$ and the differential stress (the diameter of the circle) increase, until the circle touches the red line, which is the Mohr-Coulomb failure envelope (the slope of which is $\mu_k$, the static fault friction). Note that $d \tau / d \sigma $ on the fault plane (represented by the cyan points and line) is much greater than $\mu$, so inevitably as $s_1$ increases, the $\tau / \sigma$ will hit the failure envelope.

Then, although it's not super apparent here, the fault starts slipping in an earthquake.

What happens during an earthquake is pretty complicated, and I won't try to represent it with a Mohr circle (maybe once I understand it better I will). But one thing that is certain is that $\mu_k$, the kinematic fault friction (orange line), is less than $\mu_s$, the static fault friction (red line). $\mu_k$ is almost certainly variable through the earthquake, as described by rate-state friction laws.

As the fault slips, $\tau$ is reduced as the stored elastic energy is used up in slip, heating, seismic radiation, etc. So the circle decreases, until it hits the orange line, below which the fault is stable again: $\tau$ can't overpower $\sigma$ anymore.

That is sweet! But, it's hard to see in the context of an 'earthquake cycle'.

So let's make a video that shows repeated earthquakes.

Thrust animation

Same problem as above.

#set up steps

s1a = range(20,101)
s1a.append(90)
s1a.extend( range(80,101) )
s1a.append(90)
s1a.extend( range(80,101) )
s1a.append(90)
s1a.extend( range(80,101) )

s1a = np.array(s1a)

# make animation
fig, ax = subplots()

# plot failure envelope, coeff. fric. = 0.6
plot(np.array([0,120]), (np.tan( np.pi/6.)* np.array([0,120]) 
                         + 11.8), 'r', alpha=1, label=r'$\mu_s = 0.6$')

plot(np.array([0,120]), (np.tan( np.pi/8.)* np.array([0,120]) 
                         + 11.8), 'orange', alpha=1, label=r'$\mu_k = 0.4$')
legend(loc='upper left')

axvline(color='k')
axhline(color='k')
axis('equal')
xlabel(r'$\sigma$', fontsize=18)
ylabel(r'$\tau$', fontsize=18)


mohr = Circle( (-10, 0), radius=0, fill=False,
                      color = 'b', alpha=1)

fail_pt, = ax.plot( [], [], 'co')

def init():
    fail_pt.set_data([], [])
    mohr.center = (20, 20)
    mohr.radius = 0
    ax.add_patch(mohr)

    return (mohr, fail_pt,)
    
def update(s1a):
    s2a = 20
    
    sa_ = (s1a + s2a) / 2.
    ra_ = (s1a - s2a) / 2.
    
    sa_fail = sa_ + ra_ * np.cos(2 * np.pi / 3.)
    ta_fail = ra_ * np.sin(2 * np.pi / 3.)
    
    fail_pt.set_data([sa_fail], [ta_fail])
    mohr.center = (sa_, 0)
    mohr.radius = (ra_)

    return (mohr, fail_pt,)

anim = animation.FuncAnimation(fig, update, init_func = init,
                               frames = s1a, blit=True)

#anim.save('fail_anim.mp4', dpi=200, fps=20)

#show()
display_animation(anim)

So hopefully this makes it a little more clear: stress increases until the Mohr circle reaches the failure envelope, then an earthquake happens and stress is reduced until the Mohr circle reaches the arrest envelope.

Now let's look at a normal faulting event.

Normal faulting

The setup is a little different here. Fault dips at 60 degrees. $s_1$ constant at 100, $s_2$ decreases with time(!) This was counter-intuitive to me, but it's important to remember that tectonic stress is tensional in normal faulting environments, but the lithostatic pressure is high enough that the faults are still under compression. See Mancktelow 2008 for a discussion on this; he refers to the situation as 'tectonic underpressure', because the mean normal stress (the pressure) is below lithostatic.

ns1 = 100
ns3 = np.array([100, 90, 80, 70, 60, 50, 40, 30, 20])

ns_ = (ns1 + ns3) / 2.
nr_ = (ns1 - ns3) / 2.

ns_f = ns_ + (nr_ * np.cos(2 * np.pi / 3) )
nt_f = -(nr_ * np.sin(2 * np.pi / 3) )

fig = figure()

for i, s in enumerate(ns_):
    fig.gca().add_artist( Circle( (s, 0), r_[i], fill=False, color='b',
                                 alpha=0.4) )
    plot([s, ns_f[i]], [0, nt_f[i]], 'c--', alpha=0.4)

plot(ns_f, nt_f, 'co-')

plot(np.array([0,120]), -(np.tan( np.pi/6.)* np.array([0,120])  
                         + 11.8), 'r')

plot(np.array([0,120]), -(np.tan( np.pi/8.)* np.array([0,120]) 
                         + 11.8), 'orange')

axvline(color='k')
axhline(color='k')
axis('equal')
xlabel(r'$\sigma_n$', fontsize=18)
ylabel(r'$\tau$', fontsize=18)
show()

Let's make another movie!

# set up steps, using s_2 as the steps in the animation

ns1a = 100

ns2a = range(20,101)[::-1]
ns2a.append(25)
ns2a.extend( range(20,31)[::-1])
ns2a.append(25)
ns2a.extend( range(20,31)[::-1])
ns2a.append(25)
ns2a.extend( range(20,31)[::-1])
ns2a.append(25)
ns2a.extend( range(20,31)[::-1])

ns2a = np.array(ns2a)

# make animation

fig, ax = subplots()

# plot failure envelope, coeff. fric. = 0.6
plot(np.array([0,120]), -(np.tan( np.pi/6.)* np.array([0,120])  
                         + 11.8), 'r', label=r'$\mu_s = 0.6$')

# plot arrest envelope, coeff. fric. = 0.4
plot(np.array([0,120]), -(np.tan( np.pi/8.)* np.array([0,120]) 
                         + 11.8), 'orange', label=r'$\mu_k = 0.4$')
legend(loc='lower left')

# axis and labels
axvline(color='k')
axhline(color='k')
axis('equal')
xlabel(r'$\sigma$', fontsize=18)
ylabel(r'$\tau$', fontsize=18)

# mohr cirlces
mohr = Circle( (-10, 0), radius=0, fill=False, color='b')

fail_pt, = ax.plot( [], [], 'co')

def init():
    fail_pt.set_data([], [])
    mohr.center = (100, 0)
    mohr.radius = 0
    ax.add_patch(mohr)
    
    return (mohr, fail_pt)


def update(ns2a):
    ns1a = 100
    
    nsa_ = (ns1a + ns2a) / 2.
    nra_ = (ns1a - ns2a) / 2.
    
    nsa_fail = nsa_ + (nra_ * np.cos(2 * np.pi / 3) )
    nta_fail = -(nra_ * np.sin(2 * np.pi / 3) )
    
    fail_pt.set_data([nsa_fail], [nta_fail])
    mohr.center = (nsa_, 0)
    mohr.radius = (nra_)
    
    return (mohr, fail_pt)

anim = animation.FuncAnimation(fig, update, init_func = init,
                               frames = ns2a, blit=True)

display_animation(anim)

Edifying.

To the extent that these toy models are correct, there are some interesting differences between the thrust and normal cases. One of them is that the stress drop is a lot less; $\Delta s_2$ in the normal fault model is about 10, which is half of $\Delta s_1$ in the thrust fault model. This means that for a given shear stress accumulation rate, earthquakes would be smaller and more frequent. However, continental normal faults typically have much lower slip rates than thrust faults (you heard it here first folks), so it may be difficult to find good datasets to compare this to in the field or literature.

Earlier I mentioned that Hetland has an issue with the stress evolution as depicted above. This problem is that in a nice double-couple earthquake, there is no reduction in normal stress during slip. But clearly, normal stress increases during the interseismic period, as shown by the path of the blue dot in these videos.

This leaves us with a handful of options:

1. Normal stress is actually constant during the interseismic period.
2. Normal stress decreases postseismically.
3. Normal stress actually does decrease coseismically, and we just don't see it in the radiated energy.

I am going to look at the first option. The second one seems great, but I don't know how to evaluate it. The third one is appealing to me, but the seismologists disagree.

Thrust faulting, constant $\sigma_n$

stress relationships

$\bar{s} = (s_1 + s_2) / 2 $

$ r = (s_1 - s_2) / 2 $

$ \sigma_n = \bar{s} + r \ \cos 2 \theta $

$ \tau = r \ \sin 2 \theta $

$ k = \cos 2 \theta $

$ 2 \sigma_n = (1 + k) \ s_1 + (1 - k) \ s_2 $

make some functions to help with the calcs

def get_s2(sf=None, s1=None, k=None, theta=None):
    if k == None:
        if not theta == None:
            k = np.cos( 2 * theta)
        else:
            raise Exception('need k or theta')
        
    return (2 * sf - (1 + k) * s1 ) / (1 - k)


def get_r(s1=None, s2=None):
    return (s1 - s2) / 2.


def get_s(s1=None, s2=None):
    return (s1 + s2) / 2.


def get_t_f(r=None, theta=None, s1=None, s2=None):
    if r == None:
        r = get_r(s1, s2)
    
    return r * np.sin( 2 * theta)

s_f = 40 # constant

theta = np.pi / 3.
k =  np.cos( 2 * theta)

s_1 = np.array([40, 50, 60, 70, 80, 90, 100])

s_2 = get_s2(sf = 40, s1 = s_1, k=k)

r_ = get_r(s1=s_1, s2=s_2)

t_f = get_t_f(r=r_, theta=theta)

s_ = get_s(s1 = s_1, s2 = s_2)

fig = figure()

for i, s in enumerate(s_):
    fig.gca().add_artist( Circle( (s, 0), r_[i], fill=False, 
                                 color='b', alpha=0.4) )
    
    plot([s, s_f], [0, t_f[i]], 'c--', alpha=0.4)

plot(np.ones(t_f.shape) * s_f, t_f, 'co-')

# plot failure envelope, coeff. fric. = 0.6
plot(np.array([0,120]), (np.tan( np.pi/6.)* np.array([0,120]) 
                         + 11.8), 'r', alpha=0.6)

plot(np.array([0,120]), (np.tan( np.pi/8.)* np.array([0,120]) 
                         + 11.8), 'orange', alpha=0.6)

axvline(color='k')
axhline(color='k')
axis('equal')
xlabel(r'$\sigma$', fontsize=18)
ylabel(r'$\tau$', fontsize=18)
show()

Looks a little funky to me.

Let's see how she moves.

#set up steps

s1a = range(40,101)
s1a.append(90)
s1a.extend( range(89,101) )
s1a.append(90)
s1a.extend( range(89,101) )
s1a.append(90)
s1a.extend( range(89,101) )

s1a = np.array(s1a)

# make animation
fig, ax = subplots()

# plot failure envelope, coeff. fric. = 0.6
plot(np.array([0,120]), (np.tan( np.pi/6.)* np.array([0,120]) 
                         + 11.8), 'r', alpha=1, label=r'$\mu_s = 0.6$')

plot(np.array([0,120]), (np.tan( np.pi/8.)* np.array([0,120]) 
                         + 11.8), 'orange', alpha=1, label=r'$\mu_k = 0.4$')
legend(loc='upper left')

axvline(color='k')
axhline(color='k')
axis('equal')
xlabel(r'$\sigma$', fontsize=18)
ylabel(r'$\tau$', fontsize=18)


mohr = Circle( (-10, 0), radius=0, fill=False,
                      color = 'b', alpha=1)

fail_pt, = ax.plot( [], [], 'co')

def init():
    fail_pt.set_data([], [])
    mohr.center = (40, 40)
    mohr.radius = 0
    ax.add_patch(mohr)

    return (mohr, fail_pt,)
    
def update(s1a):
    s2a = get_s2(sf = 40, s1=s1a, k=k)
    
    sa_ = get_s(s1 = s1a, s2=s2a)
    ra_ = get_r(s1 = s1a, s2=s2a)
    
    sa_fail = 40
    ta_fail = get_t_f(r = ra_, theta=theta)
    
    fail_pt.set_data([sa_fail], [ta_fail])
    mohr.center = (sa_, 0)
    mohr.radius = (ra_)

    return (mohr, fail_pt,)

anim = animation.FuncAnimation(fig, update, init_func = init,
                               frames = s1a, blit=True)

display_animation(anim)

OK, so this doesn't make a lot of sense to me. The big problem is that it's clear that $s_1$ and $s_2$ have to be related by a function that includes the fault. But these should be principal stresses and independent of the fault orientation (what if there are many faults in the immediate vicinity?). So I don't think that I believe this model.

But that still leaves us with a conundrum: How can normal stress increase during the interseismic period and not decrease coseismically?

Maybe it does decrease postseismically. And maybe I will look into that some day...

tectonics research python active faults earthquakes stress