Quasidynamic earthquake cycles on two faults
Contents
7. Quasidynamic earthquake cycles on two faults¶
7.1. Summary¶
Much like the last section, we simulate quasidynamic rupture evolution but this time on two faults, one of which is curved.
Most of the code here follows almost exactly from the previous section on strike-slip/antiplane earthquake cycles.
Realistic tectonic loading is a challenge when there is more than one fault.
The boundary integral approach used here is barely more complex with two faults than with one fault.
7.2. Loading a multi-fault system¶
In the single fault system from the last section, we were able to gloss over loading by using a slip deficit or backslip formulation. Because the physical system is invariant with respect to a rigid body displacement, any tectonic motion that has already been fully accomodated by fault slip can be ignored. Put differently, in a single fault system with a purely elastic rheology, the stresses from tectonic motion can be equivalently described by backwards slip on the fault.
(7.1)¶\[\begin{equation}
u_{\textrm{total}} = u_{\textrm{loading}} + u_{\textrm{fault}}
\end{equation}\]
There are a number of different options for loading a fault system:
backslip: as discussed, we simply apply a backwards slip rate to account for distant plate motion. Once the stresses from this backwards slip accumulate, the fault slips forward.
farfield displacement: apply a constant velocity boundary condition at a distance from the fault.
stress: Either on a lower boundary of the domain or on distant side boundaries, apply a constant traction/stress loading rate.
block model: this is a more advanced version of backslip loading where each fault-bounded subregion of the domain has its own block rotation/motion rate. These rates are assumed to be the long term slip rate. A corresponding kinematic slip deficit is applied to each fault.
Here, I do something a bit like the block model solution by just applying a backslip of 31 mm of loading per year (1 nm of loading per second) to both faults equally. This would model a real-world situation with a plate boundary that has 62 mm per year of strike-slip motion. Exactly 50% of that motion is accomodated by each of the two faults.
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from tectosaur2 import (
gauss_rule,
integrate_term,
panelize_symbolic_surface,
concat_meshes,
)
from tectosaur2.laplace2d import double_layer, hypersingular
shear_modulus = 3.2e10
quad_rule = gauss_rule(5)
t = sp.var("t")
fault_bottom = 40000
fault1 = panelize_symbolic_surface(
t, -5000 + t * 0, fault_bottom * (t + 1) * -0.5,
quad_rule,
n_panels=200
)
fault2_x = 5000 + (t * 10000 + 10000) + (t ** 2 * 5000 - 5000)
fault2 = panelize_symbolic_surface(
t, fault2_x, fault_bottom * (t + 1) * -0.5,
quad_rule,
n_panels=200
)
faults = concat_meshes((fault1, fault2))
surf_half_L = 100000
free = panelize_symbolic_surface(
t, -t * surf_half_L, 0 * t,
quad_rule,
n_panels=200
)
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 {faults.n_panels} panels with a total of {faults.n_pts} points."
)
The free surface mesh has 200 panels with a total of 1000 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(fault1.pts[:,0]/1000, fault1.pts[:,1]/1000, 'r-')
plt.plot(fault2.pts[:,0]/1000, fault2.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()
singularities = np.array(
[
[-surf_half_L, 0],
[surf_half_L, 0],
[5000, 0],
[-5000, 0]
]
)
(free_disp_to_free_disp, fault_slip_to_free_disp), report = integrate_term(
double_layer, free.pts, free, faults, tol=1e-11, 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]
)
(free_disp_to_fault_stress, fault_slip_to_fault_stress), report = integrate_term(
hypersingular,
faults.pts,
free,
faults,
tol=1e-10,
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] * faults.normals[:, :, None], axis=1
)
free_disp_to_fault_traction = np.sum(
free_disp_to_fault_stress[:,:,:,0] * faults.normals[:, :, None], axis=1
)
/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))
analytical_sd = -np.arctan(-fault_bottom / free.pts[:,0]) / np.pi
numerical_sd = Dinv.dot(-C).sum(axis=1)
rp = (faults.pts[:,1] + 1) ** 2
analytical_ft = -(shear_modulus / (2 * np.pi)) * ((1.0 / (faults.pts[:,1] + fault_bottom)) - (1.0 / (faults.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')
plt.plot(faults.pts[:,1], np.log10(np.abs(numerical_ft)), 'b-')
plt.xlim([-fault_bottom, 0])
plt.xlabel('$\mathrm{depth ~~ (m)}$')
plt.ylabel('$t$')
plt.subplot(1,2,2)
plt.title('Displacement')
plt.plot(free.pts[:,0], numerical_sd, 'b-')
plt.xlim([-100000, 100000])
plt.xlabel('$x ~~ \mathrm{(m)}$')
plt.ylabel('$u$')
plt.show()
7.3. Rate and state friction¶
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 = faults.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()
# 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
for fm in [fault1, fault2]:
mesh_L = np.max(np.abs(np.diff(fm.pts[:, 1])))
Lb = shear_modulus * fp_Dc / (sigma_n * fp_b)
hstar = (np.pi * shear_modulus * fp_Dc) / (sigma_n * (fp_b - fp_a))
print(mesh_L, Lb, np.min(hstar[hstar > 0]))
53.84693101057201 341.3333333333333 3216.990877275949
53.84693101057201 341.3333333333333 3216.990877275949
7.4. Quasidynamic earthquake cycle derivatives¶
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(faults.n_pts, init_state_scalar)
tau_amax = -qd_equation(amax, 0, Vp, init_state_scalar)
init_traction = np.full(faults.n_pts, tau_amax)
init_slip_deficit = np.zeros(faults.n_pts)
init_conditions = np.concatenate((init_slip_deficit, init_state))
class SystemState:
V_old = np.full(faults.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(SS):
"""This is just a helper function that creates some rough plots of the
current state to help with debugging"""
slip_deficit, state, tau_qs, slip_deficit_rate, dstatedt = SS
plt.figure(figsize=(15, 9))
plt.subplot(2, 3, 1)
plt.title("slip")
plt.plot(faults.pts[:, 1], slip_deficit)
plt.subplot(2, 3, 2)
plt.title("state")
plt.plot(faults.pts[:, 1], state)
plt.subplot(2, 3, 3)
plt.title("tau qs")
plt.plot(faults.pts[:, 1], tau_qs)
plt.subplot(2, 3, 4)
plt.title("slip rate")
plt.plot(faults.pts[:, 1], slip_deficit_rate)
plt.subplot(2, 3, 6)
plt.title("dstatedt")
plt.plot(faults.pts[:, 1], dstatedt)
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
7.5. 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.002737907006988508 yrs, step=0.002737907006988508
step=5000, time=161.2388012289919 yrs, step=6.068323795478569e-11
step=10000, time=161.238801483585 yrs, step=3.200375945922213e-11
step=15000, time=177.06735503165402 yrs, step=0.05782401363757997
step=20000, time=187.60702659014385 yrs, step=8.29256996752649e-11
step=25000, time=187.6070269264095 yrs, step=2.786049654906221e-10
step=30000, time=213.4153146695212 yrs, step=1.3348499105890855e-10
step=35000, time=229.49165911338014 yrs, step=1.8458825569115088e-10
step=40000, time=298.5155537603046 yrs, step=1.1030568652139828e-10
step=45000, time=298.51555399934733 yrs, step=3.251751197178755e-11
step=50000, time=298.5168884817506 yrs, step=6.51359998844683e-06
step=55000, time=323.01280183052285 yrs, step=1.2281707123916784e-10
step=60000, time=323.01280207887106 yrs, step=3.3484575524851865e-11
step=65000, time=323.2817940109216 yrs, step=0.001729718400292269
step=70000, time=351.52626665410605 yrs, step=1.0785478482910092e-08
step=75000, time=368.9516542294202 yrs, step=1.2976784052681758e-10
step=80000, time=368.951654876993 yrs, step=6.515590688770813e-11
step=85000, time=428.9915445603336 yrs, step=0.0014206634043130257
step=90000, time=446.2939275656298 yrs, step=6.588120455250637e-11
step=95000, time=446.29392801977554 yrs, step=9.456068304806992e-10
step=100000, time=469.8994523682489 yrs, step=1.0256917809688377e-10
step=105000, time=469.89945268723017 yrs, step=5.1798341561007304e-11
step=110000, time=469.9812354270467 yrs, step=0.00034767164051762437
step=115000, time=500.90614938849876 yrs, step=1.439715864624497e-10
step=120000, time=573.9064649517662 yrs, step=1.424001081887202e-10
step=125000, time=573.9064651906277 yrs, step=2.5022769435539118e-11
step=130000, time=597.5674917048766 yrs, step=0.0012973735857483035
step=135000, time=597.6230358626846 yrs, step=5.391379308333549e-11
step=140000, time=597.6230366893858 yrs, step=4.587749495737104e-09
step=145000, time=625.2112327038584 yrs, step=1.1906969997104363e-10
step=150000, time=642.3455004578122 yrs, step=1.2898210138995284e-10
step=155000, time=704.8417989867936 yrs, step=0.003113930973116129
step=160000, time=716.3612963754118 yrs, step=2.9978970144993726e-11
step=165000, time=716.3612966987445 yrs, step=3.9843018386249724e-10
step=170000, time=743.0069768079197 yrs, step=1.6222491102653862e-10
step=175000, time=743.0069772376432 yrs, step=4.6660816435353134e-11
step=180000, time=743.0070823319467 yrs, step=6.836717438689674e-07
step=185000, time=772.4229461837367 yrs, step=1.1205848921132735e-10
step=190000, time=790.5550252030326 yrs, step=9.598105764163313e-11
step=195000, time=790.5550257931098 yrs, step=9.694812119469745e-11
step=200000, time=790.5770082962302 yrs, step=0.00011880289499414454
step=205000, time=843.4093589520743 yrs, step=1.3986156636192637e-10
step=210000, time=843.4093595399053 yrs, step=5.7540281407326665e-11
step=215000, time=865.8601485323139 yrs, step=4.5113514750450234e-10
step=220000, time=903.1263077569464 yrs, step=1.2438854951289735e-10
step=225000, time=903.1263082887026 yrs, step=8.413452911659529e-11
step=230000, time=927.4523165504868 yrs, step=0.008045238903279458
step=235000, time=947.0642174585033 yrs, step=1.7528026899290684e-11
step=240000, time=947.0685830652857 yrs, step=1.768826485736451e-05
step=245000, time=972.5636479272499 yrs, step=4.078590535048743e-10
step=250000, time=988.09794803525 yrs, step=8.195863612220059e-11
step=255000, time=988.0979484325295 yrs, step=2.7077779485800782e-11
CPU times: user 2h 15min 3s, sys: 36min 52s, total: 2h 51min 55s
Wall time: 1h 48min 42s
7.6. Plotting the results¶
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()
for j in range(2):
fm = fault1 if j == 0 else fault2
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_full = y[: init_slip_deficit.shape[0]]
fault1_sd = slip_deficit_full[:fault1.n_pts]
fault2_sd = slip_deficit_full[fault1.n_pts:]
slip_deficit = fault1_sd if j == 0 else fault2_sd
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, fm.pts[:,1] / 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()
