Examples#
Below are example input files for motorcycle-3d-ratestate to simulate cycles of slow-slip events with one or two faults, respectively.
Example input file 1#
The following example is for the simulation of recurring slow-slip events along a single two-dimensional fault embedded in a three-dimensional full elastic space. The characteristic nucleation size is h* = 3750 m and the characteristic length of the velocity-weakening region is W = 10 km, leading to a Dieterich-Ruina-Rice number Ru = 2.66. The choice of frictional parameters leads to Rb = 0.285.
#!/bin/bash -e
# slow-slip events along a single fault
# h* = 3750 m; W = 10 km; Ru = 2.66; Rb = 0.285
# to avoid extremely large input files, the code allows to reuse
# the physical properties of previous patches. Normal input for
# physical properties consists of a line number followed by the
# expected physical properties. To reuse the properties of the
# previous patch, use minus the line number. We implement this
# strategy in the awk code below, where we automatically check
# whether the new parameters are different from the previous one.
selfdir=$(dirname $0)
# output directory
WDIR=$selfdir/output1
if [ ! -e "$WDIR" ]; then
echo "adding directory $WDIR"
mkdir -p "$WDIR"
fi
# number of fault patches
N1=256
N2=128
# patch size
DX=400
cat <<EOF > $WDIR/in.param
# output directory
$WDIR
# elastic moduli (lambda, rigidity)
30e3 30e3
# time interval
3.15e10
# number of faults
1
# grid dimension (N1,N2)
$N1 $N2
# sampling (dx1,dx2)
$DX $DX
# n tau0 mu0 sig a b L Vo G/(2Vs) Vl rake
$(echo "" | awk -v n1=$N1 -v n2=$N2 -v dx=$DX '
function abs(x){return (x>0)?x:-x};
function max(x,y){return (x>y)?x:y};
function boxcar(x){return (x>=-0.5 && x<=0.5)?1:0};
function heavi(x){return (x>0)?1:0};
function ramp(x){return x*boxcar(x-0.5)+heavi(x-1)};
function asinh(x){return log(x+sqrt(1+x^2))};
function sinh(x){return (exp(x)-exp(-x))/2};
BEGIN{
c=1;
# initial parameters for reuse
tau0_p=-1;mu0_p =-1;sig_p=-1;
a_p=-1;b_p=-1;L_p=-1;
Vo_p=-1;damping_p=-1;
Vl_p=-1;rake_p=-1;dirichlet_p="T";
}{
for (i2=0;i2<n2;i2++){
for (i1=0;i1<n1;i1++){
x1=(i1-n1/2)*dx;
x2=(i2-n2/2)*dx;
# normal of initial shear traction (MPa); use -1 to use steady-state value
tau0=-1;
# characteristic weakening distance (m)
L=0.05;
# rate-dependent parameter (unitless)
a=1e-2;
# state-dependent parameter (unitless)
b=(abs(x1)<37.5e3 && abs(x2)<5e3)?a+4.0e-3:a-4e-3;
# reference coefficient of friction (unitless)
mu0=0.6;
# effective normal stress (MPa)
sig=1e2;
# reference slip-rate (m/s)
Vo=1e-6;
# loading rate (m/s)
Vl=1e-9;
# rake angle for the loading rate (degree)
rake=0;
# radiation damping coefficient ( G/(2 Vs) in units of MPa/m*s )
damping=5;
# whether to apply Dirichlet boundary condition (constant velocity Vl)
dirichlet=((i2<10) || (i2>=(n2-10)) || (i1<10) || (i1>=(n1-10)))?"T":"F";
if (1==boxcar((x1+20e3)/15e3)*boxcar((x2-0e3)/15e3)){
# initial stress triggers nucleation (MPa)
tau0=mu0*sig*(1.1e-9/Vo)^(a/mu0)*(L/Vl)^(b/mu0);
}
# check whether all parameters are identical to that of the previous patch
if ((tau0_p == tau0) && (mu0_p==mu0) && (sig_p==sig) && (a_p==a) && (b_p==b) && (L_p==L) &&
(Vo_p==Vo) && (damping_p==damping) && (Vl_p==Vl) && (rake_p==rake) && (dirichlet_p==dirichlet)){
# print minus line number to save space (value of previous patch is used in code)
printf "%5d\n",-c;
} else {
# print new set of parameters
printf "%5d %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %d %10.2e %10.2e %s\n",
c, tau0, mu0, sig, a, b, L, Vo, damping, Vl, rake, dirichlet;
}
c++;
# save current parameters for reuse
tau0_p=tau0;mu0_p=mu0;sig_p=sig;
a_p=a;b_p=b;L_p=L;
Vo_p=Vo;damping_p=damping;
Vl_p=Vl;rake_p=rake;dirichlet_p=dirichlet;
}
}
}')
# number of observation patches
1
# n fault i1 i2 rate
1 1 $(echo "" | awk -v n1=$N1 -v n2=$N2 '{print n1/2+1,n2/2+1,1}')
# number of profiles
2
# n fault index direction rate
1 1 $(echo $N2 | awk '{print int($1/2)+1}') 1 20
2 1 $(echo $N1 | awk '{print int($1/2)+1}') 2 20
# number of events
0
EOF
# use --verbose=2 to output all parameters
# use --verbose=1 to output only parameters different from previous patch
mpirun -n 2 motorcycle-3d-ratestate $* \
--verbose 1 --epsilon 1e-6 --friction-law 1 \
--export-netcdf --export-state \
--maximum-step 3.15e7 --maximum-iterations 40000 \
$WDIR/in.param
Example input file 2#
This example is for cycles of slow-slip events along 2 parallel faults separated by 15 km. The characteristic nucleation size is h* = 1500 m and the velocity-weakening region is 5 km wide, leading to a Dieterich-Ruina-Rice number Ru = 3.
#!/bin/bash -e
# slow slip along 2 parallel faults
# h* = 1500 m; W = 5 km; Ru = 3; Rb = 0.285
selfdir=$(dirname $0)
WDIR=$selfdir/output2
if [ ! -e "$WDIR" ]; then
echo "adding directory $WDIR"
mkdir -p "$WDIR"
fi
# number of fault patches
N1=128
N2=128
# patch size
DX=100
cat <<EOF > $WDIR/in.param
# output directory
$WDIR
# elastic moduli (lambda, rigidity)
30e3 30e3
# time interval
3.15e10
# number of faults
2
# grid dimension (N1,N2)
$N1 $N2
# sampling (dx1,dx2)
$DX $DX
# n tau0 mu0 sig a b L Vo G/(2Vs) Vl rake
$(echo "" | awk -v n1=$N1 -v n2=$N2 -v dx=$DX '
function abs(x){return (x>0)?x:-x};
function max(x,y){return (x>y)?x:y};
function boxcar(x){return (x>=-0.5 && x<=0.5)?1:0};
function heavi(x){return (x>0)?1:0};
function ramp(x){return x*boxcar(x-0.5)+heavi(x-1)};
function asinh(x){return log(x+sqrt(1+x^2))};
function sinh(x){return (exp(x)-exp(-x))/2};
BEGIN{
c=1;
}{
for (i2=0;i2<n2;i2++){
for (i1=0;i1<n1;i1++){
x1=(i1-n1/2)*dx;
x2=(i2-n2/2)*dx;
# normal of initial shear traction (MPa); use -1 to use steady-state value
tau0=-1;
# characteristic weakening distance (m)
L=0.02;
# rate-dependent parameter (unitless)
a=1e-2;
# state-dependent parameter (unitless)
b=(abs(x1)<2.5e3 && abs(x2)<2.5e3)?a+4.0e-3:a-4e-3;
# reference coefficient of friction (unitless)
mu0=0.6;
# effective normal stress (MPa)
sig=1e2;
# reference slip-rate (m/s)
Vo=1e-6;
# loading rate (m/s)
Vl=1e-9;
# rake angle for the loading rate (degree)
rake=0;
# radiation damping coefficient ( G/(2 Vs) in units of MPa/m*s )
damping=5;
# whether to apply Dirichlet boundary condition (constant velocity Vl)
dirichlet=((i2<10) || (i2>=(n2-10)) || (i1<10) || (i1>=(n1-10)))?"T":"F";
if (1==boxcar((x1-0.5e3)/1e3)*boxcar((x2-0.5e3)/1e3)){
# initial stress triggers nucleation (MPa)
tau0=mu0*sig*(1.1e-9/Vo)^(a/mu0)*(L/Vl)^(b/mu0);
}
printf "%5d %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %d %10.2e %10.2e %s\n",
c, tau0, mu0, sig, a, b, L, Vo, damping, Vl, rake, dirichlet;
c++;
}
}
}')
# distance from fault 1
15e3
# n tau0 mu0 sig a b L Vo G/(2Vs) Vl rake
$(echo "" | awk -v n1=$N1 -v n2=$N2 -v dx=$DX '
function abs(x){return (x>0)?x:-x};
function max(x,y){return (x>y)?x:y};
function boxcar(x){return (x>=-0.5 && x<=0.5)?1:0};
function heavi(x){return (x>0)?1:0};
function ramp(x){return x*boxcar(x-0.5)+heavi(x-1)};
function asinh(x){return log(x+sqrt(1+x^2))};
function sinh(x){return (exp(x)-exp(-x))/2};
BEGIN{
c=1;
}{
for (i2=0;i2<n2;i2++){
for (i1=0;i1<n1;i1++){
x1=(i1-n1/2)*dx;
x2=(i2-n2/2)*dx;
# normal of initial shear traction (MPa); use -1 to use steady-state value
tau0=-1;
# characteristic weakening distance (m)
L=0.02;
# rate-dependent parameter (unitless)
a=1e-2;
# state-dependent parameter (unitless)
b=(abs(x1)<2.5e3 && abs(x2)<2.5e3)?a+4.0e-3:a-4e-3;
# reference coefficient of friction (unitless)
mu0=0.6;
# effective normal stress (MPa)
sig=1e2;
# reference slip-rate (m/s)
Vo=1e-6;
# loading rate (m/s)
Vl=1e-9;
# rake angle for the loading rate (degree)
rake=-1;
# radiation damping coefficient ( G/(2 Vs) in units of MPa/m*s )
damping=5;
dirichlet=((i2<10) || (i2>=(n2-10)) || (i1<10) || (i1>=(n1-10)))?"T":"F";
printf "%5d %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %d %10.2e %10.2e %s\n",
c, tau0, mu0, sig, a, b, L, Vo, damping, Vl, rake, dirichlet;
c++;
}
}
}')
# number of observation patches
2
# n nFault i1 i2 rate
1 1 $(echo "" | awk -v n1=$N1 -v n2=$N2 '{print n1/2+1,n2/2+1,1}')
2 2 $(echo "" | awk -v n1=$N1 -v n2=$N2 '{print n1/2+1,n2/2+1,1}')
# number of profiles
4
# n fault index direction rate
1 1 $(echo $N2 | awk '{print int($1/2)+1}') 1 20
2 1 $(echo $N1 | awk '{print int($1/2)+1}') 2 20
3 2 $(echo $N2 | awk '{print int($1/2)+1}') 1 20
4 2 $(echo $N1 | awk '{print int($1/2)+1}') 2 20
# number of events
0
EOF
# use --verbose=2 to output all parameters
# use --verbose=1 to output only parameters different from previous patch
mpirun -n 2 motorcycle-3d-ratestate \
--verbose 1 \
--epsilon 1e-4 \
--export-netcdf \
--maximum-step 3.15e7 \
--maximum-iterations 100000 \
--friction-law 1 \
$WDIR/in.param