Problem statement

We would like to investigate the stress interaction between a 2D plane transform fault and a 3D mantle and how it affects the seismic pattern.

using Oetqf, SpecialFunctions, Optim

Generate the meshes

Generate the mesh for the transform fault, which is suitable for using the Okada (1992) equation. The fault is 80 km long, 8 km deep, with grid sizes of 10 km and 2 km respectively, and a dip angle of 90 degrees (vertical).

mf = gen_mesh(Val(:RectOkada), 80e3, 8e3, 10e3, 2e3, 90.0);

The mesh mf is a RectOkadaMesh, which contains the fault geometry, centroid coordinates, and other properties.

Use Gmsh to generate the mantle mesh, which is suitable for using the Barbot et al. (2017) equation. The volume is 80 km long, 5 km wide, and 14 km deep, with a grid size of 4 cells along the x direction, 3 cells along the y direction, and 3 cells along the z direction. There is no refinement in the x or y direction, while cell sizes are 1.5 times progressively larger along the z axis.

gen_gmsh_mesh(Val(:BEMHex8Mesh), -40e3, -2.5e3, -8e3, 80e3, 5e3, -22e3, 4, 3, 3;
    output = joinpath(@__DIR__, "mantle.vtk"),
    rfzh = cumprod(ones(3) * 1.5), rfy = 1.0, rfyType = "Bump"
)
ma = gen_mesh(Val(:BEMHex8Mesh), joinpath(@__DIR__, "mantle.vtk"));

Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 10%] Meshing curve 2 (Line)
Info    : [ 20%] Meshing curve 3 (Line)
Info    : [ 30%] Meshing curve 4 (Line)
Info    : [ 40%] Meshing curve 6 (Extruded)
Info    : [ 50%] Meshing curve 7 (Extruded)
Info    : [ 50%] Meshing curve 8 (Extruded)
Info    : [ 60%] Meshing curve 9 (Extruded)
Info    : [ 70%] Meshing curve 11 (Extruded)
Info    : [ 80%] Meshing curve 12 (Extruded)
Info    : [ 90%] Meshing curve 16 (Extruded)
Info    : [100%] Meshing curve 20 (Extruded)
Info    : Done meshing 1D (Wall 0.000216343s, CPU 0s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Transfinite)
Info    : [ 20%] Meshing surface 13 (Extruded)
Info    : [ 40%] Meshing surface 17 (Extruded)
Info    : [ 50%] Meshing surface 21 (Extruded)
Info    : [ 70%] Meshing surface 25 (Extruded)
Info    : [ 90%] Meshing surface 26 (Extruded)
Info    : Done meshing 2D (Wall 0.00025792s, CPU 0s)
Info    : Meshing 3D...
Info    : Meshing volume 1 (Extruded)
Info    : Done meshing 3D (Wall 0.000121736s, CPU 0s)
Info    : Optimizing mesh...
Info    : Done optimizing mesh (Wall 2.445e-06s, CPU 0s)
Info    : 80 nodes 150 elements
Info    : Writing '/home/runner/work/Oetqf.jl/Oetqf.jl/docs/build/generated/mantle.vtk'...
Info    : Done writing '/home/runner/work/Oetqf.jl/Oetqf.jl/docs/build/generated/mantle.vtk'
Info    : Reading '/home/runner/work/Oetqf.jl/Oetqf.jl/docs/build/generated/mantle.vtk'...
Info    : Reading 80 points
Info    : Reading 150 cells
Info    : Done reading '/home/runner/work/Oetqf.jl/Oetqf.jl/docs/build/generated/mantle.vtk'

The mesh ma is a BEMHex8Mesh, which contains the mantle geometry, centroid coordinates, and other properties.

Compute the stress Green's functions

Assume the shear modulus and the Lamé parameter are both 3e10 Pa.

λ = μ = 3e10;

Initialize the path to save the stress Green's functions.

gffile = joinpath(@__DIR__, "gf.h5")
isfile(gffile) && rm(gffile);

Compute the stress Green's functions within the fault. We add a buffer zone of 1 times the fault length on both sides of the fault to avoid edge effects.

@time gf₁₁ = stress_greens_function(mf, λ, μ; buffer_ratio = 1)
h5write(gffile, "gf₁₁", gf₁₁); # fault -> fault

  1.066781 seconds (2.30 M allocations: 109.192 MiB, 1.53% gc time, 99.89% compilation time)

Compute the stress Green's functions from the fault to the mantle.

@time gf₁₂ = stress_greens_function(mf, ma, λ, μ; buffer_ratio = 1, qtype = "Gauss1")
h5write(gffile, "gf₁₂", gf₁₂); # fault -> mantle

  0.444816 seconds (255.74 k allocations: 12.471 MiB, 2.87% gc time, 87.50% compilation time)

Compute the stress Green's functions from the mantle to the fault and within the mantle.

@time gf₂₁ = stress_greens_function(ma, mf, λ, μ)
h5write(gffile, "gf₂₁", gf₂₁); # mantle -> fault

  4.709455 seconds (10.76 M allocations: 458.429 MiB, 1.33% gc time, 90.95% compilation time)

Compute the stress Green's functions within the mantle, using Gauss1 quadrature.

@time gf₂₂ = stress_greens_function(ma, λ, μ; qtype = "Gauss1")
h5write(gffile, "gf₂₂", gf₂₂); # mantle -> mantle

Maximum real part of eigval is: 1111120100.9504
  3.966109 seconds (1.20 M allocations: 63.237 MiB, 0.25% gc time, 26.61% compilation time)

Set up parameters

Set up the rate-and-state friction parameters on the fault. We include two velocity-weakening patches on the fault: one on the left side between -25 km and -5 km, and the other on the right side between 5 km and 25 km, and the depth between -6 km and -1 km.

cs = 3044.14 # m/s
vpl = 140e-3 / 365 / 86400 # 140 mm/yr
v0 = 1e-6
f0 = 0.6
μ = 3e10
η = μ / 2cs # radiation damping
ν = λ / 2(λ + μ)
avw = 0.015
abvw = 0.0047
Dc = 8e-3
σmax = 5e7
a = ones(mf.nx, mf.nξ) .* avw
b = ones(mf.nx, mf.nξ) .* (avw - abvw)
L = ones(mf.nx, mf.nξ) .* Dc
σ = [min(σmax, 1.5e6 + 18.0e3 * z) for z in -mf.z] # Pa
σ = repeat(σ, 1, mf.nx)' |> Matrix # Pa
left_patch = @. -25.e3 ≤ mf.x ≤ -5.e3
right_patch = @. 5.e3 ≤ mf.x ≤ 25.e3
vert_patch = @. -6.e3 ≤ mf.z ≤ -1e3
b[xor.(left_patch, right_patch), vert_patch] .= avw + abvw # assign velocity weakening
pf = RateStateQuasiDynamicProperty(a, b, L, σ, η, vpl, f0, v0)
save_property(joinpath(@__DIR__, "para-fault.bson"), pf);

Set up rheology parameters in the mantle assuming power-law viscosity with lab-derived results.

A_wet_dis = 3e1
Q_wet_dis = 480e3
V_wet_dis = 11e-6
m_wet_dis = 0
r_wet_dis = 1.2
n_wet_dis = 3.5
grain_size = 10000.0 # μm
COH = 1000 # ppm / HSi
𝙍 = 8.314 # gas constant
crust_depth = 7e3
κ = 8e-7
𝚃(z) = 1673 * erf(z / sqrt(4κ * 1e6 * 365 * 86400)) # 1 Myr OTF
𝙿(z) = 2800 * 9.8 * crust_depth + 3300 * 9.8 * (z - crust_depth)
prefactor_dis(z) = A_wet_dis / (1e6)^n_wet_dis * COH^r_wet_dis * grain_size^m_wet_dis * exp(-(Q_wet_dis + 𝙿(z) * V_wet_dis) / 𝙍 / 𝚃(z))
rel_dϵ = [0.0, -1e-12, 0.0, 0.0, 0.0, 0.0]
amplifier = 1e0
γ_dis = prefactor_dis.(-ma.cz) .* amplifier
pa = PowerLawViscosityProperty(γ_dis, ones(length(ma.cz)) * (n_wet_dis - 1), rel_dϵ) # note to save `n-1` instead of `n`, where `n` refers to the stress power
save_property(joinpath(@__DIR__, "para-mantle" * ".bson"), pa);

Set up initial conditions

We set up a uniform initial velocity field equal to the plate rate on the fault, with slight perturbations in the state variable of the left and right halves of the fault.

vinit = pf.vpl .* ones(size(pf.a))
θinit = pf.L ./ vinit
θinit[1: size(θinit, 1) >> 1, :] ./= 1.1
θinit[size(θinit, 1) >> 1 + 1: end, :] ./= 2.5
δinit = zeros(size(pf.a));

We set up initial conditions in the mantle where the initial stress matches the background strain rate through a simple 1D optimization.

ϵinit = zeros(length(pa.γ), 6)
P = map(z -> 2800 * 9.8 * crust_depth + 3300 * 9.8 * (z - crust_depth), -ma.cz) # change the depth of crust
σinit = repeat(P, 1, 6)
σinit[:,3] .= 0.0 # xz
σinit[:,5] .= 0.0 # yz
target(i) = x -> (pa.γ[i] * (sqrt(2) * x) ^ (pa.n[i]) * x - abs(pa.dϵ₀[2])) ^ 2
σxyinit = -map(i -> Optim.minimizer(optimize(target(i), 1e1, 1e14)), 1: length(pa.γ))
reldϵ = map(i -> pa.γ[i] * (sqrt(2) * abs(σxyinit[i])) ^ (pa.n[i]) * σxyinit[i], 1: length(pa.γ))
@assert all(isapprox.(reldϵ, pa.dϵ₀[2]; rtol=1e-3))
σinit[:,2] .= σxyinit;

Assemble the problem

Using the Green's functions, the properties, and the initial conditions, we assemble the problem. Notice the order of the variables here must be velocity, state variable, strain, stress, and slip.

uinit = ArrayPartition(vinit, θinit, ϵinit, σinit, δinit)
prob = assemble(gf₁₁, gf₁₂, gf₂₁, gf₂₂, pf, pa, uinit, (0.0, 0.1 * 365 * 86400));

We set up the saving scheme and solve the equation. Here, we will save (in the order of) velocity, state variable, strain rate, strain, stress, and slip, every 100 steps. All the variables are named exactly the same as in the equations.

handler(u::ArrayPartition, t, integrator) = (u.x[1], u.x[2], integrator(integrator.t, Val{1}).x[3], u.x[3], u.x[4], u.x[5])
output = joinpath(@__DIR__, "output.h5")
@time sol = wsolve(prob, VCABM5(), output, 100, handler, ["v", "θ", "dϵ", "ϵ", "σ", "δ"], "t";
    reltol=1e-6, abstol=1e-8, dtmax=0.2*365*86400, dt=1e-8, maxiters=1e9, stride=100, force=true
)

retcode: Success
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
 0.0
 3.1536e6
u: 2-element Vector{ArrayPartition{Float64, NTuple{5, Matrix{Float64}}}}:
 ([4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9; 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9; … ; 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9; 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9 4.439370877727043e-9], [1.6382337662337658e6 1.6382337662337658e6 1.6382337662337658e6 1.6382337662337658e6; 1.6382337662337658e6 1.6382337662337658e6 1.6382337662337658e6 1.6382337662337658e6; … ; 720822.857142857 720822.857142857 720822.857142857 720822.857142857; 720822.857142857 720822.857142857 720822.857142857 720822.857142857], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [2.993126315789474e8 -1.159793132258934e6 … 0.0 2.993126315789474e8; 4.865442105263158e8 -279723.0841659712 … 0.0 4.865442105263158e8; … ; 4.865442105263158e8 -279723.0841659712 … 0.0 4.865442105263158e8; 7.673915789473685e8 -235923.5595516861 … 0.0 7.673915789473685e8], [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; … ; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0])
 ([4.076715958989289e-9 4.072616164245145e-9 4.075442712126598e-9 4.088964986453331e-9; 3.4430994939670944e-9 3.417383507269353e-9 3.4394142408770543e-9 4.062644179678046e-9; … ; 7.586803484447171e-10 7.350000896538302e-10 7.495634381963362e-10 2.0502107855143374e-9; 2.0738117461621434e-9 2.038088880757506e-9 2.047039752307764e-9 2.096094491124679e-9], [1.8590418302568565e6 1.8598256460833005e6 1.8593496132564675e6 1.8570189571762218e6; 2.0017790265015583e6 2.0064139281027382e6 2.002447943468171e6 1.8612812874577832e6; … ; 2.8747044541461663e6 2.885326781900475e6 2.8786601284604673e6 2.2772683705092487e6; 2.270762754954703e6 2.279965869236065e6 2.2775966294517145e6 2.2647807208464365e6], [-1.9490353719142703e-10 -3.150225479742601e-6 … 1.117009381251456e-10 -5.7605857040498805e-12; -3.143400770248431e-10 -3.147648023609172e-6 … 4.658202337483857e-10 -4.579661948262328e-11; … ; -5.615414367667751e-10 -3.1305062949498476e-6 … -9.98692990183044e-10 -6.379063666159174e-10; -1.2895020205790048e-10 -3.1447258792037835e-6 … -6.507466105621481e-10 -2.713827511764087e-10], [2.993123106670286e8 -1.1588113187909983e6 … 109.26822681446149 2.993125101945185e8; 4.865440819229592e8 -279332.7091081294 … 107.15415541021898 4.865441499927149e8; … ; 4.8654392870839775e8 -278353.03260683204 … -211.77388305242889 4.8654394576009667e8; 7.673915013824518e8 -235463.28389308814 … -117.59525704444525 7.673914887318344e8], [0.013352977245283152 0.01334786179521655 0.01335086397226022 0.013365680395724845; 0.012402031284994014 0.012375690899781236 0.012398238171474936 0.013339120423498973; … ; 0.005481590194932253 0.00543582127551878 0.005464691937056686 0.008821136502916271; 0.008856835815734015 0.008807659239440364 0.008820414423926334 0.008889500625309089])

Analyze the results

The solution is saved in the output.h5 file, which contains the time series of velocity, state variable, strain rate, strain, stress, and slip. We can load the solution and analyze the results, for example, extracting the earthquake catalog from the velocity time series, visualizing the fault rupture, mantle strain flow, etc. Readers are encouraged to explore the figures in the Shi et al., 2022 for more insights.

gen_pvd(mf, joinpath(@__DIR__, "mantle.vtk"), output, "t", ["v", "θ"], ["dϵ"], 1: length(sol.t), joinpath(@__DIR__, "sol.pvd"))

This page was generated using Literate.jl.