# Physics-Informed Deep Learning of Rate-and-State Fault Friction

Cody Rucker<EMAIL_ADDRESS>Brittany A. Erickson

###### Abstract

Direct observations of earthquake nucleation and propagation are few and yet
the next decade will likely see an unprecedented increase in indirect, surface
observations that must be integrated into modeling efforts. Machine learning
(ML) excels in the presence of large data and is an actively growing field in
seismology. However, not all ML methods incorporate rigorous physics, and
purely data-driven models can predict physically unrealistic outcomes due to
observational bias or extrapolation. Our work focuses on the recently emergent
Physics-Informed Neural Network (PINN), which seamlessly integrates data while
ensuring that model outcomes satisfy rigorous physical constraints. In this
work we develop a multi-network PINN for both the forward problem as well as
for direct inversion of nonlinear fault friction parameters, constrained by
the physics of motion in the solid Earth, which have direct implications for
assessing seismic hazard. We present the computational PINN framework for
strike-slip faults in 1D and 2D subject to rate-and-state friction. Initial
and boundary conditions define the data on which the PINN is trained. While
the PINN is capable of approximating the solution to the governing equations
to low-errors, our primary interest lies in the network’s capacity to infer
friction parameters during the training loop. We find that the network for the
parameter inversion at the fault performs much better than the network for
material displacements to which it is coupled. Additional training iterations
and model tuning resolves this discrepancy, enabling a robust surrogate model
for solving both forward and inverse problems relevant to seismic faulting.

###### keywords:

physics-informed neural network , rate-and-state friction , earthquake ,
inverse problem , fully dynamic

[1]organization=Department of Computer Science, University of Oregon,
addressline=1477 E 13th Ave, postcode=97403, city=Eugene, country=United
States of America

[2]organization=Department of Earth Sciences, University of Oregon,
addressline=1272 University of Oregon, city=Eugene, postcode=97403,
country=United States of America

##  1 Context and Motivation

Faults are home to a vast spectrum of event types, from slow aseismic creep,
to slow-slip to megathrust earthquakes followed by postseismic afterslip. The
Cascadia subduction zone in the Pacific Northwest, for example, hosts several
types of slow earthquake processes including low (and very low) frequency
earthquakes, non-volcanic tremor (NVT) and slow-slip events (SSE) [22], but
also large, fast earthquakes, the last of which was a magnitude $\sim$9 in the
year 1700 [1]. Understanding the physical mechanisms for such diversity of
slip styles is crucial for mitigating the associated hazards but major
uncertainties remain in the depth-dependency of frictional properties at fault
zones, which affect fault locking and therefore rupture potential [3, 41].
Direct observations of earthquake nucleation and propagation are few and yet
the next decade will likely see an unprecedented increase in indirect, surface
observations that could be integrated into modeling efforts [2].

Traditional numerical approaches for solving the partial differential
equations (PDE) governing earthquake processes (e.g. finite difference
methods) have seen incredible growth in the past century, in particular in
terms of convergence theory and high-performance computing. Traditional
methods employ a mesh (either a finite number of grid points/nodes or
elements) and a range of time-integration schemes in order to obtain an
approximate solution whose accuracy depends directly on the mesh size (with
error decreasing with decreasing node spacing or element size). This mesh
dependency introduces limitations when high resolution is needed, and while
traditional methods give rise to a forward problem, solving inverse problems
require additional machinery and can be prohibitively expensive [27, 16]. In
addition, noisy or sparse data cannot be seamlessly integrated into the
computational framework of traditional methods.

Machine learning (ML), on the other hand, excels in the presence of large data
and is an actively growing field in seismology, with applications ranging from
earthquake early warning (EEW) to ground‐motion prediction [31]. However, not
all ML methods incorporate rigorous physics, and purely data-driven models can
predict physically unrealistic outcomes due to observational bias or
extrapolation [53]. A new Deep Learning technique has recently emerged called
the Physics-Informed Neural Network (PINN), which seamlessly integrates sparse
and/or noisy data while ensuring that model outcomes satisfy rigorous physical
constraints. PINNs do not outperform traditional numerical methods for forward
problems (except in high-dimensional settings) [26], but they offer advantages
over traditional numerical methods in that both forward and inverse problems
can be solved in the same computational framework. However, the majority of
PINN applications are currently limited to simple mechanical models, forward
problems and/or do not incorporate real-world observations [15, 42, 25]. Here
we introduce a new, physically-rigorous modeling framework for both forward
and inverse problems that can be integrated with observational data in order
to better understand earthquake fault processes. This Deep Learning approach
will vastly expand our computational abilities to explore and infer relevant
parameter spaces responsible for slip complexity, and the conditions that
enable the world’s largest earthquakes.

Though the PINN framework lacks the robust error analysis that comes with
traditional methods, a large number of publications have emerged since
$\sim$2017 which aim to customize PINNs through the use of different
activation functions, gradient optimization techniques, neural network
architecture, and loss function structure [7]. Careful formulation of the loss
function using the weak form of the PDE have been proposed for constructing
deep learning analogues of Ritz [52] and Galerkin [28, 29, 23] methods which
use numerical quadrature to reduce the order of the PDE resulting in a simpler
learning problem [6, 14]. In tandem, statistical learning theory has been used
to deduce global error bounds for PINNs in terms of optimization error,
generalization error, and approximation error [32]. For wide but shallow
networks utilizing hyperbolic tangent activation functions, the approximation
error has been shown to be bounded over Sobolev spaces [9]. Bounds on PINN
generalization error have been derived for linear second-order PDE [47] (later
extended to all linear problems [48]) and some specific cases like Navier-
Stokes [8]. Moreover, the abstract framework for PINNs can leverage stability
of a PDE to provide conditions under which generalization error is small
whenever training error is small for both forward and inverse problems [39,
38]. More recently, a PINN-specific optimization algorithm has achieved
markedly improved accuracy over other optimization algorithms by incorporating
a PDE energy into the backpropagation step [40]. In addition to this rapid
framework development, PINNs have been shown to perform well on a variety of
physical problems like Navier-Stokes [51, 49, 24], convection heat transfer
[5], solid mechanics, [19, 18] and the Euler equations [35].

In this work we focus specifically on rate-and-state friction [e.g. 46], an
experimentally-motivated, nonlinear friction law capable of reproducing a wide
range of observed earthquake behaviors and is used in nearly all modern
dynamic rupture and earthquake cycle simulations [20, 12]. A better
understanding of the depth-dependency of rate-and-state parameters - which
have a direct correlation to fault locking and seismic rupture potential – is
a fundamental task [3, 41]. To address this task we develop a multi-network
PINN for modeling a vertical, strike-slip frictional fault embedded in an
elastic half-space, and consider deformation in both 1D and 2D. The paper is
organized as follows: In section 2 we first provide an overview of the
physics-informed deep learning framework and PINN architecture for general
initial-boundary-value problems, in order to best describe the implementation
to our application problem. In section 3 we provide specific details of the
PINN framework applied problems, first illustrated in 1D with an example
forward problem, then further developed to include inverse problems in 2D. In
section 4 we report details of our optimal network architecture and training
methods, verifying our methods with a manufactured solution to ensure accuracy
of our inversions. We conclude with a summary and discussion of future work in
section 5.

##  2 Physics-Informed Deep Learning Framework

The physics governing motion in many applications in science and engineering
give rise to partial differential equations (PDE) where analytic solutions can
be difficult to obtain (due, e.g. to complex material properties, boundary
conditions, geometry) and we commonly turn to numerical methods. The Physics-
Informed Neural Network (PINN) is a deep learning (DL) framework for
approximating solutions to PDEs. Though this DL framework lacks the robust
mathematical theory of the traditional methods, it shows particular promise in
solving problems for which traditional numerical methods are ill-suited[43].
The DL framework produces a closed, analytic form for the solution, which is
continuous and defined at every point in the domain allowing one to evaluate
the solution “off-grid” without having to resolve the PDE (i.e. it is mesh-
free) [43]. In addition, both forward and inverse problems can be solved in
the same computational setting, as will be shown in subsequent sections.

###  2.1 Feed-forward Deep Neural Networks

PINNs are extensions of a general feed-forward neural network. We let
$\mathbf{x}\in\mathbb{R}^{n}$ and define a weighting matrix
$W\in\mathbb{R}^{m\times n}$ and bias vector $b\in\mathbb{R}^{m}$. A single
hidden layer of a neural network can be expressed as

$\displaystyle\ell(\mathbf{x};\theta)=\varphi(W\mathbf{x}+b),\quad\text{where
}\theta=(W,b),$ (1)

and $\varphi$ is a known (nonlinear) activation function. Deep neural networks
are obtained by repeated composition of hidden layers [30].

We let positive integer $L$ be the deep neural network depth (i.e. number of
hidden layers) and let $\\{\varphi_{i}\\}_{i=1}^{L}$ be a collection of
activation functions along with a sequence of trainable network parameters
$\\{\theta_{i}\\}_{i=0}^{L}$ where $\theta_{k}=(W_{k},b_{k})$ for each $0\leq
k\leq L$. Here we assume each layer consists of $m$ neurons but this
simplifying assumption may be omitted so long as one ensures that each weight
and bias are of the correct dimension for matrix multiplication and addition,
respectively. If $\mathbf{x}\in\mathbb{R}^{n}$ is the network input and
network parameters $(W_{0},b_{0})\in\mathbb{R}^{m\times
n}\times\mathbb{R}^{m}$, $(W_{k},b_{k})\in\mathbb{R}^{m\times
m}\times\mathbb{R}^{m}$ for $1\leq k<L$, and
$(W_{L},b_{L})\in\mathbb{R}^{d\times m}\times\mathbb{R}^{d}$, then the
recursive composition

$\displaystyle\ell_{0}$ $\displaystyle=\mathbf{x},$ (2a)
$\displaystyle\ell_{k}$
$\displaystyle=\varphi_{k}(W_{k}\ell_{k-1}+b_{k})\quad\text{ for }0<k<L,$ (2b)

defines a feed-forward, deep neural network
$\mathcal{N}:\mathbb{R}^{n}\to\mathbb{R}^{d}$ (of depth $L$ and width $m$)
with network parameters $\theta$ by
$\mathcal{N}(\mathbf{x};\theta)=W_{L}\ell_{L-1}+b_{L}$.

###  2.2 PINN Architecture

The PINN architecture is defined by extending the feed-forward deep neural
network to enforce physical conditions set by an initial-boundary value
problem (IBVP). For simplicity in the notation (as well as generality) we
define this extension for a general IBVP, with specific details for our
application given in the subsequent section. Consider the general operator
form of an IBVP given by

$\displaystyle\mathcal{L}\left[{\bf u};\lambda\right]$ $\displaystyle={\bf
k},\quad$ $\displaystyle\text{ in }\widehat{\Omega},$ (3a)
$\displaystyle\mathcal{B}\left[{\bf u};\lambda\right]$ $\displaystyle={\bf
g},\quad$ $\displaystyle\text{ on }\partial\widehat{\Omega},$ (3b)

where $\widehat{\Omega}\subseteq\mathbb{R}^{n}$, with boundary (including
internal interfaces, like faults) $\partial\widehat{\Omega}$. Vector ${\bf u}$
is the unknown solution and ${\bf g}$ is boundary data. The source term ${\bf
k}$ encompasses external and internal body forces, and $\mathcal{L}$,
$\mathcal{B}$ are differential and boundary operators parameterized by
$\lambda$.

The PINN first approximates the solution to the IBVP (3) using a feed-forward,
deep neural network, i.e. we assume ${\bf u}({\bf x})\approx\mathcal{N}({\bf
x};\theta)$ and consider $N_{\partial}$ many boundary points $\\{{\bf
x}^{i}_{\partial},{\bf g}^{i}\\}_{i=1}^{N_{\partial}}$ where ${\bf g}^{i}={\bf
g}({\bf x}^{i}_{\partial})$ for each $1\leq i\leq N_{\partial}$, and a set of
$N_{{\widehat{\Omega}}}$ internal collocation points $\\{{\bf
x}_{{\widehat{\Omega}}}^{i},{\bf k}^{i}\\}_{i=1}^{N_{\widehat{\Omega}}}$ with
${\bf k}^{i}={\bf k}({\bf x}^{i}_{\widehat{\Omega}})$ for $1\leq i\leq
N_{\widehat{\Omega}}$. Throughout this work, collocation points are generated
from uniform random sampling. From (3) we construct an objective function
$MSE$ (based on the mean-square error) given by

$\displaystyle
MSE(\theta)=MSE_{{\widehat{\Omega}}}(\theta)+MSE_{\partial}(\theta),$ (4)

where

$\displaystyle MSE_{{\widehat{\Omega}}}(\theta)$
$\displaystyle=\frac{1}{N_{{\widehat{\Omega}}}}\sum_{i=1}^{N_{{\widehat{\Omega}}}}\big{|}\mathcal{L}[\mathcal{N};\lambda]({\bf
x}^{i}_{{\widehat{\Omega}}};\theta)-{\bf k}^{i}\big{|}^{2},$ (5a)
$\displaystyle MSE_{\partial}(\theta)$
$\displaystyle=\frac{1}{N_{\partial}}\sum_{i=1}^{N_{\partial}}\big{|}\mathcal{B}[\mathcal{N};\lambda]({\bf
x}^{i}_{\partial};\theta)-{\bf g}^{i}\big{|}^{2},$ (5b)

are the contributions to the total loss from the PDE (3a) and boundary
conditions (3b), respectively, corresponding to misfit of the network with
data. Because our objective function (4) consists of multiple competing loss
functions it is referred to as a multi-objective loss. Solving (3) is done by
minimizing (4) with respect to the network parameters $\theta$, typically done
via an optimization algorithm that approximates the set of solution parameters
$\theta^{*}$ such that
$\theta^{*}=\operatorname*{arg\,min}_{\theta}MSE(\theta)$. Figure 1 gives a
schematic representation of the objective function (4). Note that for the rest
of this work with omit the notation specifying MSE dependency on network
parameters $\theta$ and assume this is understood.

Figure 1: A schematic of the PINN framework for solving the general boundary
value problem from equation (3). Displacement approximation network
$\mathcal{N}$ is trained on interior and boundary subdomains which are
governed by operators $\mathcal{L}$ and $\mathcal{B}$, respectively.

##  3 Learning Problems for Earthquakes on Rate-and-State Faults

In order to best illustrate the PINN computational set-up, we begin by
considering a relevant problem in 1D, in order to introduce our methodology in
a simplified framework. Here we focus on the forward problem, with details of
the inverse problem included in the subsequent section on the 2D application
problem.

###  3.1 1D Illustration

The solid Earth is governed by momentum balance and a constitutive relation
defining the material rheology; in this work we assume elastic material
properties. These assumptions give rise to the elastic wave equation, namely

$u_{tt}=c^{2}u_{xx}+s(x,t),\quad x\in[0,1],t\geq 0,$ (6)

where $u(x,t)$ is the Earth’s material displacement, $c=\sqrt{\mu/\rho}$ is
the wavespeed where $\rho$ and $\mu$ are the the density and shear modulus,
and $s(x,t)$ accounts for external and/or body forces. Assumed initial and
boundary conditions are given by

$\displaystyle u(x,0)$ $\displaystyle=u_{0}(x),$ (7a) $\displaystyle
u_{t}(x,0)$ $\displaystyle=v_{0}(x)$ (7b) $\displaystyle-\mu u_{x}$
$\displaystyle=F+g_{0}(t),\quad$ $\displaystyle x=0,$ (7c) $\displaystyle
Zu_{t}+\mu u_{x}$ $\displaystyle=g_{1}(t),\quad$ $\displaystyle x=1,$ (7d)

where $Z=\sqrt{\mu\rho}$ is the shear impedance. If boundary data
$g_{0}(t)=g_{1}(t)=0$, equation (7c) corresponds to the requirement that fault
shear stress be equal to frictional strength $F$, and (7d) allows waves to
freely exit the domain.

In this work we consider rate-and-state dependent friction (RSF), an
experimentally-motivated, nonlinear friction law, used in the majority of
modern earthquake simulations for its ability to reproduce a wide range of
observed seismic and aseismic behaviors [10, 45, 36]. In this context, the
frictional strength is given by

$F=\bar{\sigma}_{n}f(V,\psi),$ (8)

where $V$ is the slip rate (jump in velocity across a fault interface),
$\bar{\sigma}_{n}$ is the effective normal stress (normal stress minus pore
fluid pressure), and $f$ is a friction coefficient given by

$f(V,\psi)=a\ln\left(\frac{V}{V_{0}}\right)+\psi,$ (9)

where $V_{0}$ is a reference slip rate and material parameter $a$ is the
“direct effect” [50]. State variable $\psi$ evolves from some initial value
$\psi_{0}$ according to its own evolution law

$\displaystyle\frac{d\psi}{dt}$ $\displaystyle=G(V,\psi)+h(t),$ (10a)
$\displaystyle\psi(0)$ $\displaystyle=\psi_{0}$ (10b)

where we assume the aging law, i.e.
$G(V,\psi)=(bV_{0}/D_{c})\exp{(\frac{f_{0}-\psi}{b}-\frac{|V|}{V_{0}}})$,
which allows state to evolve even in the absence of slip [36], and $h(t)$
allows for the incorporation of additional source data. Here $D_{c}$ is the
characteristic slip distance, $f_{0}$ is a reference friction coefficient for
sliding at speed $V_{0}$ and material parameter $b$ captures time-dependent
“evolution” effects [50].

The initial boundary value problem (IBVP) formed by the PDE (6), specified
initial/boundary conditions (7) and state-evolution equation (10a) correspond
to specific differential and boundary operators $\mathcal{L}$ and
$\mathcal{B}$ defined in (3), for example, here
$\mathcal{L}=\partial^{2}/\partial t^{2}-c^{2}\partial^{2}/\partial x^{2}$.
Assuming all parameters of the IBVP are known, one can solve the forward
problem for the unknown material displacement $u(x,t)$ and state variable
$\psi(t)$. We define neural networks $\mathcal{N}(x,t;\theta)\approx u(x,t)$
and $\mathcal{N}^{\psi}(t;\theta_{\psi})\approx\psi(t)$ to approximate the
exact solutions. The PDE, boundary and initials conditions all give rise to a
term in the objective function. We refer to these terms as the component loss
functions and define them as

$\displaystyle MSE_{\Omega}$
$\displaystyle=\frac{1}{N_{\Omega}}\sum_{i=1}^{N_{\Omega}}\big{|}\mathcal{N}_{tt}-c^{2}\mathcal{N}_{xx}-s\big{|}^{2},$
(11a) $\displaystyle MSE_{0}$
$\displaystyle=\frac{1}{N_{0}}\sum_{i=1}^{N_{0}}\big{|}-\mu\mathcal{N}_{x}-F-g_{0}\big{|}^{2},$
(11b) $\displaystyle MSE_{1}$
$\displaystyle=\frac{1}{N_{1}}\sum_{i=1}^{N_{1}}\big{|}Z\mathcal{N}_{t}+\mu\mathcal{N}_{x}-g_{1}\big{|}^{2},$
(11c) $\displaystyle MSE_{u_{0}}$
$\displaystyle=\frac{1}{N_{u_{0}}}\sum_{i=1}^{N_{u_{0}}}\big{|}\mathcal{N}-u_{0}\big{|}^{2},$
(11d) $\displaystyle MSE_{v_{0}}$
$\displaystyle=\frac{1}{N_{v_{0}}}\sum_{i=1}^{N_{v_{0}}}\big{|}\mathcal{N}_{t}-v_{0}\big{|}^{2},$
(11e)

where $N_{\Omega}$, $N_{0}$, $N_{1}$, $N_{u_{0}}$, $N_{v_{0}}$ correspond to
the number of randomly distributed collocation points in the relevant
subdomains, i.e. the domain interior, left boundary ($x=0$), right boundary
($x=1$), and at $t=0$ (for initial displacements and velocities),
respectively. In addition to (11) we have the contributions to the total loss
from the state evolution equation, namely

$\displaystyle MSE_{\psi}$
$\displaystyle=\frac{1}{N_{\psi}}\sum_{i=1}^{N_{\psi}}\big{|}\mathcal{N}^{\psi}_{t}-G(V,\psi)-h\big{|}^{2},$
(12a) $\displaystyle MSE_{\psi_{0}}$
$\displaystyle=\frac{1}{N_{\psi_{0}}}\sum_{i=1}^{N_{\psi_{0}}}\big{|}\mathcal{N}^{\psi}-\psi_{0}|^{2},$
(12b)

so that the objective function we seek to minimize is given by
$MSE=\sum_{\xi\in\chi}MSE_{\xi}$, where subscript
$\chi=\\{\Omega,0,1,u_{0},v_{0},\psi,\psi_{0}\\}$ so that the sum involves
contributions from all subdomains. Minimizing the $MSE$ over the network
parameters of $\mathcal{N}$ and $\mathcal{N}^{\psi}$ define proxy models which
produce approximations to $u$ and $\psi$.

An important detail about the network training involved in minimizing the
objective function $MSE$ defined in the previous paragraph is that the
inclusion of the state evolution means we are now solving a coupled system of
PDE. The coupling occurs at the fault $x=0$, where two conditions (fault
friction (7d) and state evolution (10a)) are now enforced. Because of this
coupling, the training networks $\mathcal{N}$ and $\mathcal{N}^{\psi}$ both
appear in the fault loss (11b) as well as in the state evolution loss (12a).
This interconnectivity may cause training to favor accuracy in the
displacement network over accuracy in the state network. In particular, high
temporal variations in state evolution, variations in scale, and/or network
architecture may drive the training step to favor the displacement
approximation. We found it helpful to isolate the state and displacement
networks during the backpropogation step. To do this, we define two objective
functions: The state objective function consists of loss components in
equation (12) while the displacement objective function uses loss components
from equation (11). Each training iteration then requires two optimization
steps to update both state and displacement networks. This setup ensures that
each network is only updated by one associated objective function. We found
this approach to improve our network’s training speed and resulted in better
approximations of the state variable. The particulars of this training
approach are illustrated in Algorithm 1 which provides the pseudocode for
solving the 1D IBVP given by equations (6) and (7) where state evolves
according to (10a).

In order to verify our computational framework for the 1D problem, we generate
a known, manufactured solution, which enables direct comparison with the
approximate solution produced by the neural network. This approach, known as
the method of manufactured solutions [44], assumes a particular analytic
solution $u^{e}$ and derives consistent source terms and boundary data for the
IBVP, without changing the underlying PDE and boundary operators. We
manufacture $u^{e}$ and $\psi^{e}$ to be a solution to the IBVP and aging law,
respectively, such that

$\displaystyle u^{e}(x,t)$ $\displaystyle=\tanh{\big{(}0.5(x-ct+1)\big{)}},$
(13) $\displaystyle\psi^{e}(t)$ $\displaystyle=-\frac{\mu
u^{e}_{x}(0,t)}{\bar{\sigma}_{n}}-a\ln\bigg{(}\frac{2u^{e}_{t}(0,t)}{V_{0}}\bigg{)}$
(14)

This choice of manufactured solution defines the initial data $u_{0}$ and
$v_{0}$, the source term $s$, and all boundary data, namely

$\displaystyle s$ $\displaystyle=u^{e}_{tt}-c^{2}u^{e}_{xx},\quad$
$\displaystyle\text{on }\Omega\times[0,T],$ (15a) $\displaystyle g_{0}$
$\displaystyle=0,\quad$ $\displaystyle\text{at }x=0,$ (15b) $\displaystyle
g_{1}$ $\displaystyle=Zu^{e}_{t}+\mu u^{e}_{x},\quad$ $\displaystyle\text{at
}x=1,$ (15c) $\displaystyle h$
$\displaystyle=\psi^{e}_{t}-G(2u^{e}_{t},\psi^{e}),\quad$
$\displaystyle\text{at }x=0,$ (15d)

where, for the 1D problem, $V(t)=2u_{t}(0,t)$, and corresponds to a commonly
chosen exact solution [11, 21].

Parameter values used for all studies in this work are given in Table 1. We
provide Figure 2 to illustrate how the independent networks $\mathcal{N}$ and
$\mathcal{N}^{\psi}$ are used to build each component loss function and
Algorithm 1 provides an outline of the implemented code which is publicly
available on GitHub [13].

Figure 2: A schematic of the PINN framework for solving the 1D IBVP (6) (7)
with rate and state friction defined by (9) and (10a). Displacement network
$\mathcal{N}$ and state network $\mathcal{N}^{\psi}$ are trained separately by
defining two objective functions which are differentiated by red and blue
nodes, respectively, in the final layer. At each training iteration,
$\mathcal{N}$ is updated using component losses $MSE_{\Omega},MSE_{0},$ and
$MSE_{1}$ while $\mathcal{N}^{\psi}$ is updated using just $MSE_{\psi}$ as an
objective function.

1:$\mathcal{N}\leftarrow$ NeuralNetwork($x,t;\theta$) $\triangleright$
initialize network parameters $\theta$ for $\mathcal{N}$

2:$\mathcal{N}^{\psi}\leftarrow$ NeuralNetwork($t;\theta_{\psi}$)
$\triangleright$ initialize network parameters $\theta_{\psi}$ for
$\mathcal{N}^{\psi}$

3:$\mathcal{N}$opt $\leftarrow$ SetOptimizer($\theta$) $\triangleright$
Instantiate an optimizer which only tracks parameters of $\mathcal{N}$

4:$\mathcal{N}^{\psi}$opt$\leftarrow$ SetOptimizer($\theta_{\psi}$)
$\triangleright$ Instantiate an optimizer which only tracks parameters of
$\mathcal{N}^{\psi}$

5:

6:

7:$\widetilde{\mathcal{N}}\leftarrow u_{0}(x)+tv_{0}(x)+t^{2}\mathcal{N}(x,t)$
$\triangleright$ Employ hard enforcement of initial conditions

8:$\widetilde{\mathcal{N}}^{\psi}\leftarrow\psi_{0}+t\mathcal{N}^{\psi}(t)$

9:

10:for $\xi$ in $\chi$ do $\triangleright$ Loop over each subdomain and define
the relevant component loss

11: function MSEξ($x,t$)

12: output $\leftarrow$
Condition$\big{[}\xi\big{]}$($\widetilde{\mathcal{N}}(x,t)$,
$\widetilde{\mathcal{N}}^{\psi}(t)$) $\triangleright$ Conditions specified by
equations (6), (7), (8), and (10a)

13: data $\leftarrow$ SourceData$\big{[}\xi\big{]}$($x,t$) $\triangleright$
Source terms are determined in equation (15)

14: return MeanSquareError(output, data)

15:

16:for i=1 to training_iterations do

17: Loss${}_{\mathcal{N}}\leftarrow$ 0

18: Loss${}_{\mathcal{N}^{\psi}}\leftarrow$ 0

19:

20: for $\xi$ in $\chi\setminus\\{\Psi\\}$ do

21: $x_{\xi},t_{\xi}\leftarrow$rand($N_{\xi},\xi$) $\triangleright$ Generate
$N_{\xi}$ randomly sampled points from subdomain $\xi$

22: Loss${}_{\mathcal{N}}\leftarrow$ LossN \+ MSEξ($x_{\xi},t_{\xi}$)

23:

24: $x_{\Psi},t_{\Psi}\leftarrow$rand($N_{\Psi},{\Psi}$)

25: Loss${}_{\mathcal{N}^{\psi}}\leftarrow$
MSE${}_{f_{\psi}}$($x_{\Psi},t_{\Psi}$)

26:

27: $\nabla_{\psi}\leftarrow$ grad(Loss${}_{\mathcal{N}^{\psi}}$,
[$\theta_{\psi}$]) $\triangleright$ Compute network gradient

28: $\theta_{\psi}\leftarrow$ $\mathcal{N}^{\psi}$opt.Step($\nabla_{\Psi}$)
$\triangleright$ Update network weights using the optimization algorithm

29:

30: $\nabla_{\mathcal{N}}\leftarrow$ grad(LossN, [$\theta$])

31: $\theta\leftarrow$ $\mathcal{N}$opt.Step($\nabla_{\mathcal{N}}$)

Algorithm 1 Multi-network training for 1D forward problem with RSF.

Here we take $\mathcal{N}:\mathbb{R}^{2}\to\mathbb{R}$ and state variable
network $\mathcal{N}^{\psi}:\mathbb{R}\to\mathbb{R}$ to both be fully-
connected, feed-forward networks with three hidden layers, 64 neurons per
layer and use hyperbolic tangent activation functions. The loss function is
optimized using L-BFGS[4], a quasi-Newton optimization algorithm. Minimization
is done over 10 training iterations where each iteration is trained on a
random sampling of $N=100$ interior points and 50 boundary points ($N_{b}=25$
points per boundary). Network weights are initialized using uniform Xavier
initialization [17]. In Figure 3 we show learned approximations for
displacement and state compared to their respective exact solutions. These
results suggest that our neural network parameters can solve the 1D forward
problem to reasonable accuracy. To further explore performance however, and
since we are primarily interested in higher dimensional settings that enable
the inversion of depth-dependent parameters, we now move to the 2D
formulation.

Figure 3: Comparison of results from 1D illustration showing the (a)
displacement network approximation $\mathcal{N}$ with (b) manufactured
solution $u^{e}$. Additionally, the (c) state approximation network
$\mathcal{N}^{\psi}$ is plotted against the manufactured state $\psi^{e}$
along with their absolute error $|\mathcal{N}^{\psi}(t)-\psi^{e}(t)|$.
Absolute displacement error was averaged over 1000 randomly sampled points and
measured to be $|\mathcal{N}-u^{e}|_{\text{avg}}=1.57e-5$.

###  3.2 2D Application in Seismic Faulting

The 2D problem is obtained by first considering a bounded spatial domain in
$\mathbb{R}^{3}$ defined by
$(x,y,z)\in[-L_{x},L_{x}]\times[-L_{y},L_{y}]\times[0,L_{z}]$ where $z$ is
taken to be positive downward. We assume antiplane shear deformation, taking
displacements $u_{x},u_{z}$ to be zero, and assume the out-of-plane
displacement $u=u_{y}$ is independent of $y$. Momentum balance for an elastic
solid then gives rise to the elastodynamic wave equation in 2D, namely,

$\displaystyle u_{tt}=c^{2}\Delta
u+S(x,z,t),\quad(x,z)\in[-L_{x},L_{x}]\times[0,L_{z}],\,\,t\geq 0,$ (16)

where $\Delta=\left(\dfrac{\partial^{2}}{\partial
x^{2}}+\dfrac{\partial^{2}}{\partial z^{2}}\right)$ is the 2D Laplacian. Here
the shear wavespeed $c=\sqrt{\mu/\rho}$ as in the 1D case, and $S$ comprises
body forces. Assumed initial conditions are given by

$\displaystyle u(x,z,0)$ $\displaystyle=u_{0}(x,z),$ (17a) $\displaystyle
u_{t}(x,z,0)$ $\displaystyle=v_{0}(x,z).$ (17b)

We assume $z=0$ corresponds to Earth’s surface which we take to be traction
free. An RSF vertical fault is embedded at the interface $x=0$, corresponding
to interface conditions

$\displaystyle\tau^{+}$ $\displaystyle=-\tau^{-}$ (18a)
$\displaystyle\tau^{\pm}$ $\displaystyle=F(V^{\pm},\psi).$ (18b)

Here the traction $\tau$ at a boundary or interface is defined with respect to
the outward pointing unit normal $\bf n$ by $\tau={\bf n}\cdot\mu\nabla u$,
and $\tau^{\pm}=\tau(0^{\pm},z,t)$. Equation (18a) corresponds to the
requirement that the traction be “equal and opposite” across the fault and
(18b) is a requirement that fault shear stress be equal to frictional strength
$F$ as in the 1D case. As in the 1D illustration, $F$ is a nonlinear function
of the slip rate $V^{\pm}(z,t)=\dot{u}(0^{\pm},z,t)-\dot{u}(0^{\mp},z,t)$ and
an empirical state variable $\psi(z,t)$ the particular form of which is given
by (8). At the remote boundaries $(x,z)=(\pm L_{x},L_{z})$ we assume the
boundaries are non-reflecting. Next we make the assumption that the
displacement field $u(x,z,t)$ is antisymmetric about the fault interface
$x=0$, which automatically satisfies (18a), so that we may restrict our focus
to one side of the fault as in Erickson and Dunham [11] and reduce the
computational redundancy. Our final spatial domain is then
$\Omega=[0,L_{x}]\times[0,L_{z}]$ where the fault interface now becomes a
boundary. For generality (and to aid in later verification steps) we state the
final boundary conditions on $\partial\Omega$ using generic data
$g_{f},g_{s},g_{r},g_{d}$ for the fault, surface, remote, and depth
boundaries, respectively, namely

$\displaystyle\tau$ $\displaystyle=F(V)+g_{f}(z,t),\quad$ $\displaystyle x=0,$
(19a) $\displaystyle\tau$ $\displaystyle=g_{s}(x,t),\quad$ $\displaystyle
z=0,$ (19b) $\displaystyle Zu_{t}+\tau$ $\displaystyle=g_{r}(z,t),\quad$
$\displaystyle x=L_{x},$ (19c) $\displaystyle Zu_{t}+\tau$
$\displaystyle=g_{d}(z,t),\quad$ $\displaystyle z=L_{z},$ (19d)

although the Earth’s free surface and non-reflecting conditions are the far
field boundaries correspond to $g_{s}=g_{r}=g_{d}=0$. Here $\tau=\tau^{+}$,
$V=V^{+}=2u_{t}^{+}$ and $Z=\sqrt{\mu\rho}$ is the shear impedance as in the
1D case.

For this 2D application problem we are primarily interested in the RSF
parameter $a-b$, which describes the velocity-dependence of friction at steady
state. Positive values (i.e. $a>b$) correspond to stable sliding, while
negative values correspond to frictional instabilities. Knowledge of the
depth-distribution of $a-b$ is therefore of central importance to
understanding slip behavior and deformation [3]. In other words,
heterogeneities along fault interfaces can be characterized at least in part
by velocity-weakening frictional behavior ($a-b<0$) indicating that seismic
rupture may nucleate and easily propagate, while stable regions are
characterized by velocity-strengthening frictional behavior ($a-b>0$) that
inhibits the sliding at seismogenic speeds [46]. In order to focus our study
on the inference of the depth-dependency of $a-b$, we assume that the state
variable is at steady state, namely
$\psi=f_{0}+b\ln\left(\frac{V_{0}}{V}\right)$, so that the friction
coefficient reduces to the rate-dependent form

$f(V)=f_{0}+\alpha\ln\left(\frac{V}{V_{0}}\right),$ (20)

where we have introduced $\alpha=a-b$. We assume that $\alpha$ is piecewise
linear with depth, defining a shallow seismogenic zone, namely,

$\displaystyle\alpha(z)=\begin{cases}\alpha_{\text{min}}&0<z<H\\\
(z-H)*((\alpha_{\text{max}}-\alpha_{\text{min}})/D)+\alpha_{\text{min}}&H\leq
z\leq H+D\\\ \alpha_{\text{max}}&H+D<z,\end{cases}$ (21)

where $H$ defines the seismogenic depth, $\alpha_{\text{min}}$ and
$\alpha_{\text{max}}$ are constants defining the minimum and maximum values
assumed by $\alpha$, and $D$ the transition distance. In this work we consider
PINN solutions to both the forward and inverse problems. For forward problem
we solve for the unknown displacements $u$ in the IBVP (3) assuming
$\alpha(z)$ has been specified. For the inverse problems, in addition to
solving for $u$ we also infer the friction parameter $\alpha$. In this latter
case, (21) sets the data for the boundary conditions, but $\alpha$ is not
assumed to be known a priori when enforcing the frictional interface
condition, as will be described in the next section.

###  3.3 Forward and Inverse Problems for the 2D Application

Governing equation (16) along with initial and boundary conditions (17), (19)
provide specifics of the loss terms (5) that define the PINN. However, there
is more than one way to formulate the associated learning problem. First, one
can consider either a forward or inverse problem [43]. As this work is
concerned with generating approximations to both forward and inverse problems,
we will refer to the primal solution as the network approximation
$\mathcal{N}$ to the displacement $u$ specified by the IBVP (3). The forward
problem is only concerned with approximating a primal solution and requires
that all model parameters $\lambda$ are known a priori. A solution to the
inverse problem aims to approximate the primal solution while also being
tasked with learning a set of system parameters. If the $\lambda$ in (3) are
known, the forward problem is solved and the problem reduces to an
unsupervised learning task [7]. In the inverse problem, however, the PINN has
the same loss function (5) but with minor changes: Instead of knowing the
system parameters $\lambda$, we establish them as trainable networks, which we
detail shortly.

In addition to specifying whether a forward or inverse problem is being
solved, in either case one must also choose between soft or hard enforcement
of initial conditions [33, 34, 26, 7]. Soft enforcement uses loss terms to
learn boundary data whereas hard enforcement encodes boundary data into a
trial function to satisfy the conditions exactly. We proceed with details
concerning each of these below.

####  3.3.1 Forward Problem

As in the 1D case, we begin by supposing that
$\mathcal{N}(x,z,t;\theta)\approx u(x,z,t)$ for neural network $\mathcal{N}$
and solution $u$ of the IBVP (16), (17), (19). The 2D component loss functions
are therefore given by

$\displaystyle MSE_{\Omega}$
$\displaystyle=\frac{1}{N_{\Omega}}\sum_{i=1}^{N_{\Omega}}\big{|}\mathcal{N}_{tt}-c^{2}\Delta\mathcal{N}-S\big{|}^{2},$
(22a) $\displaystyle MSE_{f}$
$\displaystyle=\frac{1}{N_{f}}\sum_{i=1}^{N_{f}}\big{|}-\mu\mathcal{N}_{x}-\bar{\sigma}_{n}f-g_{f}\big{|}^{2},$
(22b) $\displaystyle MSE_{s}$
$\displaystyle=\frac{1}{N_{s}}\sum_{i=1}^{N_{s}}\big{|}-\mu\mathcal{N}_{z}-g_{s}\big{|}^{2},$
(22c) $\displaystyle MSE_{r}$
$\displaystyle=\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}\big{|}Z\mathcal{N}_{t}+\mu\mathcal{N}_{x}-g_{r}\big{|}^{2},$
(22d) $\displaystyle MSE_{d}$
$\displaystyle=\frac{1}{N_{d}}\sum_{i=1}^{N_{d}}\big{|}Z\mathcal{N}_{t}+\mu\mathcal{N}_{z}-g_{d}\big{|}^{2}$
(22e) $\displaystyle MSE_{u_{0}}$
$\displaystyle=\frac{1}{N_{u_{0}}}\sum_{i=1}^{N_{u_{0}}}\big{|}\mathcal{N}-u_{0}\big{|}^{2},$
(22f) $\displaystyle MSE_{v_{0}}$
$\displaystyle=\frac{1}{N_{v_{0}}}\sum_{i=1}^{N_{v_{0}}}\big{|}\mathcal{N}_{t}-v_{0}\big{|}^{2},$
(22g)

where, as in the 1D case, summations are over randomly distributed collocation
points in the interior and domain boundaries. The 2D objective function we
seek to minimize is given by $MSE=\sum_{\xi\in\chi}MSE_{\xi}$, where (as in
the 1D case), subscript $\chi=\\{\Omega,f,s,r,d,u_{0},v_{0}\\}$ implies that
the sum is taken over all relevant subdomains.

####  3.3.2 Inverse Problem

While the primal solution to an IBVP is desirable, replacing an assumed model
parameter (which often comes with much uncertainty) with a trainable network
could leverage real-world data to physically constrain parameter
distributions. The inverse problem is derived by making a slight modification
to the forward problem (22). As we are interested in learning the fault
friction parameter $\alpha$, we modify the loss component (22b) to include a
trainable network. Because $\alpha$ is depth-variable it must be explicitly
trained on random collocation points along the fault. So in addition to the
first network $\mathcal{N}(x,z,t;\theta)\approx u(x,z,t)$, we introduce a
secondary network $\mathcal{N}^{\alpha}(z;\theta_{\alpha})\approx\alpha(z)$.
Only one loss component needs be changed from the forward problem (22),
namely, component (22b), which enforces the fault friction condition. A
modified friction coefficient is considered, defined by

$\displaystyle\tilde{f}=f_{0}+\mathcal{N}^{\alpha}(z)\ln{\bigg{[}\frac{2\mathcal{N}_{t}(z,t)}{V_{0}}\bigg{]}},$
(23)

which gives rise to the modified fault loss

$\displaystyle\overline{MSE}_{f}$
$\displaystyle=\frac{1}{N_{f}}\sum_{i=1}^{N_{f}}\big{|}-\mu\mathcal{N}_{x}-\bar{\sigma}_{n}\tilde{f}-g_{f}\big{|}^{2}.$
(24)

Replacing loss component (22b) in the forward problem with (24) defines the
inverse problem. We give a diagram of the fault network for the inverse
problem in Figure 4 showing how the fault loss is used to update two networks
during training. The inverse problem is trained over the same subdomains as
the forward problem (22) but the depth coordinate for fault training data is
passed to both $\mathcal{N}$ and $\mathcal{N}^{\alpha}$.

Figure 4: Network diagram for enforcing the stress condition along the fault
while also learning the depth-dependent friction parameter $\alpha$. During
training, the $z$-coordinate of a training point is passed to both
$\mathcal{N}$ and $\mathcal{N}^{\alpha}$.

####  3.3.3 Soft vs. Hard Enforcement of Boundary Conditions

Both the forward and inverse problems require that we specify initial and
boundary conditions, which may be enforced in two possible ways. Soft
enforcement is done by penalizing the objective function, as the network
output need not satisfy the condition exactly (as was done in the 1D example
of the previous section). Up until now, we have presented soft enforcement, as
given in (22), which shows component losses corresponding to the boundary and
initial conditions. Soft enforcement does not give us any guarantee regarding
accuracy of the condition being enforced and each additional loss term in the
objective function increases the complexity of the optimization landscape.

Alternatively, we may enforce initial/boundary conditions in a manner that
reduces the complexity of the network training space by encoding such
conditions into the network architecture. This approach, known as hard
enforcement, requires that we specify a trial function which automatically
satisfies the initial/boundary conditions. Here we apply this technique to
enforce initial conditions, which reduces the number of constraints on the
system. Let $\mathcal{N}({\bf x},t;\theta)$ be a feed-forward neural network
and define the trial function $\mathcal{\widetilde{N}}$ to be

$\displaystyle\mathcal{\widetilde{N}}({\bf x},t;\theta)=u_{0}({\bf
x})+tv_{0}({\bf x})+t^{2}\mathcal{N}({\bf x},t;\theta)$ (25)

where ${\bf x}=(x,z)$, so that $\mathcal{\widetilde{N}}$ satisfies both
initial conditions (17) exactly and is trained by training the network
$\mathcal{N}$. If $\mathcal{\widetilde{N}}$ is used in place of $\mathcal{N}$
in the formulation of the objective function, there is no longer a need to
include loss terms for initial displacements (22f), or initial velocities
(22g) so the network can be trained on a less restrictive set of conditions.

##  4 2D Verification, Validation and Applications

When computational methods for physical problems are used to address science
questions, verification is an essential first step to ensure credible results
[20, 12]. While validation with observational data is the focus of future
work, we must first verify that our physics-informed deep learning framework
is able to solve both forward and inverse problems to reasonable accuracy.

Parameter | Value
---|---
$L_{x},L_{z}$ | 25 $\mathrm{km}$
$H$ | 12 $\mathrm{km}$
$D$ | 5 $\mathrm{km}$
$\mu$ | 32 $\mathrm{GPa}$
$\rho$ | 2.67 $\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$
$\alpha_{\text{min}}$ | -0.005
$\alpha_{\text{max}}$ | 0.015
$f_{0}$ | 0.6
$\bar{\sigma}_{n}$ | 50 $\mathrm{MPa}$
$V_{0}$ |  $10^{-6}$ $\mathrm{m}\text{/}\mathrm{s}$
$D_{c}$ | 2 $\mathrm{m}$
Table 1: Parameter values used in the manufactured solution tests.

###  4.1 Verification with the Method of Manufactured Solutions

As in the 1D case, we verify the PINN framework by generating a known,
manufactured solution $u^{e}$ which solves the IBVP (3), (17), (19). We take

$\displaystyle u^{e}(x,z,t)=\tanh{\big{(}(x+z+ct)/20\big{)}},$ (26)

which defines the initial data $u_{0}$ and $v_{0}$, the source term $S$, and
all boundary data, which for the forward problem is

$\displaystyle S$ $\displaystyle=u^{e}_{tt}-c^{2}\Delta u^{e},\quad$
$\displaystyle\text{on }\Omega\times[0,T],$ (27a) $\displaystyle g_{f}$
$\displaystyle=-\mu u^{e}_{x}-\bar{\sigma}_{n}f\big{(}2u^{e}_{t}\big{)},\quad$
$\displaystyle\text{at }x=0,$ (27b) $\displaystyle g_{s}$ $\displaystyle=-\mu
u^{e}_{z},\quad$ $\displaystyle\text{at }z=0,$ (27c) $\displaystyle g_{r}$
$\displaystyle=Zu^{e}_{t}+\mu u^{e}_{x},\quad$ $\displaystyle\text{at
}x=L_{x},$ (27d) $\displaystyle g_{d}$ $\displaystyle=Zu^{e}_{t}+\mu
u^{e}_{z},\quad$ $\displaystyle\text{at }z=L_{z}.$ (27e)

The parameters used are slightly modified from the earthquake simulations of
Harvey et al. [21], given in Table 1, and we default to this parameter set
unless stated otherwise. Note that in the case of the inverse problem,
$\alpha_{\text{min}}$ and $\alpha_{\text{max}}$ define the manufactured exact
solution $\alpha^{e}(z)$, which is used only to set the data in (27b);
$\alpha(z)$ is a learned parameter through the use of the network
$\mathcal{N}^{\alpha}$. In all scenarios we consider, we take the primal
network $\mathcal{N}:\mathbb{R}^{3}\to\mathbb{R}$ and friction network
$\mathcal{N}^{\alpha}:\mathbb{R}\to\mathbb{R}$ to both be fully-connected,
feed-forward networks with three hidden layers and 128 neurons per layer. The
wave equation (16) requires displacements to be sufficiently smooth with
respect to the input but no such requirement is imposed on the friction
parameters. Thus we use hyperbolic tangent activation functions for the
displacement network but use linear units in the friction network. We found
that the friction network performed well when we used both Rectified Linear
Unit (ReLU) and Sigmoid Linear Unit (SiLU) together. We use ReLU activation
for the outer hidden layers and a SiLU activation on the interior layer. The
code used to perform the numerical simulations in this work is developed in
Python using Pytorch and is publicly available [13]. Additionally, we provide
Algorithm 2 as a sketch of the implemented code used to solve the inverse
problem with hard enforcement of initial conditions.

1:$\mathcal{N}\leftarrow$ NeuralNetwork($x,z,t;\theta$) $\triangleright$
$\theta$ initial network parameters for $\mathcal{N}$

2:$\mathcal{N}^{\alpha}\leftarrow$ NeuralNetwork($z;\theta_{\alpha}$)
$\triangleright$ $\theta_{\alpha}$ initial network parameters for
$\mathcal{N}^{\alpha}$

3:$\theta\leftarrow[\theta,\theta_{\alpha}]$

4:

5:$\widetilde{\mathcal{N}}\leftarrow
u_{0}(x,z)+tv_{0}(x,z)+t^{2}\mathcal{N}(x,z,t)$ $\triangleright$ Employ hard
enforcement of initial conditions

6:

7:for $\xi$ in $\chi$ do $\triangleright$ Loop over each subdomain and define
the relevant component loss

8: function MSEξ(x, z, t)

9: output $\leftarrow$
Condition$\big{[}\xi\big{]}$($\widetilde{\mathcal{N}}(x,z,t)$,
$\mathcal{N}^{\alpha}(z)$)

10: data $\leftarrow$ SourceData$\big{[}\xi\big{]}$($x,z,t$)

11: return MeanSquareError(output, data)

12:

13:for i=1 to training_iterations do

14: Loss $\leftarrow$ 0

15:

16: for $\xi$ in $\chi$ do

17: ${\bf x}_{\xi}\leftarrow$rand($N_{\xi},\xi$) $\triangleright$ Generate
$N_{\xi}$ randomly sampled points from subdomain $\xi$

18: Loss $\leftarrow$ Loss + $MSE_{\xi}$(${\bf x}_{\xi}$)

19:

20: $\nabla\leftarrow$ grad(Loss, [$\theta$])

21: $\theta\leftarrow$ OptimizerStep($\nabla$) $\triangleright$ Update network
weights using the optimization algorithm

Algorithm 2 Multi-network training for solving the inverse problem

$\mathcal{N}$ and (in the case of the inverse problem) $\mathcal{N}^{\alpha}$
are trained by minimizing the MSE (the sum of component MSE given in (22))
where all loss functions are optimized using L-BFGS [4], a quasi-Newton
optimization algorithm. Minimization is done over 30 training iterations where
each iteration is trained on a random sampling of $N=400$ interior points and
400 boundary points ($N_{b}=100$ points per boundary). Additionally, if soft
enforcement of initial conditions is used, we sample 400 additional interior
points (200 for displacements and 200 for velocities). After the 30 training
iterations are complete we evaluate each of the component loss functions using
1000 randomly sampled points. Likewise we generate the same number of sample
points for $\mathcal{N}$ and $\mathcal{N}^{\alpha}$ and use their respective
manufactured solutions to compute the relative $\ell^{2}-$norm of their
errors, namely,

$\displaystyle\left\lVert\mathcal{N}-u^{e}\right\rVert^{2}_{\text{1, rel}}$
$\displaystyle=\frac{\sum_{i=1}^{N}|\mathcal{N}(x_{i},z_{i},t_{i})-u^{e}(x_{i},z_{i},t_{i})|^{2}}{\sum_{i=1}^{N}|u^{e}(x_{i},z_{i},t_{i})|^{2}},$
(28)
$\displaystyle\left\lVert\mathcal{N_{\alpha}}-\alpha^{e}\right\rVert^{2}_{\text{1,
rel}}$
$\displaystyle=\frac{\sum_{i=1}^{N_{b}}|\mathcal{N}^{\alpha}(z_{i})-\alpha^{e}(z_{i})|^{2}}{\sum_{i=1}^{N_{b}}|\alpha^{e}(z_{i})|^{2}}.$
(29)

Network weights are initialized using uniform Xavier initialization [17].
Variation due to weight initialization involving randomization is accounted
for by averaging error data across $50$ trained solutions. Table 2 presents
error data for the forward and inverse problems where both hard and soft
enforcement of initial conditions are considered for each. As the table
illustrates, for both the forward and inverse problems, hard enforcement of
initial conditions are accompanied with a smaller $\ell^{2}$-errors after 30
training iterations, which is most likely due at least in part to the fact
that hard enforcement means that the network approximation satisfies the
initial conditions exactly. Also illustrated in the table, is that all
component mean-square errors (MSE) are lower for the soft enforcement of
initial conditions compared to hard enforcement. Soft enforcement uses two
additional loss components (i.e. those enforcing initial conditions) which
results in a more biased model, or a lower variance model (i.e. the MSE is
less sensitive to training data). The trade off between bias and variance
(introduced via randomness in weight initialization and collocation points)
may account for such lower MSE when using soft enforcement (which involves
lower variances) and averaging over $50$ trained models.

Errors | Forward Problem | Inverse Problem
---|---|---
|  Soft |  Hard |  Error Ratio |  Soft |  Hard |  Error Ratio
$\left\lVert\mathcal{N}-u^{e}\right\rVert_{\text{1, rel}}$ |  1.364e-02 |  7.148e-03 |  1.908e+00 |  2.931e-02 |  1.633e-02 |  1.795e+00
$\left\lVert\mathcal{N_{\alpha}}-\alpha^{e}\right\rVert_{\text{1, rel}}$ |  |  |  |  2.977e-01 |  7.220e-02 |  4.123e+00
$MSE_{PDE}$ |  1.149e-03 |  8.883e-03 |  1.293e-01 |  1.133e-03 |  6.675e-03 |  1.697e-01
$MSE_{f}$ |  2.633e-03 |  4.115e-02 |  6.399e-02 |  8.579e-03 |  1.615e-01 |  5.311e-02
$MSE_{s}$ |  8.393e-05 |  7.868e-04 |  1.067e-01 |  1.251e-04 |  6.545e-04 |  1.911e-01
$MSE_{d}$ |  8.216e-05 |  2.115e-03 |  3.885e-02 |  4.681e-04 |  7.199e-03 |  6.502e-02
$MSE_{r}$ |  1.195e-05 |  3.233e-04 |  3.698e-02 |  1.064e-04 |  8.074e-04 |  1.318e-01
$MSE_{u_{0}}$ |  1.703e-04 |  |  |  3.722e-04 |  |
$MSE_{v_{0}}$ |  6.308e-04 |  |  |  1.375e-03 |  |
Table 2: Errors in the $\ell^{2}$-norm and MSE (loss) associated with hard and
soft enforcement of initial conditions for both the forward and inverse
problem. Each term is computed over 1000 randomly sampled points and averaged
across 50 trained solutions. Error ratios are computed by dividing soft
enforcement error by the hard enforcement error.

Next we test network performance across various subdomains and resolutions.
Let $\mathcal{N}$ be a displacement network trained to solve the inverse
problem using hard enforcement of initial conditions. This time, errors are
measured with respect to the continuous $L^{2}$-norm

$\displaystyle\left\lVert\mathcal{N}-u^{e}\right\rVert^{2}_{2,\Lambda}=\int_{\Lambda}|\mathcal{N}(x,z,t)-u^{e}(x,z,t)|^{2}\,d\lambda,$
(30)

and we consider error accumulation across $\Omega$, $[0,T]$ and
$\hat{\Omega}=\Omega\times[0,T]$, the spatial, temporal and space-time
domains, respectively. The integral in (30) is approximated using Simpson’s
composite quadrature rule which converges to the exact integral with rate
$\mathcal{O}(k^{4})$, where $k$ is the quadrature grid size (i.e. the
subinterval length). To account for variations in the trained network, we
average errors over points outside of the domain of integration. For example,
the average temporal error given by

$\displaystyle\overline{\left\lVert\mathcal{N}-u^{e}\right\rVert}_{2,[0,T]}=\frac{1}{N_{T}}\sum_{k=1}^{N_{t}}\bigg{[}\int_{[0,T]}\big{|}\mathcal{N}(x_{k},z_{k},t)-u^{e}(x_{k},z_{k},t)\big{|}^{2}\,d\lambda(t)\bigg{]}^{1/2}$
(31)

measures expected error accumulation (over time) given $N_{T}$ randomly
sampled spatial points.

Figure 5 shows displacement plots for the trained network and the exact
solution at final time $T=1$ in addition to network $L^{2}$-error
approximations for displacement in space, time and space-time (along with an
$L^{2}$-error approximation of the friction parameter) for increasingly higher
grid resolution used in the numerical quadrature. The errors remain relatively
constant with decreasing subinterval size, suggesting that the quadrature
approximation is approaching the actual $L^{2}$-error, and that the PINN
approximation maintains good accuracy even when evaluated on higher-resolution
grids. Figure 6(a) shows network component loss functions against training
iteration, revealing a non-monotonic decrease across all components. Figure
6(b) illustrates convergence of the inferred parameter $\alpha=a-b$ to the
exact distribution $\alpha^{e}$. Here we show iterations 1 and 30 merely
because after the first iteration the inferred parameter performs well and
does not vary significantly.

Figure 5: (a) 2D displacement plot for a PINN trained to solve the inverse
problem using hard enforcement of initial conditions compared to (b) the
manufactured displacements. (c)$L^{2}-$errors for displacement in space, time,
and spacetime (along with $L^{2}-$errors for the friction parameter) are
computed on a uniform grid using Simpson’s rule as a quadrature. Errors are
then recorded over several mesh refinements. Figure 6: 2D inversion results
showing (a) convergence of loss components and (b) convergence of the inferred
parameter approximation.

##  5 Summary and Future Work

We have presented a computational framework for physics-informed neural
networks (PINNs) for solving the elastodynamic wave equation with a rate-and-
state frictional fault boundary in both 1D and 2D. We consider both forward
and inverse problems, with the latter obtained by extending to a multi-network
architecture in order to learn depth-dependent friction parameters alongside
deformations in the Earth’s crust. We verified the computational framework by
applying the method of manufactured solutions to probe various error
measurements. We show that in general, hard enforcement of boundary conditions
result in trained networks that better approximate displacements and friction
parameters but tend to be worse at minimizing component loss functions when
compared to soft enforcement. We found that a PINN defined by hard enforcement
of initial conditions produces reasonable approximations for displacements and
the desired friction parameter distribution. Though the network is mesh-free,
we show that successive mesh refinements of the quadrature approximation to
the $L^{2}$-error yield near constant values, suggesting that the PINN
provides a reasonable approximation to the solution even when evaluated on
increasingly finer grids. Moreover, although the network requires sufficient
training iterations to properly learn displacements, the desired friction
parameter is learned within the first couple of iterations. This suggests that
PINNs may be a highly effective tool in inferring subsurface friction
properties along faults, constrained by both physics and observational surface
data.

While the PINN is shown to perform well when learning the state variable in
1D, and inferring depth-dependency of RSF parameter $a-b$ at steady-state in
2D, we plan to explore the capabilities of the 2D framework with non-steady-
state-evolution. This extension requires additional networks in the 2D setting
in order to approximate the state variable, and two inference networks to
capture the empirical parameters $a$ and $b$ (which become separated across
governing equations), and/or networks to learn other frictional parameters,
such as $D_{c}$, whose scaling from laboratory values to actual fault zones is
the subject of many studies [e.g. 37]. Approaching such a problem may be aided
by a better understanding of the PINN dependence on problem configuration as
well as network architecture. Additionally, it would be worthwhile to
investigate methods which hybridize PINNs with traditional numerical methods
similar to the discrete time PINN in Raissi et al. [43]. And finally, with a
PINN solution that can handle learning full rate-and-state fault friction in
2D (and eventually 3D), we would be ready to compare model outcomes against
community benchmark problems concerning dynamic rupture simulations [20] and
sequences of earthquakes and aseismic slip [12].

## References

  * Atwater et al. [2005] Atwater, B.F., Satoko, M.R., Kenji, S., Yoshinobu, T., Kazue, U., Yamaguchi, D.K., 2005\.  The Orphan Tsunami of 1700: Japanese Clues to a Parent Earthquake in North America.  2 ed., University of Washington Press.  URL: http://www.jstor.org/stable/j.ctvcwnbrv.
  * Bergen et al. [2019] Bergen, K.J., Johnson, P.A., de Hoop, M.V., Beroza, G.C., 2019\.  Machine learning for data-driven discovery in solid earth geoscience.  Science 363, eaau0323.  doi:10.1126/science.aau0323.
  * Brodsky et al. [2020] Brodsky, E.E., Mori, J.J., Anderson, L., Chester, F.M., Conin, M., Dunham, E.M., Eguchi, N., Fulton, P.M., Hino, R., Hirose, T., Ikari, M.J., Ishikawa, T., Jeppson, T., Kano, Y., Kirkpatrick, J., Kodaira, S., Lin, W., Nakamura, Y., Rabinowitz, H.S., Regalla, C., Remitti, F., Rowe, C., Saffer, D.M., Saito, S., Sample, J., Sanada, Y., Savage, H.M., Sun, T., Toczko, S., Ujiie, K., Wolfson-Schwehr, M., Yang, T., 2020\.  The state of stress on the fault before, during, and after a major earthquake.  Annual Review of Earth and Planetary Sciences 48, 49–74.  doi:10.1146/annurev-earth-053018-060507.
  * Byrd et al. [1995] Byrd, R.H., Lu, P., Nocedal, J., Zhu, C., 1995\.  A limited memory algorithm for bound constrained optimization.  SIAM Journal on scientific computing 16, 1190–1208.  doi:10.1137/0916069.
  * Cai et al. [2021] Cai, S., Wang, Z., Wang, S., Perdikaris, P., Karniadakis, G.E., 2021\.  Physics-informed neural networks for heat transfer problems.  Journal of Heat Transfer 143.  doi:10.1115/1.4050542.
  * Ciarlet [2002] Ciarlet, P.G., 2002.  The finite element method for elliptic problems. volume 40.  Siam.
  * Cuomo et al. [2022] Cuomo, S., Di Cola, V.S., Giampaolo, F., Rozza, G., Raissi, M., Piccialli, F., 2022\.  Scientific machine learning through physics–informed neural networks: Where we are and what’s next.  Journal of Scientific Computing 92, 88.  doi:10.1007/s10915-022-01939-z.
  * De Ryck et al. [2022] De Ryck, T., Jagtap, A.D., Mishra, S., 2022.  Error estimates for physics informed neural networks approximating the navier-stokes equations.  arXiv preprint arXiv:2203.09346 doi:10.48550/arXiv.2203.09346.
  * De Ryck et al. [2021] De Ryck, T., Lanthaler, S., Mishra, S., 2021.  On the approximation of functions by tanh neural networks.  Neural Networks 143, 732–750.  doi:10.1016/j.neunet.2021.08.015.
  * Dieterich [1979] Dieterich, J.H., 1979.  Modeling of rock friction 1. Experimental results and constitutive equations.  Journal of Geophysical Research, [Solid Earth] 84, 2161–2168.  doi:10.1029/JB084iB05p02161.
  * Erickson and Dunham [2014] Erickson, B.A., Dunham, E.M., 2014\.  An efficient numerical method for earthquake cycles in heterogeneous media: Alternating subbasin and surface-rupturing events on faults crossing a sedimentary basin.  Journal of Geophysical Research-Solid Earth 119, 3290–3316.  doi:10.1002/2013JB010614.
  * Erickson et al. [2020] Erickson, B.A., Jiang, J., Barall, M., Lapusta, N., Dunham, E.M., Harris, R., Abrahams, L.S., Allison, K.L., Ampuero, J.P., Barbot, S., Cattania, C., Elbanna, A., Fialko, Y., Idini, B., Kozdon, J.E., Lambert, V., Liu, Y., Luo, Y., Ma, X., Mckay, M.B., Segall, P., Shi, P., van den Ende, M., Wei, M., 2020.  The community code verification exercise for simulating sequences of earthquakes and aseismic slip (SEAS).  Seismol. Res. Lett. 91, 874–890.  doi:10.1785/0220190248.
  * [13] Erickson Laboratory, .  Earthquake PINNs.  https://github.com/Thrase/EQ_pinns.  Accessed: 2023-09-16.
  * Ern and Guermond [2004] Ern, A., Guermond, J.L., 2004\.  Theory and practice of finite elements. volume 159\.  Springer.  doi:10.1007/978-1-4757-4355-5.
  * Fukushima et al. [2023] Fukushima, R., Kano, M., Hirahara, K., 2023.  Physics-informed neural networks for fault slip monitoring: simulation, frictional parameter estimation, and prediction on slow slip events in a spring-slider system.  ESS Open Archive doi:10.22541/essoar.168988460.01601423/v1.
  * Givoli [2021] Givoli, D., 2021.  A tutorial on the adjoint method for inverse problems.  Computer Methods in Applied Mechanics and Engineering 380, 113810.  doi:10.1016/j.cma.2021.113810.
  * Glorot and Bengio [2010] Glorot, X., Bengio, Y., 2010\.  Understanding the difficulty of training deep feedforward neural networks, in: Proceedings of the thirteenth international conference on artificial intelligence and statistics, JMLR Workshop and Conference Proceedings. pp. 249–256.
  * Goswami et al. [2020] Goswami, S., Anitescu, C., Chakraborty, S., Rabczuk, T., 2020\.  Transfer learning enhanced physics informed neural network for phase-field modeling of fracture.  Theoretical and Applied Fracture Mechanics 106, 102447.  doi:10.1016/j.tafmec.2019.102447.
  * Haghighat et al. [2021] Haghighat, E., Raissi, M., Moure, A., Gomez, H., Juanes, R., 2021.  A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics.  Computer Methods in Applied Mechanics and Engineering 379, 113741.  doi:10.1016/j.cma.2021.113741.
  * Harris et al. [2009] Harris, R.A., Barall, M., Archuleta, R., Dunham, E.M., Aagaard, B., Ampuero, J.P., Bhat, H., Cruz-Atienza, V., Dalguer, L., Dawson, P., Day, S., Duan, B., Ely, G., Kaneko, Y., Kase, Y., Lapusta, N., Liu, Y., Ma, S., Oglesby, D., Olsen, K., Pitarka, A., Song, S., Templeton, E., 2009\.  The SCEC/USGS dynamic earthquake rupture code verification exercise.  Seismol. Res. Lett. 80, 119–126.  doi:10.1785/gssrl.80.1.119.
  * Harvey et al. [2023] Harvey, T.W., Erickson, B.A., Kozdon, J.E., 2023.  A high-order accurate summation-by-parts finite difference method for fully-dynamic earthquake sequence simulations within sedimentary basins.  Journal of Geophysical Research: Solid Earth , e2022JB025357doi:10.1029/2022JB025357.
  * Ide et al. [2007] Ide, S., Beroza, G.C., Shelly, D.R., Uchide, T., 2007\.  A scaling law for slow earthquakes.  Nature 447, 76–79.
  * Jagtap et al. [2020] Jagtap, A.D., Kharazmi, E., Karniadakis, G.E., 2020.  Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems.  Computer Methods in Applied Mechanics and Engineering 365, 113028.  doi:10.1016/j.cma.2020.113028.
  * Jin et al. [2021] Jin, X., Cai, S., Li, H., Karniadakis, G.E., 2021.  Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations.  Journal of Computational Physics 426, 109951.  doi:https://doi.org/10.1016/j.jcp.2020.109951.
  * Karimpouli and Tahmasebi [2020] Karimpouli, S., Tahmasebi, P., 2020\.  Physics informed machine learning: Seismic wave equation.  Geoscience Frontiers 11, 1993–2001.  doi:10.1016/j.gsf.2020.07.007.
  * Karniadakis et al. [2021] Karniadakis, G.E., Kevrekidis, I.G., Lu, L., Perdikaris, P., Wang, S., Yang, L., 2021.  Physics-informed machine learning.  Nature Reviews Physics 3, 422–440.  doi:10.1038/s42254-021-00314-5.
  * Kern [2016] Kern, M., 2016.  Numerical methods for inverse problems.  John Wiley & Sons.
  * Kharazmi et al. [2019] Kharazmi, E., Zhang, Z., Karniadakis, G.E., 2019.  Variational physics-informed neural networks for solving partial differential equations.  arXiv preprint arXiv:1912.00873 doi:10.48550/arXiv.1912.00873.
  * Kharazmi et al. [2021] Kharazmi, E., Zhang, Z., Karniadakis, G.E., 2021.  hp-vpinns: Variational physics-informed neural networks with domain decomposition.  Computer Methods in Applied Mechanics and Engineering 374, 113547.  doi:10.1016/j.cma.2020.113547.
  * Kollmannsberger et al. [2021] Kollmannsberger, S., D’Angella, D., Jokeit, M., Herrmann, L., et al., 2021\.  Deep Learning in Computational Mechanics.  Springer.
  * Kong et al. [2018] Kong, Q., Trugman, D.T., Ross, Z.E., Bianco, M.J., Meade, B.J., Gerstoft, P., 2018\.  Machine Learning in Seismology: Turning Data into Insights.  Seismological Research Letters 90, 3–14.  doi:10.1785/0220180259.
  * Kutyniok [2022] Kutyniok, G., 2022.  The mathematics of artificial intelligence.  arXiv preprint arXiv:2203.08890 doi:10.48550/arXiv.2203.08890.
  * Lagaris et al. [1998] Lagaris, I.E., Likas, A., Fotiadis, D.I., 1998.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE transactions on neural networks 9, 987–1000.  doi:10.1109/72.712178.
  * Lagaris et al. [2000] Lagaris, I.E., Likas, A.C., Papageorgiou, D.G., 2000.  Neural-network methods for boundary value problems with irregular boundaries.  IEEE Transactions on Neural Networks 11, 1041–1049.  doi:10.1109/72.870037.
  * Mao et al. [2020] Mao, Z., Jagtap, A.D., Karniadakis, G.E., 2020.  Physics-informed neural networks for high-speed flows.  Computer Methods in Applied Mechanics and Engineering 360, 112789.  doi:10.1016/j.cma.2019.112789.
  * Marone [1998] Marone, C., 1998.  Laboratory-derived friction laws and their application to seismic faulting.  Ann. Rev. Earth Pl. Sc. 26, 643–696.  doi:10.1146/annurev.earth.26.1.643.
  * Marone and Kilgore [1993] Marone, C., Kilgore, B., 1993\.  Scaling of the critical slip distance for seismic faulting with shear strain in fault zones.  Nature 362, 618–621.  doi:10.1038/362618a0.
  * Mishra and Molinaro [2022a] Mishra, S., Molinaro, R., 2022a.  Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for pdes.  IMA Journal of Numerical Analysis 42, 981–1022.  doi:10.1093/imanum/drab032.
  * Mishra and Molinaro [2022b] Mishra, S., Molinaro, R., 2022b.  Estimates on the generalization error of physics-informed neural networks for approximating pdes.  IMA Journal of Numerical Analysis doi:10.1093/imanum/drab093.
  * Müller and Zeinhofer [2023] Müller, J., Zeinhofer, M., 2023\.  Achieving high accuracy with pinns via energy natural gradients.  arXiv preprint arXiv:2302.13163 doi:10.48550/arXiv.2302.13163.
  * National Academies of Sciences, Engineering, and Medicine and others [2020] National Academies of Sciences, Engineering, and Medicine and others, 2020.  A vision for NSF Earth sciences 2020-2030: Earth in time.  National Academies Press, Washington, DC.
  * Okazaki et al. [2022] Okazaki, T., Ito, T., Hirahara, K., Ueda, N., 2022\.  Physics-informed deep learning approach for modeling crustal deformation.  Nature Communications 13, 7092\.  doi:10.1038/s41467-022-34922-1.
  * Raissi et al. [2019] Raissi, M., Perdikaris, P., Karniadakis, G.E., 2019.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational physics 378, 686–707.  doi:10.1016/j.jcp.2018.10.045.
  * Roache [1998] Roache, P., 1998.  Verification and Validation in Computational Science and Engineering.  1 ed., Hermosa Publishers, Albuquerque, NM.
  * Ruina [1983] Ruina, A., 1983.  Slip instability and state variable friction laws.  J. Geophys. Res. Solid Earth 88, 10359–10370.  doi:10.1029/JB088iB12p10359.
  * Scholz [2019] Scholz, C.H., 2019.  The Mechanics of Earthquakes and Faulting.  3 ed., Cambridge University Press.  doi:10.1017/9781316681473.
  * Shin and Em Karniadakis [2020] Shin, YeonjongDarbon, J., Em Karniadakis, G., 2020.  On the convergence of physics informed neural networks for linear second-order elliptic and parabolic type pdes.  Communications in Computational Physics 28, 2042–2074.  doi:10.4208/cicp.OA-2020-0193.
  * Shin et al. [2023] Shin, Y., Zhang, Z., Karniadakis, G.E., 2023.  Error estimates of residual minimization using neural networks for linear pdes.  Journal of Machine Learning for Modeling and Computing 4, 73–101.  doi:10.1615/JMachLearnModelComput.2023050411.
  * Sun et al. [2020] Sun, L., Gao, H., Pan, S., Wang, J.X., 2020.  Surrogate modeling for fluid flows based on physics-constrained deep learning without simulation data.  Computer Methods in Applied Mechanics and Engineering 361, 112732.  doi:10.1016/j.cma.2019.112732.
  * van den Ende et al. [2018] van den Ende, M., Chen, J., Ampuero, J.P., Niemeijer, A., 2018\.  A comparison between rate-and-state friction and microphysical models, based on numerical simulations of fault slip.  Tectonophysics 733, 273–295.  doi:10.1016/j.tecto.2017.11.040.
  * Wang and Cai [2020] Wang, BoZhang, W., Cai, W., 2020\.  Multi-scale deep neural network (mscalednn) methods for oscillatory stokes flows in complex domains.  Communications in Computational Physics 28, 2139–2157.  doi:10.4208/cicp.OA-2020-0192.
  * Yu et al. [2018] Yu, B., et al., 2018.  The deep ritz method: a deep learning-based numerical algorithm for solving variational problems.  Communications in Mathematics and Statistics 6, 1–12.  doi:10.1007/s40304-018-0127-z.
  * Zhao et al. [2023] Zhao, J., Ling, H., Liu, J., Wang, J., Burke, A.F., Lian, Y., 2023.  Machine learning for predicting battery capacity for electric vehicles.  eTransportation 15, 100214\.  doi:10.1016/j.etran.2022.100214.