Skip to content

Example 1D

Note

This example is from Benchmark Problem 1 (hence referred as BP1).

Define Parameters

First, we load the package

1
2
using JuEQ
using Plots

Instead of using SI unit, we refactor ours into the follow:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
ms2mmyr = 365 * 86400 * 1e3 # convert velocity from m/s to mm/yr
ρ = 2670.0 # density [kg/m³]
vs = 3464.0 # shear wave velocity [m/s]
σ = 500.0 # effective normal stress [bar]
a0 = 0.010 # frictional paramter `a` in vw zone
amax = 0.025 # frictional paramter `a` in vs zone
b0 = 0.015 # frictional paramter `b`
L = 8.0 # critical distance [mm]
vpl = 1e-9 * ms2mmyr # plate rate [mm/yr]
vinit = 1e-9 * ms2mmyr # initial velocity [mm/yr]
v0 = 1e-6 * ms2mmyr # reference velocity [mm/yr]
f0 = 0.6 # reference frictional coefficient
H = 15.0 # vw zone [km]
h = 3.0 # vw-vs changing zone [km]
Wf = 40.0 # fault depth [km]
Δz = 100.0e-3 # grid size interval [km]
tf = 400.0; # simulation time [yr]

1
400.0

Warning

Make sure your units are consistent across the whole variable space. Pontenial imporvement may incoporate Unitful.jl package.

Then we arrive at some parameters that are implicit by above:

1
2
3
4
μ = vs^2 * ρ / 1e5 / 1e6 # shear modulus [bar·km/mm]
λ = μ # poisson material
η = μ / 2(vs * 1e-3 * 365 * 86400)
ngrid = round(Int, Wf / Δz); # number of grids

1
400

Now, we start to construct our model using parameters above. First, we create a 'fault' by specifying fault type and depth:

Tip

Here, we do not need to provide dip for strike-slip fault as it automatically choose 90. See fault.

Construct Model

1
fa = fault(StrikeSlipFault, Wf);

1
PlaneFaultDomain{StrikeSlipFault,1,Float64}(90.0, (40.0,))

Next, we generate the grid regarding the fault we just created by giving number of grids:

Note

This package use ξ for denoting downdip coordinate and x for along-strike one. See discretize.

1
gd = discretize(fa; nξ=ngrid);

1
JuEQ.BEMGrid_1D{Array{Float64,1},Array{Array{Float64,1},1},Float64,Int64}([-0.05, -0.15, -0.25, -0.35, -0.45, -0.55, -0.65, -0.75, -0.85, -0.95  …  -39.05, -39.15, -39.25, -39.35, -39.45, -39.55, -39.65, -39.75, -39.85, -39.95], 0.1, 400, Array{Float64,1}[[-0.1, 0.0], [-0.2, -0.1], [-0.3, -0.2], [-0.4, -0.3], [-0.5, -0.4], [-0.6, -0.5], [-0.7, -0.6], [-0.8, -0.7], [-0.9, -0.8], [-1.0, -0.9]  …  [-39.1, -39.0], [-39.2, -39.1], [-39.3, -39.2], [-39.4, -39.3], [-39.5, -39.4], [-39.6, -39.5], [-39.7, -39.6], [-39.8, -39.7], [-39.9, -39.8], [-40.0, -39.9]], [-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, -0.0, -0.0, -0.0, -0.0], [-0.05, -0.15, -0.25, -0.35, -0.45, -0.55, -0.65, -0.75, -0.85, -0.95  …  -39.05, -39.15, -39.25, -39.35, -39.45, -39.55, -39.65, -39.75, -39.85, -39.95], [-500.0, 500.0])

Next, we construct the required frictional parameter profile:

1
2
3
4
z = -gd.ξ
az = fill(a0, size(z))
az[z .≥ (H + h)] .= amax
az[H .< z .< H + h] = a0 .+ (amax - a0) / (h / Δz) * collect(1: Int(h / Δz));

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
30-element Array{Float64,1}:
 0.0105
 0.011
 0.0115
 0.012
 0.0125
 0.013000000000000001
 0.0135
 0.014
 0.0145
 0.015
 ⋮
 0.020999999999999998
 0.0215
 0.022
 0.0225
 0.023
 0.0235
 0.024
 0.0245
 0.025

Then, we provide the required initial condition satisfying uniform slip distribution over the depth:

1
2
3
4
5
τ0 = σ * amax * asinh(vinit / 2v0 * exp((f0 + b0 * log(v0 / vinit)) / amax)) + η * vinit
τz = fill(τ0, size(z))
θz = @. L / v0 * exp(az / b0 * log(2v0 / vinit * sinh((τz - η * vinit) / az / σ)) - f0 / b0)
vz = fill(vinit, size(z))
u0 = hcat(vz, θz);

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
400×2 Array{Float64,2}:
 31.536  0.000253678
 31.536  0.000253678
 31.536  0.000253678
 31.536  0.000253678
 31.536  0.000253678
 31.536  0.000253678
 31.536  0.000253678
 31.536  0.000253678
 31.536  0.000253678
 31.536  0.000253678
  ⋮
 31.536  0.253678
 31.536  0.253678
 31.536  0.253678
 31.536  0.253678
 31.536  0.253678
 31.536  0.253678
 31.536  0.253678
 31.536  0.253678
 31.536  0.253678

Let's simulate only the first 200 years:

1
tspan = (0., 200.);

1
(0.0, 200.0)

Finally, we provide the material properties w.r.t. our 'fault', 'grid' as well as other necessary parameters predefined using the same grid size & dimension:

1
mp = properties(;fault=fa, grid=gd, parameters=[:a=>az, :b=>b0, :L=>L, :σ=>σ, :η=>η, :k=>[:λ=>λ, :μ=>μ], :vpl=>vpl, :f0=>f0, :v0=>v0]);

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
[ Info: Calculating stiffness tensor ...
[ Info: Fault material properties establised.
JuEQ.PlaneMaterialProperties{1,Float64,Array{Float64,1},Array{Float64,2}}
  dims: Tuple{Int64}
  a: Array{Float64}((400,)) [0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01  …  0.025, 0.025, 0.025, 0.025, 0.025, 0.025, 0.025, 0.025, 0.025, 0.025]
  b: Array{Float64}((400,)) [0.015, 0.015, 0.015, 0.015, 0.015, 0.015, 0.015, 0.015, 0.015, 0.015  …  0.015, 0.015, 0.015, 0.015, 0.015, 0.015, 0.015, 0.015, 0.015, 0.015]
  L: Array{Float64}((400,)) [8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0  …  8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0]
  k: Array{Float64}((400, 400)) [1.35974 -0.815844 … -6.39586e-6 -6.36374e-6; -0.815844 1.98134 … -6.39609e-6 -6.36397e-6; … ; -6.39586e-6 -6.39609e-6 … 2.03961 -0.679871; -6.36374e-6 -6.36397e-6 … -0.679871 2.03961]
  σ: Array{Float64}((400,)) [500.0, 500.0, 500.0, 500.0, 500.0, 500.0, 500.0, 500.0, 500.0, 500.0  …  500.0, 500.0, 500.0, 500.0, 500.0, 500.0, 500.0, 500.0, 500.0, 500.0]
  η: Array{Float64}((400,)) [1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9  …  1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9, 1.4664e-9]
  vpl: Float64 31.536
  f0: Float64 0.6
  v0: Float64 31536.0

Tip

Check properties for extended options.

Check our profile now:

1
plot([mp.a, mp.b], z, label=["a", "b"], yflip=true, ylabel="Depth (km)")

Documenter.Documents.RawHTML("<?xml version=\"1.0\" encoding=\"utf-8\"?>\n\n\n \n \n \n\n\n \n \n \n\n\n\n \n \n \n\n\n\n \n \n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n0.012\n\n\n0.015\n\n\n0.018\n\n\n0.021\n\n\n0.024\n\n\n0\n\n\n10\n\n\n20\n\n\n30\n\n\n40\n\n\nDepth (km)\n\n\n\n\n\n\n\na\n\n\n\nb\n\n\n")

We then contruct the ODEProblem as following by stating which state evolution law to use and frcitonal law form, plus initial condition and simulation time:

1
prob = EarthquakeCycleProblem(gd, mp, u0, tspan; se=DieterichStateLaw(), fform=RForm());

1
2
3
ODEProblem with uType Array{Float64,2} and tType Float64. In-place: true
timespan: (0.0, 200.0)
u0: [31.536 0.000253678; 31.536 0.000253678; … ; 31.536 0.253678; 31.536 0.253678]

Solve Model

We then solve the ODEs:

1
sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6);

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 4814-element Array{Float64,1}:
   0.0
   6.483553518719132e-5
   9.44532839048002e-5
   0.00014950199694886354
   0.0001969017343314171
   0.0002589370267499311
   0.00032362086183610665
   0.0004000242025107322
   0.0004844340682573306
   0.0005812591543481377
   ⋮
 198.66317580661035
 198.81599601090952
 198.9735692145639
 199.13681420817997
 199.30704525039675
 199.48496682770454
 199.66952071102318
 199.8589025257933
 200.0
u: 4814-element Array{Array{Float64,2},1}:
 [31.536 0.000253678; 31.536 0.000253678; … ; 31.536 0.253678; 31.536 0.253678]
 [22.4215 0.000318453; 22.4215 0.000318453; … ; 31.536 0.253678; 31.536 0.253678]
 [19.6236 0.000348045; 19.6236 0.000348045; … ; 31.536 0.253678; 31.536 0.253678]
 [15.747 0.000403048; 15.747 0.000403048; … ; 31.536 0.253678; 31.536 0.253678]
 [13.3297 0.000450411; 13.3297 0.000450411; … ; 31.536 0.253678; 31.536 0.253678]
 [10.9855 0.000512401; 10.9855 0.000512401; … ; 31.536 0.253678; 31.536 0.253678]
 [9.1923 0.000577041; 9.1923 0.000577041; … ; 31.536 0.253678; 31.536 0.253678]
 [7.62907 0.000653395; 7.62907 0.000653395; … ; 31.536 0.253678; 31.536 0.253678]
 [6.35871 0.000737754; 6.35871 0.000737754; … ; 31.536 0.253678; 31.536 0.253678]
 [5.28542 0.000834524; 5.28542 0.000834524; … ; 31.536 0.253678; 31.536 0.253678]
 ⋮
 [0.000251795 3.02087; 0.00025182 3.02087; … ; 39.5126 0.201016; 36.8001 0.216185]
 [0.000236784 3.17367; 0.000236808 3.17367; … ; 39.2862 0.202129; 36.6489 0.217045]
 [0.00022298 3.33123; 0.000223003 3.33123; … ; 39.0508 0.203309; 36.4921 0.217951]
 [0.000210202 3.49446; 0.000210224 3.49446; … ; 38.8067 0.204558; 36.3297 0.218908]
 [0.000198281 3.66468; 0.000198302 3.66468; … ; 38.5537 0.205881; 36.1616 0.219915]
 [0.000187127 3.84258; 0.000187147 3.84258; … ; 38.2928 0.207276; 35.9886 0.220973]
 [0.000176761 4.02712; 0.000176781 4.02712; … ; 38.0279 0.208724; 35.813 0.222065]
 [0.000167214 4.21648; 0.000167232 4.21648; … ; 37.7635 0.210201; 35.638 0.223175]
 [0.000160724 4.35757; 0.000160741 4.35757; … ; 37.5721 0.21129; 35.5113 0.223989]

Tip

For details of solving options, see here.

Tip

Raise the accuracy option if you get instability when solving these ODEs.

Results

The first event happens at around 196 year:

1
2
maxv = max_velocity(sol)
plot(sol.t, log10.(maxv / ms2mmyr), xlabel="Time (year)", ylabel="Max Velocity (log10 (m/s))", xlims=(190, 200), label="")

Documenter.Documents.RawHTML("<?xml version=\"1.0\" encoding=\"utf-8\"?>\n\n\n \n \n \n\n\n \n \n \n\n\n\n \n \n \n\n\n\n \n \n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n190.0\n\n\n192.5\n\n\n195.0\n\n\n197.5\n\n\n200.0\n\n\n-8\n\n\n-6\n\n\n-4\n\n\n-2\n\n\n0\n\n\nTime (year)\n\n\nMax Velocity (log10 (m/s))\n\n\n\n")

Note

Click here for the slip evolution over 3000 years simulation. It may need some time to load the page.

This page was generated using Literate.jl.