---
Date: Aug 31, 2023
Title: motorcycle-3d-ratestate documentation
---

## 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
~~~