6. Quasidynamic earthquake cycles: antiplane¶

6.1. Summary¶

Solving the SCEC SEAS problem BP1-QD for earthquake cycle simulation.
Computing traction given fault slip by rearranging the discretized integral equation.
Solving the quasidynamic rate-and-state friction equations with a protected Newton solver.
Time-stepping over several earthquake ruptures with an adaptive time-stepping routine.
Comparing against results from SCEC SEAS benchmark participants.

6.2. Setting up the geometry.¶

Previously, we built the infrastructure for solving for displacement and stress everywhere in the domain given slip on a fault in an antiplane half space. This time, we’ll go a step further and simulate rate and state friction and slip on that fault. We’ll consider the quasidynamic approximation where we use a zero-th order approximation to the stress changes caused by fault slip. We’ll be solving SCEC SEAS problem BP1-QD. Click the link!

For a bit more background on the ideas here, you could take a look at my previous explanation of modeling a quasidynamic spring-block-slider.

import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from tectosaur2 import (
    gauss_rule,
    refine_surfaces,
    integrate_term,
    panelize_symbolic_surface,
)
from tectosaur2.laplace2d import double_layer, hypersingular

To start, we’ll build a heavily refined fault mesh and surface mesh. The surface mesh will be quite long to make sure that we’re capturing the half space assumption.

surf_half_L = 10000000
fault_bottom = 40000
max_panel_length = 200
shear_modulus = 3.2e10

quad_rule = gauss_rule(5)
t = sp.var("t")

panel_bounds = np.stack((edges[:-1],edges[1:]), axis=1)

fault = panelize_symbolic_surface(
    t, t * 0, fault_bottom * (t + 1) * -0.5,
    quad_rule,
    n_panels=400
)

free = refine_surfaces(
    [
        (t, -t * surf_half_L, 0 * t)  # free surface
    ],
    (quad_rule),
    control_points = [
        # nearfield surface panels and fault panels will be limited to 200m
        # at 200m per panel, we have ~40m per solution node because the panels
        # have 5 nodes each
        (0, 0, 1.5 * fault_bottom, 200),
        # farfield panels will be limited to 200000 m per panel at most
        (0, 0, surf_half_L, 200000),
    ]
)
print(
    f"The free surface mesh has {free.n_panels} panels with a total of {free.n_pts} points."
)
print(
    f"The fault mesh has {fault.n_panels} panels with a total of {fault.n_pts} points."
)

The free surface mesh has 958 panels with a total of 4790 points.
The fault mesh has 400 panels with a total of 2000 points.

plt.plot(free.pts[:,0]/1000, free.pts[:,1]/1000, 'k-')
plt.plot(fault.pts[:,0]/1000, fault.pts[:,1]/1000, 'r-')
plt.xlabel(r'$x ~ \mathrm{(km)}$')
plt.ylabel(r'$y ~ \mathrm{(km)}$')
plt.axis('scaled')
plt.xlim([-100, 100])
plt.ylim([-80, 20])
plt.show()

And, to start off the integration, we’ll construct the operators necessary for solving for free surface displacement from fault slip.

singularities = np.array(
    [
        [-surf_half_L, 0],
        [surf_half_L, 0],
        [0, 0],
        [0, -fault_bottom],
    ]
)

(free_disp_to_free_disp, fault_slip_to_free_disp), report = integrate_term(
    double_layer, free.pts, free, fault, singularities=singularities, return_report=True
)
fault_slip_to_free_disp = fault_slip_to_free_disp[:, 0, :, 0]
free_disp_solve_mat = (
    np.eye(free_disp_to_free_disp.shape[0]) + free_disp_to_free_disp[:, 0, :, 0]
)

6.3. Solving for fault traction given fault slip¶

The previous section should have been familiar, but now we need to go a step further. The goal is to construct a single matrix that computes the fault tractions given an input fault slip distribution. This matrix will be used during each time step of our fault simulation.

But because we’re using fullspace Green’s function, this matrix needs to include both the direct fault slip effect and the indirect effect that is mediated through the free surface. So, we’re actually going to be composing several matrices together to get our final total_fault_slip_to_fault_stress matrix:

Combining fault_slip_to_free_disp and free_disp_solve_mat, we can compute the surface displacement given a fault slip distribution.
We will construct free_disp_to_fault_stress for computing the fault stress caused by a surface displacement field.
We also construct the direct fault slip effect matrix, fault_slip_to_fault_stress.

Renaming these wordy matrix names to some simple letters:

A = fault_slip_to_fault_traction
B = free_disp_to_fault_traction
C = fault_slip_to_free_disp
D = free_disp_solve_mat

we can write the full system:

(6.1)¶\[\begin{bmatrix} A & B \\ C & D \end{bmatrix}\]

The second equation expands to the same linear system that we used before to solve for surface displacement: \(u = -D^{-1}Cs\). Plugging that in, we get:

(6.2)¶\[\begin{equation} t = As - BD^{-1}Cs = (A - BD^{-1}C)s \end{equation}\]

Thus, our goal matrix operator is \(A - BD^{-1}C\) and we will name this total_fault_slip_to_fault_traction. You might recognize this as the Schur complement of the matrix system above. We performed one of step of block matrix Gaussian elimination!

(free_disp_to_fault_stress, fault_slip_to_fault_stress), report = integrate_term(
    hypersingular,
    fault.pts,
    free,
    fault,
    tol=1e-12,
    singularities=singularities,
    return_report=True,
)
fault_slip_to_fault_stress *= shear_modulus
free_disp_to_fault_stress *= shear_modulus

fault_slip_to_fault_traction = np.sum(
    fault_slip_to_fault_stress[:,:,:,0] * fault.normals[:, :, None], axis=1
)
free_disp_to_fault_traction = np.sum(
    free_disp_to_fault_stress[:,:,:,0] * fault.normals[:, :, None], axis=1
)

/Users/tbent/Dropbox/active/eq/tectosaur2/tectosaur2/integrate.py:201: UserWarning: Some integrals failed to converge during adaptive integration. This an indication of a problem in either the integration or the problem formulation.
  warnings.warn(
/Users/tbent/Dropbox/active/eq/tectosaur2/tectosaur2/integrate.py:208: UserWarning: Some expanded integrals reached maximum expansion order. These integrals may be inaccurate.
  warnings.warn(

A = -fault_slip_to_fault_traction
B = -free_disp_to_fault_traction
C = fault_slip_to_free_disp
Dinv = np.linalg.inv(free_disp_solve_mat)

total_fault_slip_to_fault_traction = A - B.dot(Dinv.dot(C))

I’ll just do a quick check to confirm that this operator is accurate:

analytical_sd = -np.arctan(-fault_bottom / free.pts[:,0]) / np.pi
numerical_sd = Dinv.dot(-C).sum(axis=1)
rp = (fault.pts[:,1] + 1) ** 2
analytical_ft = -(shear_modulus / (2 * np.pi)) * ((1.0 / (fault.pts[:,1] + fault_bottom)) - (1.0 / (fault.pts[:,1] - fault_bottom)))
numerical_ft = total_fault_slip_to_fault_traction.sum(axis=1)
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.title('Stress error')
plt.plot(fault.pts[:,1], np.log10(np.abs((numerical_ft - analytical_ft)/(analytical_ft))), 'b-')
plt.xlim([-fault_bottom, 0])
plt.xlabel('$\mathrm{depth ~~ (m)}$')
plt.ylabel('$\log_{10}(|E_{\sigma}|)$')
plt.subplot(1,2,2)
plt.title('Displacement error')
plt.plot(free.pts[:,0], np.log10(np.abs(numerical_sd - analytical_sd)), 'b-')
plt.xlim([-100000, 100000])
plt.xlabel('$x ~~ \mathrm{(m)}$')
plt.ylabel('$\log_{10}(|E_u|)$')
plt.show()

/var/folders/mt/cmys2v_143q1kpcrdt5wcdyr0000gn/T/ipykernel_79373/3764098280.py:15: RuntimeWarning: divide by zero encountered in log10
  plt.plot(free.pts[:,0], np.log10(np.abs(numerical_sd - analytical_sd)), 'b-')

6.4. Rate and state friction¶

Okay, now that we’ve constructed the necessary boundary integral operators, we get to move on to describing the frictional behavior on the fault.

This isn’t intended to be an introducion to earthquake physics or rate-and-state friction. Instead, we’re focusing on the numerical implementation. If you’d like to learn more about rate-and-state friction, I’d recommend… (SUGGEST A BOOK OR PAPER HERE)

siay = 31556952 # seconds in a year
density = 2670  # rock density (kg/m^3)
cs = np.sqrt(shear_modulus / density)  # Shear wave speed (m/s)
eta = shear_modulus / (2 * cs)  # The radiation damping coefficient (kg / (m^2 * s))
Vp = 1e-9  # Rate of plate motion
sigma_n = 50e6  # Normal stress (Pa)

# parameters describing "a", the coefficient of the direct velocity strengthening effect
a0 = 0.01
amax = 0.025
fy = fault.pts[:, 1]
H = 15000
h = 3000
fp_a = np.where(
    fy > -H, a0, np.where(fy > -(H + h), a0 + (amax - a0) * (fy + H) / -h, amax)
)

# state-based velocity weakening effect
fp_b = 0.015

# state evolution length scale (m)
fp_Dc = 0.008

# baseline coefficient of friction
fp_f0 = 0.6

# if V = V0 and at steady state, then f = f0, units are m/s
fp_V0 = 1e-6

plt.figure(figsize=(3, 5))
plt.plot(fp_a, fy/1000, label='a')
plt.plot(np.full(fy.shape[0], fp_b), fy/1000, label='b')
plt.xlim([0, 0.03])
plt.ylabel('depth')
plt.legend()
plt.show()

First, a Python implementation of the quasidynamic shear stress equations. There are two components to a standard definition of rate and state friction:

(6.3)¶\[\begin{align} \tau_{qd} = F(V, \theta)\\ \frac{d\theta}{dt} = G(V, \theta) \end{align}\]

where \(F(V, \theta)\) is the strength of friction, depending on slip rate and state and \(G(V, \theta)\) gives the state evolution. For \(G\), we will use the aging law, given by the aging_law function below.

We also need to relate the quasidynamic shear stress to the quasistatic shear stress:

(6.4)¶\[\begin{equation} \tau_{qd} = \tau_{qs} - \eta * V \end{equation}\]

which results in the full quasidynamic stress balance equation:

(6.5)¶\[\begin{equation} \tau_{qs} - \eta * V - F(V, \theta) = 0 \end{equation}\]

The specific form we use here is implemented below in the qd_equation function.

At each time step, we’ll need to compute the velocity, \(V\) and the state derivative \(d\theta/dt\). Computing the state derivative is easy because it’s directly given by the aging law. But, computing the velocity is harder because its only given implicitly in the frictional strength equation. The solve_friction function below implements Newton’s method with a backtracking line search to solve the friction equation for the current velocity. The qd_equation_dV function implements the derivative of the stress balance with respect to velocity. This derivative is necessary for Newton’s method.

# The regularized form of the aging law. 
def aging_law(fp_a, V, state):
    return (fp_b * fp_V0 / fp_Dc) * (np.exp((fp_f0 - state) / fp_b) - (V / fp_V0))

def qd_equation(fp_a, shear_stress, V, state):
    # The regularized rate and state friction equation
    F = sigma_n * fp_a * np.arcsinh(V / (2 * fp_V0) * np.exp(state / fp_a))
    
    # The full shear stress balance:
    return shear_stress - eta * V - F

def qd_equation_dV(fp_a, V, state):
    # First calculate the derivative of the friction law with respect to velocity
    # This is helpful for equation solving using Newton's method
    expsa = np.exp(state / fp_a)
    Q = (V * expsa) / (2 * fp_V0)
    dFdV = fp_a * expsa * sigma_n / (2 * fp_V0 * np.sqrt(1 + Q * Q))

    # The derivative of the full shear stress balance.
    return -eta - dFdV


def solve_friction(fp_a, shear, V_old, state, tol=1e-10):
    V = V_old
    max_iter = 150
    for i in range(max_iter):
        # Newton's method step!
        f = qd_equation(fp_a, shear, V, state)
        dfdv = qd_equation_dV(fp_a, V, state)
        step = f / dfdv

        # We know that slip velocity should not be negative so any step that
        # would lead to negative slip velocities is "bad". In those cases, we
        # cut the step size in half iteratively until the step is no longer
        # bad. This is a backtracking line search.
        while True:
            bad = step > V
            if np.any(bad):
                step[bad] *= 0.5
            else:
                break

        # Take the step and check if everything has converged.
        Vn = V - step
        if np.max(np.abs(step) / Vn) < tol:
            break
        V = Vn
        if i == max_iter - 1:
            raise Exception("Failed to converge.")

    # Return the solution and the number of iterations required.
    return Vn, i

We’ll also check \(h^*\) which is the minimum length scale of an instability and \(L_b\), the length scale of the rupture process zone. Both these length scales need to be well resolved by the fault discretization. Here we have approximately eight point within the process zone and almost 80 points within \(h^*\).

mesh_L = np.max(np.abs(np.diff(fault.pts[:, 1])))
Lb = shear_modulus * fp_Dc / (sigma_n * fp_b)
hstar = (np.pi * shear_modulus * fp_Dc) / (sigma_n * (fp_b - fp_a))
mesh_L, Lb, np.min(hstar[hstar > 0])

(26.92346550528964, 341.3333333333333, 3216.990877275949)

6.5. Quasidynamic earthquake cycle derivatives¶

Let’s set of the last few pieces to do a full earthquake cycle simulation:

initial conditions.
a function to put the pieces together and calculate the full system state at each time step, including the time derivatives of slip and frictional state.
the time stepping algorithm itself.

First, initial conditions. This initial state is exactly as specified in the BP-1 description linked at the top:

We solve for the steady frictional state at each point using scipy.optimize.fsolve. This is the initial state.
We identify the value of shear stress that will result in steady plate rate slip rates in the deeper portions of the fault: tau_amax. This is the initial shear stress.
The initial slip deficit is zero.

from scipy.optimize import fsolve
import copy

init_state_scalar = fsolve(lambda S: aging_law(fp_a, Vp, S), 0.7)[0]
init_state = np.full(fault.n_pts, init_state_scalar)
tau_amax = -qd_equation(amax, 0, Vp, init_state_scalar)

init_traction = np.full(fault.n_pts, tau_amax)
init_slip_deficit = np.zeros(fault.n_pts)
init_conditions = np.concatenate((init_slip_deficit, init_state))

Next, solving for system state. This ties the pieces together by:

Solving for the quasistatic shear stress using the boundary integral matrices derived at the beginning of this notebook.
Solving the rate and state friction equations for the slip rate.
Calculating the state evolution using the aging law.

The middle lines of SystemState.calc do these three steps. There’s a bunch of other code surrounding those three lines in order to deal with invalid inputs and transform from slip to slip deficit.

class SystemState:

    V_old = np.full(fault.n_pts, Vp)
    state = None

    def calc(self, t, y, verbose=False):
        # Separate the slip_deficit and state sub components of the 
        # time integration state.
        slip_deficit = y[: init_slip_deficit.shape[0]]
        state = y[init_slip_deficit.shape[0] :]

        # If the state values are bad, then the adaptive integrator probably
        # took a bad step.
        if np.any((state < 0) | (state > 2.0)):
            return False

        # The big three lines solving for quasistatic shear stress, slip rate
        # and state evolution
        tau_qs = init_traction - total_fault_slip_to_fault_traction.dot(slip_deficit)
        V = solve_friction(fp_a, tau_qs, self.V_old, state)[0]
        dstatedt = aging_law(fp_a, V, state)
        if not np.all(np.isfinite(V)):
            return False
        self.V_old = V

        slip_deficit_rate = Vp - V
        out = slip_deficit, state, tau_qs, V, slip_deficit_rate, dstatedt
        self.state = out
        return self.state


def plot_system_state(t, SS, xlim=None):
    """This is just a helper function that creates some rough plots of the
    current state to help with debugging"""
    slip_deficit, state, tau_qs, V, slip_deficit_rate, dstatedt = SS
    
    slip = Vp * t - slip_deficit
    
    fd = -np.linalg.norm(fault.pts, axis=1)

    plt.figure(figsize=(15, 9))
    plt.suptitle(f't={t/siay}')
    plt.subplot(3, 3, 1)
    plt.title("slip")
    plt.plot(fd, slip)
    plt.xlim(xlim)

    plt.subplot(3, 3, 2)
    plt.title("slip deficit")
    plt.plot(fd, slip_deficit)
    plt.xlim(xlim)

    # plt.subplot(3, 3, 2)
    # plt.title("slip deficit rate")
    # plt.plot(fd, slip_deficit_rate)
    # plt.xlim(xlim)

    # plt.subplot(3, 3, 2)
    # plt.title("strength")
    # plt.plot(fd, tau_qs/sigma_qs)
    # plt.xlim(xlim)

    plt.subplot(3, 3, 3)
    plt.title("V")
    plt.plot(fd, V)
    plt.xlim(xlim)

    plt.subplot(3, 3, 5)
    plt.title(r"$\tau_{qs}$")
    plt.plot(fd, tau_qs)
    plt.xlim(xlim)

    plt.subplot(3, 3, 6)
    plt.title("state")
    plt.plot(fd, state)
    plt.xlim(xlim)

    plt.subplot(3, 3, 9)
    plt.title("dstatedt")
    plt.plot(fd, dstatedt)
    plt.xlim(xlim)
    plt.tight_layout()


    plt.show()


def calc_derivatives(state, t, y):
    """
    This helper function calculates the system state and then extracts the
    relevant derivatives that the integrator needs. It also intentionally
    returns infinite derivatives when the `y` vector provided by the integrator
    is invalid.
    """
    if not np.all(np.isfinite(y)):
        return np.inf * y
    state = state.calc(t, y)
    if not state:
        return np.inf * y
    derivatives = np.concatenate((state[-2], state[-1]))
    return derivatives

6.6. Integrating through time¶

%%time
from scipy.integrate import RK45

# We use a 5th order adaptive Runge Kutta method and pass the derivative function to it
# the relative tolerance will be 1e-11 to make sure that even 
state = SystemState()
derivs = lambda t, y: calc_derivatives(state, t, y)
integrator = RK45
atol = Vp * 1e-6
rtol = 1e-11
rk = integrator(derivs, 0, init_conditions, 1e50, atol=atol, rtol=rtol)

# Set the initial time step to one day.
rk.h_abs = 60 * 60 * 24

# Integrate for 1000 years.
max_T = 1000 * siay

n_steps = 1000000
t_history = [0]
y_history = [init_conditions.copy()]
for i in range(n_steps):
    
    # Take a time step and store the result
    if rk.step() != None:
        raise Exception("TIME STEPPING FAILED")
        break
    t_history.append(rk.t)
    y_history.append(rk.y.copy())

    # Print the time every 5000 steps
    if i % 5000 == 0:
        print(f"step={i}, time={rk.t / siay} yrs, step={(rk.t - t_history[-2]) / siay}")
    if rk.t > max_T:
        break

y_history = np.array(y_history)
t_history = np.array(t_history)

step=0, time=0.0021812680462764775 yrs, step=0.0021812680462764775
step=5000, time=140.7230427507717 yrs, step=0.02945013453871661
step=10000, time=196.52497632980683 yrs, step=6.699937178573698e-11
step=15000, time=196.5249766930108 yrs, step=7.13511577745264e-11
step=20000, time=196.52498488866365 yrs, step=7.762718868876156e-08
step=25000, time=272.17487700592767 yrs, step=1.3297123854634314e-10
step=30000, time=272.1748773997298 yrs, step=7.935965282334024e-11
step=35000, time=272.1748778087133 yrs, step=4.735589336411811e-11
step=40000, time=292.5353440506115 yrs, step=0.025529773554574715
step=45000, time=350.0284817740895 yrs, step=7.585404744348211e-11
step=50000, time=350.0284822570143 yrs, step=1.143552651498551e-10
step=55000, time=350.0284828072335 yrs, step=1.5738959326121706e-09
step=60000, time=428.3603641171298 yrs, step=1.396734966774442e-06
step=65000, time=428.3603895166313 yrs, step=8.38323217562627e-11
step=70000, time=428.36038998662326 yrs, step=6.908460257203191e-11
step=75000, time=432.77635516149826 yrs, step=0.023396389027961796
step=80000, time=506.74075025565014 yrs, step=6.93868099323645e-11
step=85000, time=506.74075072729175 yrs, step=1.3369653621114137e-10
step=90000, time=506.7407511239481 yrs, step=1.6760420204045886e-10
step=95000, time=571.1772883200007 yrs, step=0.03115928106939193
step=100000, time=585.1249321168854 yrs, step=8.195863612220059e-11
step=105000, time=585.1249326058276 yrs, step=7.277153236808961e-11
step=110000, time=585.1651243047899 yrs, step=0.00020757875702647396
step=115000, time=663.509382683838 yrs, step=7.482654241835127e-11
step=120000, time=663.5093831323696 yrs, step=1.1858616819451146e-10
step=125000, time=663.5093835144739 yrs, step=7.228800059155745e-11
step=130000, time=708.7786006395438 yrs, step=0.0306217858381618
step=135000, time=741.8938514583174 yrs, step=9.815695063602784e-11
step=140000, time=741.8938519540883 yrs, step=8.812366627298559e-11
step=145000, time=741.8944130341275 yrs, step=4.4345806513532365e-06
step=150000, time=820.2783209816516 yrs, step=6.962857582063059e-11
step=155000, time=820.2783214059152 yrs, step=1.1133319154652912e-10
step=160000, time=820.2783218379793 yrs, step=4.2550796334829805e-11
step=165000, time=846.4938094052974 yrs, step=0.026173168524190365
step=170000, time=898.6627909507023 yrs, step=7.45847765300852e-11
step=175000, time=898.662791441536 yrs, step=9.126662282044461e-11
step=180000, time=898.6627924229714 yrs, step=6.237439034320691e-09
step=185000, time=977.047260512177 yrs, step=1.2691500304527786e-09
step=190000, time=977.0472609600538 yrs, step=9.042044221151333e-11
step=195000, time=977.0472614254056 yrs, step=6.068323795478569e-11
step=200000, time=984.7989573452244 yrs, step=0.02885660199627129
CPU times: user 1h 28min 3s, sys: 28min 22s, total: 1h 56min 26s
Wall time: 1h 15min 56s

6.7. Plotting the results¶

Now that we’ve solved for 1000 years of fault slip evolution, let’s plot some of the results. I’ll start with a super simple plot of the maximum log slip rate over time.

derivs_history = np.diff(y_history, axis=0) / np.diff(t_history)[:, None]
max_vel = np.max(np.abs(derivs_history), axis=1)
plt.plot(t_history[1:] / siay, np.log10(max_vel))
plt.xlabel('$t ~~ \mathrm{(yrs)}$')
plt.ylabel('$\log_{10}(V)$')
plt.show()

And next, we’ll make the classic plot showing the spatial distribution of slip over time:

the blue lines show interseismic slip evolution and are plotted every fifteen years
the red lines show evolution during rupture every three seconds.

plt.figure(figsize=(10, 4))
last_plt_t = -1000
last_plt_slip = init_slip_deficit
event_times = []
for i in range(len(y_history) - 1):
    y = y_history[i]
    t = t_history[i]
    slip_deficit = y[: init_slip_deficit.shape[0]]
    should_plot = False

    # Plot a red line every three second if the slip rate is over 0.1 mm/s.
    if (
        max_vel[i] >= 0.0001 and t - last_plt_t > 3
    ):
        if len(event_times) == 0 or t - event_times[-1] > siay:
            event_times.append(t)
        should_plot = True
        color = "r"
    
    # Plot a blue line every fifteen years during the interseismic period
    if t - last_plt_t > 15 * siay:
        should_plot = True
        color = "b"

    if should_plot:
        # Convert from slip deficit to slip:
        slip = -slip_deficit + Vp * t
        plt.plot(slip, fy / 1000.0, color + "-", linewidth=0.5)
        last_plt_t = t
        last_plt_slip = slip
plt.xlim([0, np.max(last_plt_slip)])
plt.ylim([-40, 0])
plt.ylabel(r"$\textrm{z (km)}$")
plt.xlabel(r"$\textrm{slip (m)}$")
plt.tight_layout()
plt.savefig("halfspace.png", dpi=300)
plt.show()

And a plot of recurrence interval:

plt.title("Recurrence interval")
plt.plot(np.diff(event_times) / siay, "k-*")
plt.xticks(np.arange(0, 10, 1))
plt.yticks(np.arange(75, 80, 0.5))
plt.xlabel("Event number")
plt.ylabel("Time between events (yr)")
plt.show()

6.8. Comparison against SCEC SEAS results¶

Finally, we will compare against the results reported as part of the SEAS project. Here, I’ve grabbed the time series for a point at 7500m depth on the fault from results produced by Junle Jiang.

To compare, I grab the fault history for just the first event and then interpolate the values to a depth of exactly 7500m.

jiang_data = np.loadtxt("jiang7500.txt")
jiang_slip_rate = 10 ** jiang_data[:, 2]
jiang_stress = jiang_data[:, 3]

t_start_idx = np.argmax(max_vel > 1e-4)
t_end_idx = np.argmax(max_vel[t_start_idx:] < 1e-6)
n_steps = t_end_idx - t_start_idx
t_chunk = t_history[t_start_idx : t_end_idx]

shear_chunk = []
slip_rate_chunk = []
for i in range(n_steps):
    system_state = SystemState().calc(t_history[t_start_idx + i], y_history[t_start_idx + i])
    slip_deficit, state, tau_qs, V, slip_deficit_rate, dstatedt = system_state
    shear_chunk.append((tau_qs - eta * V))
    slip_rate_chunk.append(V)
shear_chunk = np.array(shear_chunk)
slip_rate_chunk = np.array(slip_rate_chunk)

fault_idx = np.argmax((-7480 > fy) & (fy > -7500))
VAvg = np.mean(slip_rate_chunk[:, fault_idx:(fault_idx+2)], axis=1)
SAvg = np.mean(shear_chunk[:, fault_idx:(fault_idx+2)], axis=1)

t_align = t_chunk[np.argmax(VAvg > 0.1)]
jiang_t_align = np.argmax(jiang_slip_rate > 0.1)

for lims in [(-1, 1), (-15, 30)]:
    plt.figure(figsize=(12, 8))
    plt.subplot(2, 1, 1)
    plt.plot(t_chunk - t_align, SAvg / 1e6, "k-o", markersize=0.5, linewidth=0.5, label='here')
    plt.plot(
        jiang_data[:, 0] - jiang_data[jiang_t_align, 0],
        jiang_stress,
        "b-*",
        markersize=0.5,
        linewidth=0.5,
        label='Jiang'
    )
    plt.legend()
    plt.xlim(lims)
    plt.xlabel("Time (s)")
    plt.ylabel("Shear Stress (MPa)")
    # plt.show()

    plt.subplot(2, 1, 2)
    plt.plot(t_chunk - t_align, VAvg, "k-o", markersize=0.5, linewidth=0.5, label='here')
    plt.plot(
        jiang_data[:, 0] - jiang_data[jiang_t_align, 0],
        jiang_slip_rate[:],
        "b-*",
        markersize=0.5,
        linewidth=0.5,
        label='Jiang'
    )
    plt.legend()
    plt.xlim(lims)
    plt.xlabel("Time (s)")
    plt.ylabel("Slip rate (m/s)")
    plt.tight_layout()
    plt.show()

Why don’t the solutions match more precisely? I think the most likely answer is that the Jiang result is using a 3D model that is very long in the \(z\) direction but not truly infinite like the antiplane model here.

Integral equation tutorials

Quasidynamic earthquake cycles: antiplane

Contents