# Adjoint-based inversion for stress and frictional parameters in earthquake
modeling

Vidar Stiernström, Martin Almquist11footnotemark: 1, Eric M. Dunham Department
of Information Technology, Uppsala University, SE-751 05 Uppsala, Sweden,
<EMAIL_ADDRESS><EMAIL_ADDRESS>of Geophysics,
Stanford University, Stanford CA 94305, USA<EMAIL_ADDRESS>
###### Abstract

We present an adjoint-based optimization method to invert for stress and
frictional parameters used in earthquake modeling. The forward problem is
linear elastodynamics with nonlinear rate-and-state frictional faults. The
misfit functional quantifies the difference between simulated and measured
particle displacements or velocities at receiver locations. The misfit may
include windowing or filtering operators. We derive the corresponding adjoint
problem, which is linear elasticity with linearized rate-and-state friction
with time-dependent coefficients derived from the forward solution. The
gradient of the misfit is efficiently computed by convolving forward and
adjoint variables on the fault. The method thus extends the framework of full-
waveform inversion to include frictional faults with rate-and-state friction.
In addition, we present a space-time dual-consistent discretization of a
dynamic rupture problem with a rough fault in antiplane shear, using high-
order accurate summation-by-parts finite differences in combination with
explicit Runge–Kutta time integration. The dual consistency of the
discretization ensures that the discrete adjoint-based gradient is the exact
gradient of the discrete misfit functional as well as a consistent
approximation of the continuous gradient. Our theoretical results are
corroborated by inversions with synthetic data. We anticipate that adjoint-
based inversion of seismic and/or geodetic data will be a powerful tool for
studying earthquake source processes; it can also be used to interpret
laboratory friction experiments.

##  1 Introduction

Earthquake cycle and dynamic rupture simulations with fault dynamics governed
by rate-and-state friction have emerged as promising tools to better
understand the processes governing earthquake nucleation and seismicity.
Dynamic source models are complementary to kinematic source models that
express the seismic wavefield and solid displacements through a convolution of
the specified space-time history of slip with elastic Green’s functions.
Kinematic models are widely used in slip inversions, which can be set up as
linear least squares problems due to the linear relation between slip and the
wavefield. However, kinematic inversions only reveal how faults slip, but not
why they slip in a certain way. Furthermore, the solutions to the inverse
problem can violate certain physical constraints, such as dissipation rather
than creation of mechanical energy during frictional sliding. Dynamic source
models go beyond kinematic source models through the introduction of friction
laws that are based on laboratory experiments. These friction laws involve
nonlinear relations between fault shear and normal tractions, slip velocity,
and one or more state variables that capture the dependence of frictional
strength on the slip history of the interface. After specification of the
initial stresses and frictional parameters, dynamic rupture and earthquake
cycle models provide both the slip history and wavefield (or quasi-static
solid displacement field) as part of the solution to the nonlinear problem.
Connecting cycle and dynamic rupture simulations to real data permits
determination of the stresses and frictional parameters, which are otherwise
difficult or impossible to constrain. This requires augmenting dynamic source
models with inversion capabilities. The inversion is a PDE-constrained
optimization problem where we seek parameter values that minimize the misfit
between model output and data. We propose a gradient-based optimization where
the gradient of the misfit is computed with only two simulations: one of the
forward problem and one of the adjoint problem. Gradient-based optimization
with the adjoint method is the preferred approach in full waveform inversions
(FWI) in seismic tomography Tromp et al. (2005); Fichtner (2011).

To date, dynamic rupture inversions have been done either by grid-search
algorithms or Bayesian methods that require running thousands to order of a
million forward models (Peyrat and Olsen (2004); Gallovič et al. (2019a, b);
Premus et al. (2022)). Friction inversion for earthquake cycle models has been
posed as a sequential data assimilation problem that is solved, for example,
using the ensemble Kalman filter van Dinther et al. (2019); Hirahara and
Nishikiori (2019). Frictional properties are treated as part of the state
vector and updated in time. However, these parameters should not be time-
dependent and we seek an alternative formulation that respects this.

In this work, we present an adjoint-based optimization framework for inverting
rate-and-state frictional parameters and initial stresses from geodetic and
seismic observations. The method presented in this work extends FWI to include
nonlinear rate-and-state friction laws. Similar to FWI, the adjoint equations
are the same linear elasticity equations as the forward equations, but with
fault dynamics governed by what resembles linearized rate-and-state friction.
The adjoint friction law and state evolution equation include time-dependent
coefficients that are functions of slip velocity and state from the forward
problem. The framework is presented for dynamic elasticity but may be
straightforwardly adapted to quasi-static or quasi-dynamic settings, to be
used in earthquake cycle or aseismic slip simulations. We note that a special
case of this method has been derived for quasi-dynamic cycle models, for the
spatially discretized problem with a boundary element discretization of the
linear elastic slip-stress change relation, by Kano et al. Kano et al. (2013,
2015, 2020).

To demonstrate the adjoint-based optimization framework we consider a 2D
dynamic rupture problem with a rough fault in antiplane shear, discretized
using finite difference methods and explicit Runge–Kutta time stepping
methods. High-order finite difference methods satisfying a summation-by-parts
(SBP) property have successfully been applied to linear elasticity with rate-
and-state friction in e. g. , Kozdon et al. (2013); Erickson and Dunham
(2014); O’Reilly et al. (2015); Duru and Dunham (2016); Duru et al. (2019);
Harvey et al. (2023), and extended to poro- and viscoelasticity in
Torberntsson et al. (2018); Allison and Dunham (2018); Heimisson et al.
(2019). In this work, we utilize the recently developed boundary-optimized
second derivative operators of Stiernström et al. (2023), providing increased
accuracy for problems where boundary or interface effects are of particular
importance. The operators are combined with the non-stiff interface treatment
based on characteristic variables presented in Erickson et al. (2022). The
fault interface conditions are imposed weakly, using simultaneous
approximation terms (SAT). The resulting spatial SBP-SAT discretization allows
for efficient time integration by means of explicit methods. We prove that
this SBP-SAT discretization is dual consistent Hicken and Zingg (2011); Berg
and Nordström (2012); Hicken and Zingg (2014), meaning that the spatial
discretization of the misfit functional and the forward problem results in a
stable and consistent discretization of the adjoint problem. Furthermore, the
gradient of the semi-discrete problem is the exact gradient (to machine
precision) of the semi-discrete misfit functional as well as a consistent
approximation of the continuous gradient. We extend dual consistency to the
fully discrete problem by employing the standard fourth-order Runge–Kutta time
stepping and using the associated Runge–Kutta quadrature to discretize the
time integral in the misfit functional Sanz-Serna (2016). While our
discretization is based on SBP finite differences, dual consistency should
straightforwardly extend to other methods satisfying SBP or corresponding
discrete integration-by-parts properties. This includes discontinuous- and
continuous finite element methods as well as finite volume methods.

Our theoretical findings are corroborated by a series of numerical experiments
on model problems using synthetic data. First, the discrete adjoint-based
gradient is compared to a brute-force finite-difference approximation of the
gradient. Second, we perform inversions for frictional parameters in an
inverse crime setting, as well as using high-resolution synthetic data,
demonstrating the capabilities of the method.

The article is developed as follows: In Section 2 we establish the notation
used in the continuous and discrete analysis. In Section 3 the equations of
linear elasticity and rate-and-state friction are presented. In Section 4 the
inverse problem is formulated as a PDE-constrained optimization problem. The
adjoint-based gradient of the misfit functional is derived in Section 5. We
demonstrate how filtered and windowed residuals in particle displacement and
velocity enter the adjoint equations and discuss well-posedness of the adjoint
problem. In Section 6 we present the dual-consistent SBP-SAT discretization
for dynamic rupture simulations in antiplane shear and use it to invert
synthetic data, demonstrating that the method indeed is capable of
reconstructing frictional parameters. Finally, the article is concluded in
Section 7.

##  2 Notation

###  2.1 Integrals, inner products and index notation

The $L^{2}$-inner product for scalar functions $u$, $v$ on $\Omega$ is denoted

$\left(u,v\right)_{\Omega}:=\int_{\Omega}uv\mathrm{d}\Omega.$ (2.1)

Similarly, define the bilinear form for integrating over a surface
$\Gamma\in\partial\Omega$ as

$\left(u,v\right)_{\Gamma}:=\int_{\Gamma}uv\mathrm{d}\Gamma.$ (2.2)

The same notation is used when integrating over a space-time domain, i. e. ,
for the domain $\Omega\times\mathcal{T}$ with $\mathcal{T}=[0,T]$, we write

$\left(u,v\right)_{\Omega\times\mathcal{T}}:=\int_{0}^{T}\left(u,v\right)_{\Omega}\mathrm{d}t,$
(2.3)

and correspondingly for $\left(u,v\right)_{\Gamma\times\mathcal{T}}$.
Throughout the article, we will make use of index notation. Here the Einstein
summation convention is applied to sum over repeated indices $I,J,K,L$ used to
denote spatial components. For instance, (2.3) is extended to vector-valued
functions $\bar{u}$ and $\bar{v}$ as

$\left(u_{I},v_{I}\right)_{\Omega\times
T}:=\sum_{I=1}^{d}\left(u_{I},v_{I}\right)_{\Omega\times\mathcal{T}},$ (2.4)

where $d$ is the dimension of $\Omega$. Furthermore, we will use the notation
$\partial_{I}:=\frac{\partial}{\partial x_{I}}$, and the summation convention
then applies to derivatives. For instance, the variable coefficient Laplace
operator is $\partial_{I}\mu\partial_{I}$.

###  2.2 Grid functions and discrete operators

Boldface symbols will be used to denote discrete quantities, e. g. , the
discrete grid function approximating the function $u(\bar{x})$,
$\bar{x}\in\Omega$ is denoted $\mathbf{u}$ and is defined on the grid
$\mathbf{\Omega}$. Vector-valued functions and grid functions are denoted with
a bar, e. g. , the grid function restriction of $\bar{u}$ is denoted
$\bar{\mathbf{u}}$. Scalar functions acting as diagonal operators are denoted
with double bars, e. g. , the operator approximating $f(\bar{x})$ is
$\bar{\bar{\mathbf{f}}}:=\text{diag}(\mathbf{f})$.

$\mathbf{H}$ with appropriate subscripts will be used to denote quadrature
rules. For instance, a quadrature on $\mathbf{\Omega}$ is denoted
$\mathbf{H}_{\Omega}$. The notation for inner products is analogous to the
continuous setting. For instance

$\left(\mathbf{u},\mathbf{v}\right)_{\mathbf{\Omega}}:=\mathbf{u}^{T}\mathbf{H}_{\Omega}\mathbf{v},$
(2.5)

For numerical integration along a boundary grid segment $\bm{\Gamma}$, we
first introduce the boundary restriction operators $\mathbf{e}_{\Gamma}$,
where $\mathbf{e}_{\Gamma}^{T}\mathbf{u}\approx u(\bar{x}),\bar{x}\in\Gamma$.
A discrete bilinear form approximating (2.2) is then defined as

$\left(\mathbf{u},\mathbf{v}\right)_{\bm{\Gamma}}:=(\mathbf{e}_{\Gamma}^{T}\mathbf{u})^{T}\mathbf{H}_{\Gamma}(\mathbf{e}_{\Gamma}^{T}\mathbf{v}),$
(2.6)

where $\mathbf{H}_{\Gamma}$ is a quadrature rule on $\bm{\Gamma}$. A bilinear
form over the entire boundary grid $\bm{\partial\Omega}$ is defined by summing
over all boundary grid segments, i. e. ,

$\left(\mathbf{u},\mathbf{v}\right)_{\bm{\partial\Omega}}:=\sum_{\bm{\Gamma}\subset\bm{\partial\Omega}}\left(\mathbf{u},\mathbf{v}\right)_{\bm{\Gamma}}.$
(2.7)

In a semi-discrete setting, where time is left continuous, we will also use
semi-discrete inner products, e. g. ,
$\left(\mathbf{u},\mathbf{v}\right)_{\mathbf{\Omega}\times\mathcal{T}}:=\int_{0}^{T}\mathbf{u}^{T}\mathbf{H}_{\mathbf{\Omega}}\mathbf{v}\mathrm{d}t$.
The summation convention for $I,J,K,L$ also applies to discrete relations. For
instance, for vector-valued grid functions $\bar{\mathbf{v}}$ and
$\bar{\mathbf{u}}$ the extension of (2.5) is

$\left(\mathbf{u}_{I},\mathbf{v}_{I}\right)_{\mathbf{\Omega}}:=\sum_{I=1}^{d}\left(\mathbf{u}_{I},\mathbf{v}_{I}\right)_{\mathbf{\Omega}}.$
(2.8)

##  3 Linear elasticity with rate-and-state friction

Consider a linear elastic medium $\Omega=\Omega_{-}\cup\Omega_{+}$, separated
by a frictional fault $\Gamma=\Omega_{-}\cap\Omega_{+}$. Let $\bar{u}$ denote
the displacement vector. From Hooke’s law, the changes in the stress tensor
$\sigma$ and traction vector $\bar{T}$ for elastic deformations about a
prestressed reference configuration are

$\sigma_{IJ}=C_{IJKL}\partial_{K}u_{L},$ (3.1)

and

$T_{J}=n_{I}\sigma_{IJ},$ (3.2)

where $\bar{n}$ denotes the outward unit normal and $C_{IJKL}$ is the elastic
stiffness tensor. (We note that these stress changes are sometimes denoted as
$\Delta\sigma_{IJ}$.) The corresponding shear traction is the projection of
$\bar{T}$ onto the fault plane, given by

$\bar{\tau}=\bar{T}-(\bar{n}\cdot\bar{T})\bar{n}.$ (3.3)

Henceforth, superscripts $+$ and $-$ will be used to distinguish between
fields on the $+$ and $-$ sides of the fault. The total stress on the fault is
the sum of prestress and changes in stress due to the displacement $\bar{u}$.
The compressive normal stress used for evaluating frictional strength is

$\sigma_{n}^{tot}=\max(0,\sigma_{n}^{0}-\bar{n}^{-}\cdot\bar{T}^{-}),$ (3.4)

which is set to zero in tension, a commonly used procedure that avoids having
to handle inequality conditions associated with fault opening Harris et al.
(2009, 2018); Erickson et al. (2023). The shear traction vector is

$\bar{\tau}^{tot}=\bar{\tau}^{0}+\bar{\tau}^{-}-\eta\bar{V},$ (3.5)

where $\eta\bar{V}$ is the radiation damping term, with slip velocity
$\bar{V}$ defined as the jump in particle velocity across the fault, i. e. ,

$\bar{V}=[\\![\dot{\bar{u}}]\\!]=\dot{\bar{u}}^{+}-\dot{\bar{u}}^{-}.$ (3.6)

The radiation damping term is only used in quasi-dynamic models, where the
coefficient $\eta$ is the impedance of shear waves radiating away from the
fault Rice (1993). In the fully dynamic setting, $\eta$ is set to zero.
Although we only perform experiments with $\eta=0$ in this paper, we include
the radiation damping term in the derivation for completeness.

Force balance on the fault and the usual condition of no penetration or
separation are stated as

$\displaystyle\bar{T}^{+}$ $\displaystyle=-\bar{T}^{-}$ (3.7)
$\displaystyle\bar{n}^{+}\cdot\bar{u}^{+}$
$\displaystyle=-\bar{n}^{-}\cdot\bar{u}^{-}.$ (3.8)

Equating shear traction with frictional shear strength yields

$\bar{\tau}^{tot}=\sigma_{n}^{tot}f\frac{\bar{V}}{\lvert\bar{V}\rvert},$ (3.9)

where $f$ is the friction coefficient. Solving (3.5) and (3.9) for
$\bar{\tau}^{\pm}$ leads to

$\bar{\tau}^{\pm}=\mp\bar{F}(\bar{V},\Psi),$ (3.10)

where the rate-and-state friction law is given by

$\bar{F}(\bar{V},\Psi)=\sigma_{n}^{tot}f(\lvert\bar{V}\rvert,\Psi)\frac{\bar{V}}{\lvert\bar{V}\rvert}-\bar{\tau}^{0}+\eta\bar{V}.$
(3.11)

Here, $f(\lvert\bar{V}\rvert,\Psi)$ is a rate-and-state friction coefficient,
dependent on the dimensionless state variable $\Psi$, a measure of the
interface contact strength related to the past history of sliding. It is
governed by a state evolution equation

$\dot{\Psi}=G(\bar{V},\Psi).$ (3.12)

Let $\bar{Q}(\bar{x},t)$ denote external forcings. The governing equations of
the forward problem are given by

$\begin{array}[]{lll}\vskip 3.0pt plus 1.0pt minus
1.0pt\rho\ddot{u}_{J}=\partial_{I}C_{IJKL}\partial_{K}u_{L}+Q_{J},&\bar{x}\in\Omega,&t\in\mathcal{T},\\\
\vskip 3.0pt plus 1.0pt minus
1.0pt\bar{u}=\bar{u}_{0},\quad\dot{\bar{u}}=\bar{v}_{0},&\bar{x}\in\Omega,&t=0,\\\
\vskip 3.0pt plus 1.0pt minus
1.0ptL\bar{u}=\bar{g},&\bar{x}\in\partial\Omega,&t\in\mathcal{T},\\\ \vskip
3.0pt plus 1.0pt minus
1.0pt\displaystyle\bar{\tau}^{\pm}=\mp\bar{F}(\bar{V},\Psi),&\bar{x}\in\Gamma,&t\in\mathcal{T},\\\
\vskip 3.0pt plus 1.0pt minus
1.0pt\bar{n}^{+}\cdot\bar{T}^{+}=\bar{n}^{-}\cdot\bar{T}^{-},&\bar{x}\in\Gamma,&t\in\mathcal{T},\\\
\vskip 3.0pt plus 1.0pt minus
1.0pt\bar{n}^{+}\cdot\bar{u}^{+}=-\bar{n}^{-}\cdot\bar{u}^{-},&\bar{x}\in\Gamma,&t\in\mathcal{T},\\\
\vskip 3.0pt plus 1.0pt minus
1.0pt\dot{\Psi}=G(\bar{V},\Psi),&\bar{x}\in\Gamma,&t\in\mathcal{T},\\\ \vskip
3.0pt plus 1.0pt minus 1.0pt\Psi=\Psi_{0},&\bar{x}\in\Gamma,&t=0,\end{array}$
(3.13)

where the first equation is the momentum balance equation, $\bar{u}_{0}$ and
$\bar{v}_{0}$ are initial data, and $L\bar{u}=\bar{g}(\bar{x},t)$ denotes
boundary conditions on exterior boundaries, for a boundary operator $L$. The
data $\bar{u}_{0}$, $\bar{v}_{0}$, and $\bar{g}$ are assumed to be independent
of the parameter $p$ that we later invert for. We will consider $L$ such that
for homogeneous boundary data $\bar{g}=\bar{0}$, the boundary conditions are
of the form

$\tau_{J}=-\alpha_{JL}\dot{u}_{L},$ (3.14)

where $\alpha_{JL}=\alpha_{LJ}$ is positive semi-definite. Note that this
includes traction conditions ($\alpha_{JL}=0$), displacement conditions
($\alpha_{JL}\rightarrow\infty$), and characteristic outflow conditions
Petersson and Sjögreen (2009).

##  4 The inverse problem

Now consider $N_{rec}$ receivers positioned at $\bar{x}=\bar{x}_{r}^{(k)}$,
$k=1,\dots,N_{rec}$, each with a time series of measurements
$\bar{m}_{data}^{(k)}(t)$ of either particle displacements ($\bar{m}=\bar{u}$)
or velocities ($\bar{m}=\dot{u}$). We seek parameter values that minimize the
residual $\bar{r}$, defined as the difference between model predictions at the
receiver locations (obtained by solving (3.13)) and the measured data, i. e. ,

$\bar{r}^{(k)}(t)=\bar{m}(\bar{x}_{r}^{(k)},t)-\bar{m}_{data}^{(k)}(t).$ (4.1)

In practice, the residual is often filtered or windowed. We therefore consider
applying linear operators $W^{(k)}$ to obtain the adjusted residuals

$\bar{R}^{(k)}(t)=W^{(k)}\left[\bar{r}^{(k)}\right](t).$ (4.2)

The notation $W^{(k)}\left[\bar{r}^{(k)}\right]$ means that the linear
operator $W^{(k)}$ acts on $\bar{r}^{(k)}$. The result is the time-dependent
adjusted residual. We define the misfit functional $\mathcal{F}$ as the sum of
the $L^{2}$-norms of the adjusted residuals:

$\mathcal{F}=\frac{1}{2}\sum_{k=1}^{N_{rec}}\int\limits_{0}^{T}\lvert\bar{R}^{(k)}(t)\rvert^{2}\,\mathrm{d}t.$
(4.3)

Weighting of the different terms can be included in the $W^{(k)}$ operator.

For a model parameter $p$ in (3.13), the inverse problem is the PDE-
constrained optimization problem given by

$\min_{p}\mathcal{F}\text{ subject to \eqref{eq:elastic_3d_fwd}}.$ (4.4)

To solve (4.4) using gradient-based optimization algorithms, the misfit
gradient $\frac{\delta\mathcal{F}}{\delta p}$ is required. However, taking the
functional derivative of (4.3) directly requires computing
$\frac{\delta\bar{m}}{\delta p}$ (see Section 5.2). In a discrete setting,
where $\mathbf{p}$ is represented by $N$ degrees of freedom, even the simplest
first-order difference approximation of $\frac{\delta\mathcal{F}}{\delta p}$
would require (3.13) to be solved numerically $2N$ times. This is of course
not feasible for large-scale 3D computations. We therefore seek an alternative
expression through the adjoint-state framework.

##  5 Adjoint equations and misfit gradient

We begin this section by presenting the adjoint equations and the adjoint-
based gradient of the misfit (4.3) in the form of a theorem. To this end,
introduce the following adjoint variables:

$\begin{array}[]{ll}\vskip 3.0pt plus 1.0pt minus
1.0pt\bar{u}^{\dagger}&\text{(adjoint displacement)},\\\ \vskip 3.0pt plus
1.0pt minus 1.0pt\Psi^{\dagger}&\text{(adjoint state variable)},\\\ \vskip
3.0pt plus 1.0pt minus
1.0pt\bar{V}^{\dagger}:=-[\\![\dot{\bar{u}}^{\dagger}]\\!]&\text{(adjoint slip
velocity)},\\\ \vskip 3.0pt plus 1.0pt minus
1.0ptT^{\dagger}_{J}:=n_{I}C_{IJKL}\partial_{K}u^{\dagger}_{L}&\text{(adjoint
traction)},\\\ \vskip 3.0pt plus 1.0pt minus
1.0pt\bar{\tau}^{\dagger}=\bar{T}^{\dagger}-(\bar{n}\cdot\bar{T}^{\dagger})\bar{n}&\text{(adjoint
shear traction)}.\end{array}$ (5.1)

The negative sign in $\bar{V}^{\dagger}$ is for notational convenience. As
shortly explained, the adjoint problem will later be phrased in reversed time
in which case $\bar{V}^{\dagger}$ is defined analogously to (3.6).
Furthermore, let $\bar{Q}^{\dagger}(\bar{x},t)$ be the adjoint source term. It
consists of the residuals (4.2) acting as point forces at the receiver
locations. For a complete description, see Section 5.2. For now, it is
sufficient that $\bar{Q}^{\dagger}$ satisfies the relation

$\left(Q_{J}^{\dagger},\frac{\delta\dot{u}_{J}}{\delta
p}\right)_{\Omega\times\mathcal{T}}=\frac{\delta\mathcal{F}}{\delta p}.$ (5.2)

We now define the adjoint equations to the inverse problem (4.4) as

$\begin{array}[]{lll}\vskip 3.0pt plus 1.0pt minus
1.0pt\rho\ddot{u}^{\dagger}_{J}=\partial_{I}C_{IJKL}\partial_{K}u^{\dagger}_{L}+Q^{\dagger}_{J},&\bar{x}\in\Omega,&t\in\mathcal{T},\\\
\vskip 3.0pt plus 1.0pt minus
1.0pt\bar{u}^{\dagger}=\bar{0},\quad\dot{\bar{u}}^{\dagger}=\bar{0},&\bar{x}\in\Omega,&t=T,\\\
\vskip 3.0pt plus 1.0pt minus
1.0ptL^{\dagger}\bar{u}^{\dagger}=0,&\bar{x}\in\partial\Omega,&t\in\mathcal{T},\\\
\vskip 3.0pt plus 1.0pt minus
1.0pt\displaystyle\bar{\tau}^{\dagger\pm}=\mp\bar{F}^{\dagger}(\bar{V}^{\dagger},\Psi^{\dagger}),&\bar{x}\in\Gamma,&t\in\mathcal{T},\\\
\vskip 3.0pt plus 1.0pt minus
1.0pt\bar{n}^{+}\cdot\bar{T}^{\dagger+}=\bar{n}^{-}\cdot\bar{T}^{\dagger-},&\bar{x}\in\Gamma,&t\in\mathcal{T},\\\
\vskip 3.0pt plus 1.0pt minus
1.0pt\bar{n}^{+}\cdot\bar{u}^{\dagger+}=-\bar{n}^{-}\cdot\bar{u}^{\dagger-},&\bar{x}\in\Gamma,&t\in\mathcal{T},\\\
\vskip 3.0pt plus 1.0pt minus
1.0pt-\dot{\Psi}^{\dagger}=G^{\dagger}(\bar{V}^{\dagger},\Psi^{\dagger}),&\bar{x}\in\Gamma,&t\in\mathcal{T},\\\
\vskip 3.0pt plus 1.0pt minus
1.0pt\Psi^{\dagger}=0,&\bar{x}\in\Gamma,&t=T.\end{array}$ (5.3)

The adjoint problem is equipped with terminal conditions (imposed at time
$t=T$) and is naturally solved in reversed time. It is common to introduce the
change of variables $t^{\dagger}=T-t$. We utilize this when discretizing (5.3)
(see Section 6.3). Here, for notational convenience, we keep the original time
$t$ in the derivation of the adjoint-based gradient. Note that
$\frac{\partial}{\partial t}=-\frac{\partial}{\partial t^{\dagger}}$, and
therefore

$-\dot{\Psi}^{\dagger}=\frac{\partial\Psi^{\dagger}}{\partial
t^{\dagger}},\quad\bar{V}^{\dagger}=-[\\![\dot{\bar{u}}^{\dagger}]\\!]=\frac{\partial}{\partial
t^{\dagger}}{[\\![\bar{u}^{\dagger}]\\!]}.$ (5.4)

Notably, the boundary operator is self-adjoint in reversed time, i. e. ,
$L^{\dagger}=L$, and the adjoint boundary conditions can be formulated as

$\tau^{\dagger}_{J}=\alpha_{JL}\dot{u}^{\dagger}_{L}=-\alpha_{JL}\frac{\partial
u^{\dagger}_{L}}{\partial t^{\dagger}},$ (5.5)

(cf. (3.14)). As mentioned in Section 3, (5.5) includes homogeneous
displacement, traction, and characteristic outflow conditions. The relations
(5.4) - (5.5), show that the adjoint problem (5.3), when considered in
reversed time, is of the same form as the forward problem (3.13).

The adjoint friction law and state evolution equation in (5.3) are given by

$F^{\dagger}_{J}(\bar{V}^{\dagger},\Psi^{\dagger})=\frac{\partial
F_{I}}{\partial V_{J}}V^{\dagger}_{I}+\frac{\partial G}{\partial
V_{J}}\Psi^{\dagger},$ (5.6)

and

$G^{\dagger}(\bar{V}^{\dagger},\Psi^{\dagger})=\frac{\partial
F_{J}}{\partial\Psi}V^{\dagger}_{J}+\frac{\partial
G}{\partial\Psi}\Psi^{\dagger}.$ (5.7)

Note that (5.6) - (5.7) resemble the equations of linearized rate-and-state
friction. However, the coefficients are time-dependent and are actually
functions of the forward variables, e. g. $\frac{\partial F_{I}}{\partial
V_{J}}=\frac{\partial F_{I}}{\partial V_{J}}(\bar{V},\Psi)$.

With the adjoint problem defined, we are ready to state the first major result
of this paper.

######  Theorem 5.1.

Let $V^{\dagger*}$, $\Psi^{\dagger}$, satisfy (5.3). Further, let $p$ be a
parameter in $\bar{F}$ or $G$, and let $\Psi_{0}$ be the initial state. Then,
the gradient of the misfit in the inverse problem (4.4) is given by

$\displaystyle\frac{\delta\mathcal{F}}{\delta p}$
$\displaystyle=-\left(V^{\dagger}_{J},\frac{\partial F_{J}}{\partial
p}\right)_{\mathcal{T}}-\left(\Psi^{\dagger},\frac{\partial G}{\partial
p}\right)_{\mathcal{T}},$ (5.8)
$\displaystyle\frac{\delta\mathcal{F}}{\delta\Psi_{0}}$
$\displaystyle=\Psi_{0}^{\dagger},$

where $\Psi^{\dagger}_{0}=\Psi^{\dagger}(t=0)$.

###### Proof.

See section 5.1. ∎

The parameter $p$ in Theorem 5.1 can be any parameter that appears in the
rate-and-state friction law $\bar{F}$ and/or state evolution $G$, including
rate-state parameters and initial stresses.

###  5.1 Derivation of the adjoint-based gradient

To prove Theorem 5.1, we formulate the Lagrangian cost functional

$\displaystyle\mathcal{L}=\mathcal{F}$
$\displaystyle+\left(\dot{u}^{\dagger}_{J},\rho\ddot{u}_{J}-\partial_{I}C_{IJKL}\partial_{K}u_{L}-Q_{J}\right)_{\Omega\times\mathcal{T}}$
(5.9)
$\displaystyle+\left(\Psi^{\dagger},\dot{\Psi}-G\right)_{\Gamma\times\mathcal{T}}.$

Note that for any $\bar{u}$ and $\Psi$ satisfying (3.13),
$\mathcal{L}=\mathcal{F}$ and therefore $\frac{\delta\mathcal{L}}{\delta
p}=\frac{\delta\mathcal{F}}{\delta p}$. Moreover, since initial- and boundary
data in (3.13) are independent of $p$ we may set $\bar{u}_{0}=0$,
$\bar{v}_{0}=0$ and $\bar{g}=0$ for the remaining part of the derivation.

Before considering the first variation $\delta\mathcal{L}$, we rewrite the
cost functional by performing partial integrations to shift derivatives from
$\bar{u}$ to $\bar{u}^{\dagger}$. Integrating by parts in time yields

$\left(\dot{u}^{\dagger}_{J},\rho\ddot{u}_{J}\right)_{\Omega\times\mathcal{T}}=-\left(\rho\ddot{u}^{\dagger}_{J},\dot{u}_{J}\right)_{\Omega\times\mathcal{T}}.$
(5.10)

where terms at $t=0$ and $t=T$ vanish due to the initial and terminal
conditions $\dot{\bar{u}}=0$ and $\dot{\bar{u}}^{\dagger}=0$ in (3.13) and
(5.3). Next, consider the spatial elastic operator and integrate by parts in
space to obtain

$\displaystyle-\left(\dot{u}^{\dagger}_{J},\partial_{I}C_{IJKL}\partial_{K}u_{L}\right)_{\Omega\times\mathcal{T}}=$
$\displaystyle-\left(\dot{u}^{\dagger}_{J},\tau_{J}\right)_{\partial\Omega\times\mathcal{T}}+\left(\dot{\tau}^{\dagger}_{J},u_{J}\right)_{\partial\Omega\times\mathcal{T}}$
(5.11)
$\displaystyle-\left(\dot{u}^{\dagger+}_{J},\tau^{+}_{J}\right)_{\Gamma\times\mathcal{T}}+\left(\dot{\tau}^{\dagger+}_{J},u^{+}_{J}\right)_{\Gamma\times\mathcal{T}}$
$\displaystyle-\left(\dot{u}^{\dagger-}_{J},\tau^{-}_{J}\right)_{\Gamma\times\mathcal{T}}+\left(\dot{\tau}^{\dagger-}_{J},u^{-}_{J}\right)_{\Gamma\times\mathcal{T}}$
$\displaystyle-\left(\partial_{I}C_{IJKL}\partial_{K}\dot{u}^{\dagger}_{L},u_{J}\right)_{\Omega\times\mathcal{T}},$

where we used the major symmetry of the stiffness tensor: $C_{IJKL}=C_{KLIJ}$.
Integrating by parts in time leads to

$\displaystyle\left(\dot{\tau}^{\dagger}_{J},u_{J}\right)_{\partial\Omega\times\mathcal{T}}$
$\displaystyle=-\left(\tau^{\dagger}_{J},\dot{u}_{J}\right)_{\partial\Omega\times\mathcal{T}}$
(5.12)
$\displaystyle\left(\dot{\tau}^{\dagger\pm}_{J},u^{\pm}_{J}\right)_{\Gamma\times\mathcal{T}}$
$\displaystyle=-\left(\tau^{\dagger\pm}_{J},\dot{u}^{\pm}_{J}\right)_{\Gamma\times\mathcal{T}}$

and

$\displaystyle-\left(\partial_{I}C_{IJKL}\partial_{K}\dot{u}^{\dagger}_{L},u_{J}\right)_{\Omega\times\mathcal{T}}=\left(\partial_{I}C_{IJKL}\partial_{K}u^{\dagger}_{L},\dot{u}_{J}\right)_{\Omega\times\mathcal{T}},$
(5.13)

where, again, terms at the initial and final times vanish due to the initial
and terminal conditions $\bar{u}=0$, $\bar{u}^{\dagger}=0$. We have now
derived

$\displaystyle\left(\dot{u}^{\dagger}_{J},\rho\ddot{u}_{J}-\partial_{I}C_{IJKL}\partial_{K}u_{L}\right)_{\Omega\times\mathcal{T}}=$
$\displaystyle-\left(\rho\ddot{u}^{\dagger}_{J}-\partial_{I}C_{IJKL}\partial_{K}u^{\dagger}_{L},\dot{u}_{J}\right)_{\Omega\times\mathcal{T}}$
(5.14)
$\displaystyle\underbrace{-\left(\dot{u}^{\dagger}_{J},\tau_{J}\right)_{\partial\Omega\times\mathcal{T}}-\left(\tau^{\dagger}_{J},\dot{u}_{J}\right)_{\partial\Omega\times\mathcal{T}}}_{BT_{ext}}$
$\displaystyle\underbrace{-\left(\dot{u}^{\dagger+}_{J},\tau^{+}_{J}\right)_{\Gamma\times\mathcal{T}}-\left(\tau^{\dagger+}_{J},\dot{u}^{+}_{J}\right)_{\Gamma\times\mathcal{T}}}_{BT_{fault}^{+}}$
$\displaystyle\underbrace{-\left(\dot{u}^{\dagger-}_{J},\tau^{-}_{J}\right)_{\Gamma\times\mathcal{T}}-\left(\tau^{\dagger-}_{J},\dot{u}^{-}_{J}\right)_{\Gamma\times\mathcal{T}}}_{BT_{fault}^{-}}.$

To simplify this expression, we first note that the exterior boundary terms
$BT_{ext}$ vanish. To see this, we use that the forward and adjoint boundary
conditions satisfy (3.14) and (5.5), and $\alpha_{JL}=\alpha_{LJ}$, such that

$\displaystyle BT_{ext}$
$\displaystyle=\left(\dot{u}^{\dagger}_{J},\alpha_{JL}\dot{u}_{L}\right)_{\partial\Omega\times\mathcal{T}}-\left(\tau^{\dagger}_{J},\dot{u}_{J}\right)_{\partial\Omega\times\mathcal{T}}$
(5.15)
$\displaystyle=\left(\alpha_{JL}\dot{u}^{\dagger}_{L}-\tau^{\dagger}_{J},\dot{u}_{J}\right)_{\partial\Omega\times\mathcal{T}}=0.$

This naturally holds for self-adjoint boundary operators $L=L^{\dagger}$.
Second, for the fault terms $BT_{fault}^{\pm}$, we use the fault conditions
$\bar{\tau}^{+}=-\bar{\tau}^{-}=\bar{F}$,
$\bar{\tau}^{\dagger+}=-\bar{\tau}^{\dagger-}=\bar{F}^{\dagger}$ in (3.13) and
(5.3), together with the definitions of slip velocity
$\bar{V}=[\\![\dot{\bar{u}}]\\!]$,
$\bar{V}^{\dagger}=-[\\![\dot{\bar{u}}^{\dagger}]\\!]$, to obtain

$BT_{fault}^{+}+BT_{fault}^{-}=-\left(V^{\dagger}_{J},F_{J}\right)_{\Gamma\times\mathcal{T}}+\left(V_{J},F^{\dagger}_{J}\right)_{\Gamma\times\mathcal{T}}.$
(5.16)

The last term in the Lagrangian cost functional can be rewritten, using
integration by parts in time,

$\left(\Psi^{\dagger},\dot{\Psi}-G\right)_{\Gamma\times\mathcal{T}}=-\left(\dot{\Psi}^{\dagger},\Psi\right)_{\Gamma\times\mathcal{T}}-\left(\Psi^{\dagger},G\right)_{\Gamma\times\mathcal{T}}+\left(\Psi^{\dagger}_{0},\Psi_{0}\right)_{\Gamma},$
(5.17)

where $\Psi^{\dagger}_{0},\Psi_{0}$ are the state variables at $t=0$. The
corresponding term at $t=T$ vanishes due to the terminal condition
$\Psi(t=T)=0$ in (5.3). We have derived

$\displaystyle\mathcal{L}=\mathcal{F}$
$\displaystyle-\left(\rho\ddot{u}^{\dagger}_{J}-\partial_{I}C_{IJKL}\partial_{K}u^{\dagger}_{L},\dot{u}_{J}\right)_{\Omega\times\mathcal{T}}-\left(\dot{u}^{\dagger}_{J},Q_{J}\right)_{\Omega\times\mathcal{T}}$
(5.18)
$\displaystyle-\left(V^{\dagger}_{J},F_{J}\right)_{\Gamma\times\mathcal{T}}+\left(V_{J},F^{\dagger}_{J}\right)_{\Gamma\times\mathcal{T}}$
$\displaystyle-\left(\dot{\Psi}^{\dagger},\Psi\right)_{\Gamma\times\mathcal{T}}-\left(\Psi^{\dagger},G\right)_{\Gamma\times\mathcal{T}}+\left(\Psi^{\dagger}_{0},\Psi_{0}\right)_{\Gamma}.$

We are now ready to derive the first variation of $\mathcal{L}$. First,
$\delta Q_{J}=0$ since external forces do not depend on frictional parameters
or initial stresses. Second, by the chain rule,

$\delta F_{J}=\frac{\partial F_{J}}{\partial V_{I}}\delta V_{I}+\frac{\partial
F_{J}}{\partial\Psi}\delta\Psi+\frac{\partial F_{J}}{\partial p}\delta p,$
(5.19) $\delta G=\frac{\partial G}{\partial V_{I}}\delta V_{I}+\frac{\partial
G}{\partial\Psi}\delta\Psi+\frac{\partial G}{\partial p}\delta p.$ (5.20)

The above considerations together with the relation (5.2) lead to

$\displaystyle\delta\mathcal{L}=$
$\displaystyle-\left(\rho\ddot{u}^{\dagger}_{J}-\partial_{I}C_{IJKL}\partial_{K}u^{\dagger}_{L}-Q_{J}^{\dagger},\delta\dot{u}_{J}\right)_{\Omega\times\mathcal{T}}$
(5.21) $\displaystyle-\left(V^{\dagger}_{J},\frac{\partial F_{J}}{\partial
V_{I}}\delta V_{I}+\frac{\partial
F_{J}}{\partial\Psi}\delta\Psi+\frac{\partial F_{J}}{\partial p}\delta
p\right)_{\Gamma\times\mathcal{T}}+\left(F^{\dagger}_{J},\delta
V_{J}\right)_{\Gamma\times\mathcal{T}}$
$\displaystyle-\left(\dot{\Psi}^{\dagger},\delta\Psi\right)_{\Gamma\times\mathcal{T}}-\left(\Psi^{\dagger},\frac{\partial
G}{\partial V_{I}}\delta V_{I}+\frac{\partial
G}{\partial\Psi}\delta\Psi+\frac{\partial G}{\partial p}\delta
p\right)_{\Gamma\times\mathcal{T}}$
$\displaystyle+\left(\Psi^{\dagger}_{0},\delta\Psi_{0}\right)_{\Gamma}.$

The volume term vanishes due to the adjoint momentum balance in (5.3).
Gathering the remaining terms leads to

$\displaystyle\delta\mathcal{L}$
$\displaystyle=\left(F^{\dagger}_{J}-\frac{\partial F_{I}}{\partial
V_{J}}V^{\dagger}_{I}-\frac{\partial G}{\partial V_{J}}\Psi^{\dagger},\delta
V_{J}\right)_{\Gamma\times\mathcal{T}}$ (5.22)
$\displaystyle-\left(\dot{\Psi}^{\dagger}+\frac{\partial
F_{J}}{\partial\Psi}V^{\dagger}_{J}+\frac{\partial
G}{\partial\Psi}\Psi^{\dagger},\delta\Psi\right)_{\Gamma\times\mathcal{T}}$
$\displaystyle-\left(\Psi^{\dagger},\frac{\partial G}{\partial p}\delta
p\right)_{\Gamma\times\mathcal{T}}-\left(V^{\dagger}_{J},\frac{\partial
F_{J}}{\partial p}\delta p\right)_{\Gamma\times\mathcal{T}}$
$\displaystyle+\left(\Psi^{\dagger}_{0},\delta\Psi_{0}\right)_{\Gamma}.$

The first two terms vanish due to the adjoint friction law (5.6) and state
evolution equation (5.7). We have arrived at

$\delta\mathcal{L}=-\left(\Psi^{\dagger},\frac{\partial G}{\partial p}\delta
p\right)_{\Gamma\times\mathcal{T}}-\left(V^{\dagger}_{J},\frac{\partial
F_{J}}{\partial p}\delta
p\right)_{\Gamma\times\mathcal{T}}+\left(\Psi^{\dagger}_{0},\delta\Psi_{0}\right)_{\Gamma}.$
(5.23)

Since $\frac{\delta\mathcal{F}}{\delta p}=\frac{\delta\mathcal{L}}{\delta p}$
and $\frac{\delta}{\delta p}\left(v,w\delta
p\right)_{\Gamma\times\mathcal{T}}=\left(v,w\right)_{\mathcal{T}}$, (5.8)
follows. This proves Theorem 5.1.

######  Remark.

Note that in the Lagrangian cost functional (5.9) we have taken the inner
product with adjoint particle velocity $\dot{\bar{u}}^{\dagger}$, rather than
particle displacement $\bar{u}^{\dagger}$. This choice significantly
simplifies the adjoint friction and state evolution equations. If one instead
were to use adjoint displacement, the equations would, in addition to the
current terms, depend also on forward particle accelerations, adjoint slip and
adjoint state rate.

###  5.2 Misfits and adjoint sources

This section details how misfits of particle displacement and velocity enter
the adjoint problem (5.3), i. e. , how the different types of misfits
translate to the adjoint source term $\bar{Q}^{\dagger}$.

Consider the adjusted residuals $\bar{R}^{(k)}$ in (4.2). We will refer to the
linear operators $W^{(k)}$ as filters, although they could equally well be
windowing operators. It is understood that the $W^{(k)}$ may be different for
each component of $\bar{r}^{(k)}$, but to simplify the presentation we avoid
explicit notation for this. Furthermore, we first consider a single receiver
and will temporarily drop the superscript. Since $W$ is linear, the first
variation of $R_{J}$ can be expressed as

$\displaystyle\delta R_{J}$ $\displaystyle=W[\delta r_{J}]=W[\delta
m_{J}]\hat{\delta}_{k},$ (5.24)

where $\hat{\delta}_{k}(\bar{x}):=\hat{\delta}(\bar{x}-\bar{x}_{r}^{(k)})$ is
the shifted Dirac delta function. Let $W^{\dagger}$ denote the Hilbert adjoint
of $W$ with respect to the $L^{2}$-inner product, such that

$\int\limits_{0}^{T}q(t)W[s](t)\,\mathrm{d}t=\int\limits_{0}^{T}W^{\dagger}[q](t)s(t)\,\mathrm{d}t\quad\forall
q,s\in L^{2}(\mathcal{T}).$ (5.25)

Note that windowing operators of the form $W[r](t)=w(t)r(t)$, for some
function $w$, are self-adjoint ($W^{\dagger}=W$), while a general filter might
not be self-adjoint.

Now consider the first variation $\delta F$ of the misfit functional (4.3)
with a single receiver. Carrying out the chain rule for the first variation
yields

$\displaystyle\delta\mathcal{F}$
$\displaystyle=\int\limits_{0}^{T}\int\limits_{\Omega}R_{J}\delta
R_{J}\,\mathrm{d}\bar{x}\,\mathrm{d}t=\int\limits_{0}^{T}\int\limits_{\Omega}R_{J}W[\delta
m_{J}]\hat{\delta}_{k}\,\mathrm{d}\bar{x}\,\mathrm{d}t$ (5.26)
$\displaystyle=\int\limits_{0}^{T}\int\limits_{\Omega}W^{\dagger}[R_{J}]\hat{\delta}_{k}\delta
m_{J}\,\mathrm{d}\bar{x}\,\mathrm{d}t=\left(W^{\dagger}[R_{J}]\hat{\delta}_{k},\delta
m_{J}\right)_{\Omega\times\mathcal{T}}.$

For a velocity misfit where $\bar{m}=\dot{\bar{u}}$, we note that
$Q_{J}^{\dagger}=W^{\dagger}[R_{J}]\hat{\delta}_{k}$ satisfies (5.2), which
was used to derive the adjoint-based gradient. In the case of a displacement
misfit, i. e. , $\bar{m}=\bar{u}$, introduce the time-integrated quantity

$\hat{R}_{J}(t)=\int\limits_{0}^{t}W^{\dagger}[R_{J}](t^{\prime})\,\mathrm{d}t^{\prime}+\hat{R}_{J,0},$
(5.27)

where the constant $\hat{R}_{J,0}$ is selected such that $\hat{R}_{J}(T)=0$.
By (5.27), $\hat{R}_{J}$ satisfies the ODE

$\dot{\hat{R}}_{J}=W^{\dagger}[R_{J}].$ (5.28)

Using integration by parts in time, it therefore follows that

$\delta\mathcal{F}=\left(W^{\dagger}[R_{J}]\hat{\delta}_{k},\delta
u_{J}\right)_{\Omega\times\mathcal{T}}=-\left(\hat{R}_{J}\hat{\delta}_{k},\delta\dot{u}_{J}\right)_{\Omega\times\mathcal{T}},$
(5.29)

where the term at $t=0$ vanishes because of the initial condition $u_{J}=0$
(or in general the $p$-independecy of the inital data) and the term at $t=T$
vanishes because $\hat{R}_{J}(T)=0$. Considering multiple receivers again, the
adjoint source function at receiver $k$ is thus

$S_{J}^{(k)}(t)=\begin{cases}-\hat{R}_{J}^{(k)}(t)&m_{J}=u_{J}\quad\text{(displacement)},\\\
W^{(k)\dagger}[R^{(k)}_{J}](t)&m_{J}=\dot{u}_{J}\quad\text{(velocity)}.\end{cases}$
(5.30)

Note that both types of signals involve the application of an adjoint filter
operator. The complete adjoint source term is obtained by summing over all
receivers,

$Q_{J}^{\dagger}(\bar{x},t)=\sum_{k=1}^{N_{rec}}S_{J}^{(k)}(t)\delta(\bar{x}-\bar{x}_{r}^{(k)}).$
(5.31)

In this work, the misfit functionals are considered to measure either
displacement or velocity exclusively. It is straightforward to extend the
derivation to a linear combination of displacement and velocity misfits. We
briefly comment on extending the misfit to other types of measurements.
Measurements to consider include strain rate from fiber optics cables,
pressure from hydrophones, or snap shot measurements of the displacement of
Earth’s surface through InSAR. If the resulting misfit can be cast in the form
of (5.2), the expression for the adjoint-based gradient in Theorem 5.1 follows
directly. The difference will only be in the adjoint source term. In the case
of an InSAR measurement, we anticipate that the adjoint source term consists
of a surface forcing in space, multiplied by a Heaviside function in time
activating at the snapshot time.

###  5.3 Well-posedness

We now proceed to analyze and discuss the well-posedness of the adjoint
problem (5.3). Since the adjoint equations are similar to the forward
equation, the conditions for well-posedness are also similar. In Duru et al.
(2019) the right-hand-side function $G$ in the state evolution equation (3.12)
is analyzed in terms of well-posedness of the forward problem. State evolution
equations resulting in a well-posed forward problem are termed admissible. The
forward state evolution equation is admissible if $\frac{\partial
G}{\partial\Psi}$ is bounded. For the adjoint-state evolution equation (5.7),
differentiating $G^{\dagger}$ yields

$\frac{\partial G^{\dagger}}{\partial\Psi^{\dagger}}=\frac{\partial
G}{\partial\Psi}.$ (5.32)

Hence, $\frac{\partial G^{\dagger}}{\partial\Psi^{\dagger}}$ is bounded
precisely if the forward state evolution equation is admissible. For the
forward problem, the rate of change of mechanical energy is (see, Duru et al.
(2019))

$\displaystyle\frac{\mathrm{d}E}{\mathrm{d}t}$
$\displaystyle=-\left(V_{J},\tau^{0}_{J}+\tau^{-}_{J}\right)_{\Gamma}=-\left(V_{J},\tau^{0}_{J}+F_{J}\right)_{\Gamma}$
(5.33)
$\displaystyle=-\left(V_{J},\sigma_{n}^{tot}f(\lvert\bar{V}\rvert,\Psi)\frac{V_{J}}{\lvert\bar{V}\rvert}+\eta
V_{J}\right)_{\Gamma}\leq 0.$

Since there is no prestress in the adjoint problem, we can infer that the
adjoint energy rate is

$\frac{\mathrm{d}E^{\dagger}}{\mathrm{d}t^{\dagger}}=-\left(V^{\dagger}_{J},F_{J}^{\dagger}\right)_{\Gamma}=-\left(V_{J}^{\dagger},\frac{\partial
F_{I}}{\partial
V_{J}}V^{\dagger}_{I}\right)_{\Gamma}-\left(V_{J}^{\dagger},\frac{\partial
G}{\partial V_{J}}\Psi^{\dagger}\right)_{\Gamma}.$ (5.34)

For the adjoint problem to be well-posed, the energy rate must be bounded from
above. An important step in bounding the energy rate is to prove that
$\frac{\partial F_{I}}{\partial V_{J}}$ is positive definite, i. e. ,

$V^{\dagger}_{J}\frac{\partial F_{I}}{\partial
V_{J}}V^{\dagger}_{I}>0\quad\forall\bar{V}^{\dagger}.$ (5.35)

To prove (5.35), we differentiate $F_{I}$ to obtain

$\displaystyle\frac{\partial F_{I}}{\partial V_{J}}$
$\displaystyle=\sigma_{n}^{tot}\left(f_{V}(\lvert\bar{V}\rvert,\Psi)\frac{V_{J}}{\lvert\bar{V}\rvert}\frac{V_{I}}{\lvert\bar{V}\rvert}+f(\lvert\bar{V}\rvert,\Psi)\frac{\delta_{IJ}\lvert\bar{V}\rvert^{2}-V_{I}V_{J}}{\lvert\bar{V}\rvert^{3}}\right)+\eta\delta_{IJ}$
(5.36)
$\displaystyle=\frac{\sigma_{n}^{tot}}{\lvert\bar{V}\rvert}\left(f\delta_{IJ}+(f_{V}\lvert\bar{V}\rvert-f)\frac{V_{I}V_{J}}{\lvert\bar{V}\rvert^{2}}\right)+\eta\delta_{IJ}.$

Next, note that
$V_{I}V^{\dagger}_{I}=\gamma\lvert\bar{V}\rvert\lvert\bar{V}^{\dagger}\rvert$,
where $-1\leq\gamma\leq 1$. Therefore, we obtain

$\displaystyle V_{J}^{\dagger}\frac{\partial F_{I}}{\partial
V_{J}}V^{\dagger}_{I}$
$\displaystyle=\frac{\sigma_{n}^{tot}}{\lvert\bar{V}\rvert}\left(f\lvert\bar{V}^{\dagger}\rvert^{2}+(f_{V}\lvert\bar{V}\rvert-f)\frac{\left(V_{I}V^{\dagger}_{I}\right)^{2}}{\lvert\bar{V}\rvert^{2}}\right)+\eta\lvert
V^{\dagger}\rvert^{2}$ (5.37)
$\displaystyle=\frac{\sigma_{n}^{tot}}{\lvert\bar{V}\rvert}\left(f+(f_{V}\lvert\bar{V}\rvert-f)\gamma^{2}\right)\lvert
V^{\dagger}\rvert^{2}+\eta\lvert V^{\dagger}\rvert^{2}$
$\displaystyle\geq\sigma_{n}^{tot}\min\left(\frac{f}{\lvert\bar{V}\rvert},f_{V}\right)\lvert\bar{V}^{\dagger}\rvert^{2}.$

The direct effect is assumed to be rate-strengthening, i.e., $f_{V}>0$, and
the friction coefficient is strictly positive, i. e. , $f>0$. Both
$\frac{f}{\lvert\bar{V}\rvert}$ and $f_{V}$ may tend to zero as
$\lvert\bar{V}\rvert\rightarrow\infty$, but since the forward problem is well-
posed $\lvert\bar{V}\rvert<\infty$. Therefore, (5.37) shows that
$\frac{\partial F_{I}}{\partial V_{J}}$ is positive definite.

Next, defining $G_{J}^{\prime}=\frac{\partial G}{\partial V_{J}}$, then
$\lvert\bar{G^{\prime}}\rvert<\infty$, again due to the well-posedness of the
forward problem. From (5.37) it follows that

$\displaystyle F_{J}^{\dagger}V^{\dagger}_{J}$
$\displaystyle=V_{J}^{\dagger}\frac{\partial F_{I}}{\partial
V_{J}}V^{\dagger}_{I}+V_{J}^{\dagger}\frac{\partial G}{\partial
V_{J}}\Psi^{\dagger}$ (5.38)
$\displaystyle\geq\sigma_{n}^{tot}\min\left(\frac{f}{\lvert\bar{V}\rvert},f_{V}\right)\lvert\bar{V}^{\dagger}\rvert^{2}-\lvert\bar{G^{\prime}}\rvert\lvert\bar{V}^{\dagger}\rvert\lvert\Psi^{\dagger}\rvert.$

Completing the squares leads to

$F_{J}^{\dagger}V^{\dagger}_{J}\geq-\frac{\lvert\bar{G^{\prime}}\rvert^{2}}{4\sigma_{n}^{tot}\min\left(\frac{f}{\lvert\bar{V}\rvert},f_{V}\right)}(\Psi^{\dagger})^{2}.$
(5.39)

Defining

$\beta=\frac{\lvert\bar{G^{\prime}}\rvert^{2}}{4\sigma_{n}^{tot}\min\left(\frac{f}{\lvert\bar{V}\rvert},f_{V}\right)}\geq
0,$ (5.40)

we may write down the energy estimate

$\frac{\mathrm{d}E^{\dagger}}{\mathrm{d}t^{\dagger}}\leq\left(\Psi^{\dagger},\beta\Psi^{\dagger}\right)_{\Gamma}.$
(5.41)

It remains to bound $\Psi^{\dagger}$, which proves to be more difficult.
Returning to the adjoint state evolution equation (5.7), it is clear that if
$\frac{\partial G}{\partial\Psi}<0$, then the term $\frac{\partial
G}{\partial\Psi}\Psi^{\dagger}$ in $G^{\dagger}$ features exponential decay in
$\Psi^{\dagger}$. This is true for typical state evolution equations, e. g. ,
the slip and aging law Duru et al. (2019). It is therefore very likely that
$\Psi^{\dagger}$ indeed is bounded. However, for a strict bound we must
further bound $\frac{\partial F_{J}}{\partial\Psi}V^{\dagger}_{J}$, which has
so far eluded us. Although very likely to be the case, we conclude that well-
posedness of the adjoint problem (5.3) remains an open problem. The numerical
studies presented herein are performed using the typical rate-and-state
friction coefficient (6.2) together with the slip law (6.4). No blow-ups of
the adjoint solution have been observed, and the results therefore indicate
that the adjoint equations are well-posed, at least in this setting.

##  6 Dynamic rupture in antiplane shear

We now proceed to demonstrate the method applied to dynamic rupture
simulations in antiplane shear. Consider the domain $\Omega$ with a frictional
fault along $\Gamma$, as illustrated in Figure 6.1. In the semi-discrete
setting, the discretizations will be carried out blockwise in $\Omega_{\pm}$.
For this reason, quantities on the respective domain will be denoted with the
subscripts $+$ and $-$, or compactly when the equations apply to both domains
individually by $\pm$. Letting $(x_{1},x_{2})=(x,y)\in\Omega$ be the in-plane
coordinates, the governing equations (3.13) for particle displacement $u$ in
the antiplane direction $x_{3}=z$ reduce to

$\begin{array}[]{lll}\rho_{\pm}\ddot{u}_{\pm}=\partial_{I}\mu_{\pm}\partial_{I}u_{\pm},&\bar{x}\in\Omega_{\pm},&t\in\mathcal{T},\\\
u_{\pm}=u_{0\pm},\quad\dot{u}_{\pm}=v_{0\pm},&\bar{x}\in\Omega_{\pm},&t=0,\\\
\tau_{\pm}+Z\dot{u}_{\pm}=0,&\bar{x}\in\partial\Omega_{\pm}\setminus\Gamma,&t\in\mathcal{T},\\\
\tau_{\pm}=\mp F(V,\Psi),&\bar{x}\in\Gamma,&t\in\mathcal{T},\\\
\dot{\Psi}=G(V,\Psi),&\bar{x}\in\Gamma,&t\in\mathcal{T},\\\
\Psi=\Psi_{0},&\bar{x}\in\Gamma,&t=0.\end{array}$ (6.1)

On the exterior boundaries $\partial\Omega_{\pm}\setminus\Gamma$
characteristic outflow conditions Engquist and Majda (1977) are imposed. On
the fault, we consider a typical rate-and-state friction coefficient Rice et
al. (2001)

$f(|V|,\Psi)=a\sinh^{-1}\left(\frac{|V|}{2V_{0}}e^{\frac{\Psi}{a}}\right),$
(6.2)

where $a$ is the dimensionless direct effect parameter and $V_{0}$ an
arbitrarily chosen reference velocity. In addition, an external loading
$\tau_{L}(\bar{x})$, $\bar{x}\in\Gamma$ is added to the fault, such that the
friction law reads

$F(V,\Psi)=\sigma_{n}f(|V|,\Psi)\frac{V}{|V|}-\tau_{0}-\tau_{L}.$ (6.3)

The external loading could of course be included directly in $\tau_{0}$, but
is kept separate here since we later will perform inversions for $\tau_{0}$
with a known $\tau_{L}$. State evolution is governed by the slip law Marone
(1998); Lapusta et al. (2000)

$G(V,\Psi)=-\frac{|V|}{D_{c}}\left(f(|V|,\Psi)-f_{ss}(|V|)\right),$ (6.4)

where $f_{ss}(|V|)=f_{0}-(b-a)\ln{\frac{|V|}{V_{0}}}$ is the steady state
friction coefficient Rice et al. (2001). Here $b$ is the state evolution
parameter, $D_{c}$ the state evolution distance and $f_{0}$ is the reference
coefficient for steady sliding at velocity $V_{0}$.

In the following sections, we introduce a dual consistent Hicken and Zingg
(2014) spacetime discretization based on SBP difference operators and
Runge–Kutta time integration.

Figure 6.1: Domain with a fault separating $\Omega^{+}$ and $\Omega^{-}$

###  6.1 SBP finite difference discretization

Discretizations of (6.1) using SBP finite differences, with weak enforcement
of fault and boundary conditions through the simultaneous-approximation-term
(SAT) method have been presented in e. g. Erickson and Dunham (2014); Duru et
al. (2019); Erickson et al. (2022); Harvey et al. (2023). Here, we utilize the
boundary-optimized SBP difference operators presented in Stiernström et al.
(2023) to discretize the scalar wave operator on the physical grid
$\mathbf{\Omega}$. The physical grid is obtained through transfinite
interpolation of the coordinate mapping between $\Omega$ and a two-block
Cartesian reference domain. See e. g. Almquist and Dunham (2020, 2021);
Stiernström et al. (2023) for further details on the finite difference
discretization of the scalar wave operator and the discretization of
$\mathbf{\Omega}$.

Let $\mathbf{D}_{II}(\bm{\mu}_{\pm})\approx\partial_{I}\mu_{\pm}\partial_{I}$
be SBP finite difference approximation of the variable coefficient Laplace
operator on the physical domains $\mathbf{\Omega}_{\pm}$. Then, the
discretization of the scalar wave equation in (6.1) is given by

$\begin{array}[]{ll}\bar{\bar{\bm{\rho}}}_{\pm}\ddot{\mathbf{u}}_{\pm}=\mathbf{D}_{II}(\bm{\mu}_{\pm})\mathbf{u}_{\pm}+\text{SAT}_{\pm},&t\in\mathcal{T},\\\
\mathbf{u}_{\pm}=\mathbf{u}_{0\pm},\quad\dot{\mathbf{u}}_{\pm}=\mathbf{v}_{0\pm},&t=0,\\\
\end{array}$ (6.5)

where $\text{SAT}_{\pm}$ weakly enforces boundary and interface conditions.
Imposing the non-linear rate-and-state friction law through a standard
traction SAT (e. g. as in Almquist and Dunham (2020)) may lead to significant
stiffness in dynamic rupture problems. To resolve the issue, we use the non-
stiff SAT presented in Erickson et al. (2022); Harvey et al. (2023), enforcing
boundary and interface conditions through characteristic variables.

Consider first a single block $\Omega$ (either $\Omega_{+}$ or $\Omega_{-}$).
On a boundary segment $\bm{\zeta}$ (either an exterior boundary or the fault
interface), the non-stiff SAT is given by

$\text{SAT}(\mathbf{u},\mathbf{u}^{*},\bm{\tau}^{*})=\mathbf{H}^{-1}_{\Omega}\left(\mathbf{e}\mathbf{H}_{\zeta}(\bm{\tau}^{*}-\bm{\tau})-\mathbf{T}^{T}\mathbf{H}_{\zeta}(\mathbf{u}^{*}-\mathbf{u})\right),$
(6.6)

where $\mathbf{u}$, $\mathbf{u}^{*}$, and $\bm{\tau}^{*}$ all are evaluated on
$\bm{\zeta}$. Here, $\mathbf{H}_{\Omega}$ and $\mathbf{H}_{\Gamma}$ are the
SBP quadratures on $\mathbf{\Omega}$ and $\bm{\zeta}$, while $\mathbf{e}$ is
the boundary restriction operator. Moreover, $\mathbf{T}$ is the boundary
traction operator such that $\bm{\tau}:=\mathbf{T}\mathbf{u}\approx\mu
n_{I}\partial_{I}u(\bar{x})$, $\bar{x}\in\zeta$, while $\mathbf{u}^{*}$ and
$\bm{\tau}^{*}$ are target values (or fluxes) for displacement and traction.
The target values are defined to preserve the outgoing grid-based
characteristic variables,

$\mathbf{w}^{*}:=\bm{Z}\dot{\mathbf{u}}^{*}-\bm{\tau}^{*}=\bm{Z}\dot{\mathbf{u}}-\tilde{\bm{\tau}}:=\mathbf{w},\\\
$ (6.7)

where $\bm{Z}=\sqrt{\bm{\rho}\bm{\mu}}$ is the shear impedance and
$\tilde{\bm{\tau}}$ is the modified traction

$\tilde{\bm{\tau}}=\bm{\tau}+\gamma(\mathbf{u}^{*}-\mathbf{u}).$ (6.8)

The parameter $\gamma$ is the semi-bounded Dirichlet penalty parameter of
Almquist and Dunham (2020). Note that the units of $\gamma$ are such that
$\gamma(\mathbf{u}^{*}-\mathbf{u})$ has the units of traction. Moreover,
$\mathbf{u}^{*}$ and $\bm{\tau}^{*}$ should satisfy the boundary or fault
conditions.

On an exterior boundary, the boundary conditions are conveniently expressed by
introducing a reflection coefficient $\mathbf{R}$, such that we may write the
ingoing characteristic as
$\bm{Z}\dot{\mathbf{u}}^{*}+\bm{\tau}^{*}=\mathbf{R}\mathbf{w}^{*}$. This
together with (6.7) results in

$\displaystyle\bm{\tau}^{*}$
$\displaystyle=\frac{\mathbf{R}-\mathbf{I}}{2}(\bm{Z}\dot{\mathbf{u}}-\tilde{\bm{\tau}}),$
(6.9) $\displaystyle\dot{\mathbf{u}}^{*}$
$\displaystyle=\frac{\mathbf{R}+\mathbf{I}}{2}(\dot{\mathbf{u}}-\bm{Z}^{-1}\tilde{\bm{\tau}}).$
(6.10)

Here $\mathbf{I}$ is an identity matrix with dimensions matching the number of
boundary points. Then, displacement (Dirichlet), traction (Neumann) or
characteristic boundary conditions are imposed by setting

$\begin{array}[]{ll}\mathbf{R}=-\mathbf{I}&\text{(displacement)},\\\
\mathbf{R}=\mathbf{I}&\text{(traction)},\\\
\mathbf{R}=\mathbf{0}&\text{(characteristic outflow)},\\\ \end{array}$ (6.11)

along the corresponding boundary segment. On the outer boundaries,
characteristic outflow conditions are specified, and we therefore set
$\mathbf{R}=0$.

On the fault, the target tractions instead satisfy the fault condition

$\bm{\tau}_{\pm}^{*}=\mp\mathbf{F}(\mathbf{V}^{*},\bm{\Psi}),$ (6.12)

where $\mathbf{F}$ is considered to apply point-wise. From (6.7) and (6.12)
the target slip velocity $\mathbf{V}^{*}$ can be shown to satisfy

$\bm{\kappa}\mathbf{V}^{*}+\mathbf{F}(\mathbf{V}^{*},\bm{\Psi})=-\bm{\tau}_{\ell},$
(6.13)

where $\bm{\kappa}=\frac{\bm{Z}_{+}\bm{Z}_{-}}{\bm{Z}_{+}+\bm{Z}_{-}}$ and
$\bm{\tau}_{\ell}=\frac{\bm{Z}_{+}\mathbf{w}_{-}-\bm{Z}_{-}\mathbf{w}_{+}}{\bm{Z}_{+}+\bm{Z}_{-}}$.
For non-linear $F$, solving (6.13) requires employing a root-finding
algorithm. Note that $\bm{\Psi}$ is held fixed while solving (6.13). Let
$\tilde{\mathbf{V}}^{*}(\bm{\tau}_{\ell})$ denote the solution to

$\mathbf{F}(\tilde{\mathbf{V}}^{*},\bm{\Psi})=-\bm{\tau}_{\ell},$ (6.14)

which can be solved analytically with our friction coefficient (6.2). As noted
in Erickson et al. (2022) one can show that
$\mathbf{V}^{*}\in[-\tilde{\mathbf{V}}^{*}(\bm{\tau}_{\ell}),0]$ for
$\bm{\tau}_{\ell}\geq 0$ and
$\mathbf{V}^{*}\in[0,-\tilde{\mathbf{V}}^{*}(\bm{\tau}_{\ell})]$ for
$\bm{\tau}_{\ell}<0$. Since $F_{V}>0$ for admissible friction laws, there
exists a unique solution $\mathbf{V}^{*}$.

From (6.7) it directly follows that $\dot{\mathbf{u}}^{*}_{\pm}$ on the fault
is given by

$\dot{\mathbf{u}}^{*}_{\pm}=\dot{\mathbf{u}}_{\pm}-\bm{Z}_{\pm}^{-1}(\tilde{\bm{\tau}}_{\pm}-\bm{\tau}_{\pm}^{*}).$
(6.15)

Note that the characteristic boundary and fault conditions require the target
displacements $\mathbf{u}^{*}_{\pm}$ to be regarded as additional unknown grid
functions, evolved together with $\mathbf{u}_{\pm}$ and $\bm{\Psi}$ through
(6.10) and (6.15). For notational purposes we introduce

$\dot{\mathbf{u}}^{*}_{\pm}=\mathbf{L}^{*}(\dot{\mathbf{u}}_{\pm},\tilde{\bm{\tau}}_{\pm})=\begin{cases}\frac{\mathbf{R}+\mathbf{I}}{2}(\dot{\mathbf{u}}_{\pm}-\bm{Z}_{\pm}^{-1}\tilde{\bm{\tau}}_{\pm})&\text{(exterior
boundaries)},\\\
\dot{\mathbf{u}}_{\pm}-\bm{Z}_{\pm}^{-1}(\tilde{\bm{\tau}}_{\pm}-\bm{\tau}_{\pm}^{*})&\text{(fault)},\end{cases}$
(6.16)

with $\bm{\tau}_{\pm}^{*}$ given by (6.12). Similarly, the linear operator for
the right-hand-side of (6.5) is denoted

$\mathbf{L}(\mathbf{u}_{\pm},\mathbf{u}^{*}_{\pm},\bm{\tau}^{*}_{\pm}):=\mathbf{D}_{II}(\bm{\mu}_{\pm})\mathbf{u}_{\pm}+\sum_{\bm{\zeta}}\text{SAT}(\mathbf{u}_{\pm},\mathbf{u}^{*}_{\pm},\bm{\tau}^{*}_{\pm}),$
(6.17)

The semi-discretization of the state evolution equation (6.4) is

$\begin{array}[]{ll}\dot{\bm{\Psi}}=\mathbf{G}(\mathbf{V}^{*},\bm{\Psi}),&t\in\mathcal{T},\\\
\bm{\Psi}=\bm{\Psi}_{0},&t=0,\\\ \end{array}$ (6.18)

where $\mathbf{G}$ is evaluated pointwise. Note that $\mathbf{V}^{*}$ is used
in evolving $\bm{\Psi}$. This will be of importance when deriving the gradient
to the discrete optimization problem in Section 6.2.

Finally, the SBP-SAT semi-discretization of (6.1) is given by

$\begin{array}[]{ll}\bar{\bar{\bm{\rho}}}_{\pm}\ddot{\mathbf{u}}_{\pm}=\mathbf{L}(\mathbf{u}_{\pm},\mathbf{u}^{*}_{\pm},\bm{\tau}^{*}_{\pm}),&t\in\mathcal{T},\\\
\dot{\mathbf{u}}^{*}_{\pm}=\mathbf{L}^{*}(\dot{\mathbf{u}}_{\pm},\tilde{\bm{\tau}}_{\pm})&t\in\mathcal{T},\\\
\mathbf{u}_{\pm}=\mathbf{u}_{0\pm},\quad\dot{\mathbf{u}}_{\pm}=\mathbf{v}_{0\pm},&t=0,\\\
\mathbf{u}^{*}_{\pm}=\mathbf{e}^{T}_{\partial\Omega_{\pm}}\mathbf{u}_{0\pm},&t=0,\\\
\dot{\bm{\Psi}}=\mathbf{G}(\mathbf{V}^{*},\bm{\Psi}),&t\in\mathcal{T},\\\
\bm{\Psi}=\bm{\Psi}_{0},&t=0.\end{array}$ (6.19)

In Duru et al. (2019); Erickson et al. (2022) it is shown (by using the SBP
properties in Appendix A) that the scheme is energy stable and that the semi-
discrete energy rate is a consistent approximation of the continuous energy
rate.

With the forward problem discretized in space, the semi-discretization of the
optimization problem (4.4) reads

$\min_{\mathbf{p}}\bm{\mathcal{F}}\quad\text{subject to
\eqref{eq:anti_plane_shear_discr_fwd}},$ (6.20)

where the semi-discrete misfit is given by

$\bm{\mathcal{F}}=\frac{1}{2}\sum_{k=1}^{N_{rec}}\int_{0}^{T}|\mathbf{r}^{(k)}(t)|^{2}\mathrm{d}t,$
(6.21)

and the semi-discrete residuals are

$\mathbf{r}^{(k)}(t)=\left(\left(\hat{\bm{\delta}}_{k},\mathbf{m}(t)\right)_{\mathbf{\Omega}}-m_{data}(t)\right),$
(6.22)

for a measured field $\mathbf{m}$. We refrain from including filtering
operators $W$ in the semi-discrete presentation since no filtering is used in
the numerical experiments performed in Section 6.5. The discrete Dirac delta
$\hat{\bm{\delta}}_{k}\approx\hat{\delta}(\bar{x}-\bar{x}^{(k)}_{r})$
presented in Petersson et al. (2016) is used to restrict the measured field to
the receiver location $\mathbf{x}_{r}$. The number of moment conditions used
when constructing $\hat{\bm{\delta}}_{k}$ is the same as the order of the SBP
difference operators. No smoothness conditions are used.

###  6.2 Adjoint scheme and semi-discrete gradient

Analogously to the forward problem, the adjoint equations (5.3) reduce to

$\begin{array}[]{lll}\rho_{\pm}\ddot{u}^{\dagger}_{\pm}=\partial_{I}\mu_{\pm}\partial_{I}u^{\dagger}_{\pm}+Q^{\dagger}_{\pm},&\bar{x}\in\Omega_{\pm},&t\in\mathcal{T},\\\
u^{\dagger}_{\pm}=0,\quad\dot{u}^{\dagger}_{\pm}=0,&\bar{x}\in\Omega_{\pm},&t=T,\\\
\tau^{\dagger}_{\pm}-Z\dot{u}^{\dagger}_{\pm}=0,&\bar{x}\in\partial\Omega_{\pm}\setminus\Gamma,&t\in\mathcal{T},\\\
\tau^{\dagger}_{\pm}=\mp
F^{\dagger}(V^{\dagger},\Psi^{\dagger}),&\bar{x}\in\Gamma,&t\in\mathcal{T},\\\
-\dot{\Psi}^{\dagger}=G^{\dagger}(V^{\dagger},\Psi^{\dagger}),&\bar{x}\in\Gamma,&t\in\mathcal{T},\\\
\Psi^{\dagger}=0,&\bar{x}\in\Gamma,&t=T.\end{array}$ (6.23)

In the semi-discrete setting, we start by defining

$\mathbf{w}^{\dagger*}:=-\bm{Z}\dot{\mathbf{u}}^{\dagger*}-\bm{\tau}^{\dagger*}=-\bm{Z}\dot{\mathbf{u}}^{\dagger}-\tilde{\bm{\tau}}^{\dagger}:=\mathbf{w}^{\dagger},$
(6.24)

which are the outgoing characteristic adjoint variables on a boundary segment
after time reversal, where $\tilde{\bm{\tau}}^{\dagger}$ is defined
analogously to (6.8). On an exterior boundary, $\bm{\tau}^{\dagger*}$ and
$\dot{\mathbf{u}}^{\dagger*}$ then satisfies

$\displaystyle\bm{\tau}^{\dagger*}$
$\displaystyle=\frac{\mathbf{R}-\mathbf{I}}{2}(-\bm{Z}\dot{\mathbf{u}}^{\dagger}-\tilde{\bm{\tau}}^{\dagger}),$
(6.25) $\displaystyle-\dot{\mathbf{u}}^{\dagger*}$
$\displaystyle=\frac{\mathbf{R}+\mathbf{I}}{2}(-\dot{\mathbf{u}}^{\dagger}-\bm{Z}^{-1}\tilde{\bm{\tau}}^{\dagger}).$
(6.26)

Similarly, on the fault

$\displaystyle\bm{\tau}_{\pm}^{\dagger*}$
$\displaystyle=\mp\mathbf{F}^{\dagger}(\mathbf{V}^{\dagger*},\bm{\Psi}^{\dagger})=\mp(\bar{\bar{\mathbf{F}}}_{V}\mathbf{V}^{\dagger*}+\bar{\bar{\mathbf{G}}}_{V}\bm{\Psi}^{\dagger}),$
(6.27) $\displaystyle-\dot{\mathbf{u}}^{\dagger*}_{\pm}$
$\displaystyle=-\dot{\mathbf{u}}^{\dagger}_{\pm}-\bm{Z}_{\pm}^{-1}(\tilde{\bm{\tau}}^{\dagger}_{\pm}-\bm{\tau}_{\pm}^{\dagger*}),$
(6.28)

where $\mathbf{F}_{V}(\mathbf{V}^{*},\bm{\Psi})\approx\frac{\partial
F}{\partial V}$ and
$\mathbf{G}_{V}(\mathbf{V}^{*},\bm{\Psi})\approx\frac{\partial G}{\partial V}$
are evaluated point-wise. The adjoint target slip velocity now satisfies

$\bm{\kappa}\mathbf{V}^{\dagger*}+\mathbf{F}^{\dagger}(\mathbf{V}^{\dagger*},\bm{\Psi}^{\dagger})=-\bm{\tau}_{\ell}^{\dagger}=-\frac{\bm{Z}_{+}\mathbf{w}_{-}^{\dagger}-\bm{Z}_{-}\mathbf{w}_{+}^{\dagger}}{\bm{Z}_{+}+\bm{Z}_{-}}.$
(6.29)

Note that $F^{\dagger}$ is linear in $V^{\dagger}$, such that (6.29) reduces
to a time-dependent linear equation in $\mathbf{V}^{\dagger*}$. The semi-
discrete adjoint state evolution equation is simply

$\begin{array}[]{ll}-\dot{\bm{\Psi}}^{\dagger}=\mathbf{G}^{\dagger}(\mathbf{V}^{\dagger*},\bm{\Psi}^{\dagger})=(\bar{\bar{\mathbf{G}}}_{\Psi}\bm{\Psi}^{\dagger}+\bar{\bar{\mathbf{F}}}_{\Psi}\mathbf{V}^{\dagger*}),&t\in\mathcal{T},\\\
\bm{\Psi}^{\dagger}=0,&t=T.\\\ \end{array}$ (6.30)

Again, $\mathbf{G}_{\Psi}(\mathbf{V}^{*},\bm{\Psi})\approx\frac{\partial
G}{\partial\Psi}$ and
$\mathbf{F}_{\Psi}(\mathbf{V}^{*},\bm{\Psi})\approx\frac{\partial
F}{\partial\Psi}$ are evaluated point-wise.

Adopting the notation for the operators $\mathbf{L}$ and $\mathbf{L}^{*}$ in
(6.17) and (6.16), the adjoint equations are then discretized according to
(6.19), resulting in

$\begin{array}[]{ll}\bar{\bar{\bm{\rho}}}_{\pm}\ddot{\mathbf{u}}_{\pm}^{\dagger}=\mathbf{L}(\mathbf{u}^{\dagger}_{\pm},\mathbf{u}^{\dagger*}_{\pm},\bm{\tau}^{\dagger*}_{\pm})+\mathbf{Q}^{\dagger}_{\pm},&t\in\mathcal{T},\\\
-\dot{\mathbf{u}}^{\dagger*}_{\pm}=\mathbf{L}^{*}(-\dot{\mathbf{u}}^{\dagger}_{\pm},\tilde{\bm{\tau}}^{\dagger}_{\pm})&t\in\mathcal{T},\\\
\mathbf{u}^{\dagger}_{\pm}=0,\quad\dot{\mathbf{u}}^{\dagger}_{\pm}=0,&t=T,\\\
\mathbf{u}^{\dagger*}_{\pm}=0,&t=T,\\\
-\dot{\bm{\Psi}}^{\dagger}=\mathbf{G}^{\dagger}(\mathbf{V}^{\dagger*},\bm{\Psi}^{\dagger}),&t\in\mathcal{T},\\\
\bm{\Psi}^{\dagger}=0,&t=T.\end{array}$ (6.31)

The adjoint residual source term is given by
$\mathbf{Q}^{\dagger}_{\pm}=\sum_{k=1}^{N_{rec}}\mathbf{S}^{(k)}(t)\hat{\bm{\delta}}_{k}$
where

$\mathbf{S}^{(k)}(t)=\begin{cases}-\hat{\mathbf{r}}^{(k)}(t),&\mathbf{m}=\mathbf{u},\\\
\mathbf{r}^{(k)}(t),&\mathbf{m}=\dot{\mathbf{u}},\end{cases}$ (6.32)

with
$\hat{\mathbf{r}}^{(k)}(t)=\int_{0}^{t}\mathbf{r}(t^{\prime})dt^{\prime}+\hat{\mathbf{r}}^{(k)}_{0}$,
with $\hat{\mathbf{r}}_{0}^{(k)}$ chosen such that
$\hat{\mathbf{r}}^{(k)}(T)=0$. See Appendix B for a derivation of the residual
source term.

We may now state the gradient to (6.20). It is given by the following Lemma
which is the semi-discrete counterpart to Theorem 5.1 for antiplane shear.

######  Lemma 6.1.

Let $\mathbf{V}^{\dagger*}$, $\bm{\Psi}^{\dagger}$, satisfy (6.31). Further,
let $\mathbf{p}$ be the fault grid function of a parameter in $F$ or $G$, and
let $\bm{\Psi}_{0}$ be the initial state. If $\mathbf{V}^{*}$ is computed
exactly from (6.13), then the gradient of the misfit in the semi-discrete
inverse problem (6.20) is given by

$\displaystyle\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}}$
$\displaystyle=-\int_{0}^{T}\mathbf{H}_{\Gamma_{i}}\left(\bar{\bar{\mathbf{F}}}_{p}\mathbf{V}^{\dagger*}+\bar{\bar{\mathbf{G}}}_{p}\bm{\Psi}^{\dagger}\right)\mathrm{d}t,$
(6.33) $\displaystyle\frac{\partial\bm{\mathcal{F}}}{\partial\bm{\Psi}_{0}}$
$\displaystyle=\mathbf{H}_{\Gamma_{i}}\bm{\Psi}^{\dagger}_{0},$

where $\mathbf{F}_{p}(\mathbf{V}^{*},\bm{\Psi})\approx\frac{\partial
F}{\partial p}$ and
$\mathbf{G}_{p}(\mathbf{V}^{*},\bm{\Psi})\approx\frac{\partial G}{\partial p}$
are evaulated point-wise and
$\bm{\Psi}^{\dagger}_{0}=\bm{\Psi}^{\dagger}(t=0)$.

###### Proof.

See Appendix C. ∎

Note that (6.33) is a consistent approximation of (5.8). In proving Lemma 6.1
we further show that the forward scheme (6.19) is dual consistent Hicken and
Zingg (2011); Berg and Nordström (2012); Hicken and Zingg (2014). The
condition on the accuracy of $\mathbf{V}^{*}$ is discussed in the remark at
the end of Appendix C. In summary, errors in $\mathbf{V}^{*}$ translate to an
error in the gradient.

###  6.3 Time integration and the discrete gradient

To integrate the forward and adjoint schemes (6.19), (6.31) in time, the
standard fourth-order accurate Runge–Kutta method (RK4) will be used. To this
end, introduce $\mathbf{v}_{\pm}=\dot{\mathbf{u}}_{\pm}$ such that (6.19) may
be rewritten to first-order in time as

$\begin{array}[]{ll}\frac{\mathrm{d}{\mathbf{U}}(t)}{\mathrm{d}t}=\mathbf{A}(t,\mathbf{U}),&t\in\mathcal{T},\\\
\mathbf{U}=\mathbf{U}_{0},&t=0,\end{array}$ (6.34)

where

$\begin{array}[]{ll}\mathbf{A}(t,\mathbf{U})=\begin{bmatrix}\mathbf{A}_{-}(t,\mathbf{U}_{-})\\\
\mathbf{A}_{+}(t,\mathbf{U}_{+})\\\
\mathbf{G}(\mathbf{V}^{*},\bm{\Psi})\end{bmatrix},&\mathbf{U}=\begin{bmatrix}\mathbf{U}_{-}\\\
\mathbf{U}_{+}\\\ \bm{\Psi}\end{bmatrix},\end{array}$ (6.35)

and

$\begin{array}[]{ll}\mathbf{A}_{\pm}(t,\mathbf{U}_{\pm})=\begin{bmatrix}\mathbf{v}_{\pm}\\\
\bar{\bar{\bm{\rho}}}^{-1}_{\pm}(\mathbf{L}(\mathbf{u}_{\pm},\mathbf{u}^{*}_{\pm},\bm{\tau}^{*}_{\pm})+\mathbf{Q}_{\pm}(t))\\\
\mathbf{L}^{*}(\mathbf{v}_{\pm},\tilde{\bm{\tau}}_{\pm})\\\
\end{bmatrix},&\mathbf{U}_{\pm}=\begin{bmatrix}\mathbf{u}_{\pm}\\\
\mathbf{v}_{\pm}\\\ \mathbf{u}^{*}_{\pm}\end{bmatrix}.\end{array}$ (6.36)

The ODE system (6.34) is then integrated in time using RK4, where in each step
(6.13) is solved for $\mathbf{V}^{*}$ using bisection.

To integrate the adjoint scheme (6.19) the reversed time $t^{\dagger}=T-t$ is
introduced. Then, the adjoint ODE system is

$\begin{array}[]{ll}\frac{\mathrm{d}{\mathbf{U}}^{\dagger}(t^{\dagger})}{\mathrm{d}t^{\dagger}}=\mathbf{A}^{\dagger}(t^{\dagger},\mathbf{U}^{\dagger}),&t^{\dagger}\in\mathcal{T},\\\
\mathbf{U}^{\dagger}=\mathbf{U}_{0}^{\dagger},&t^{\dagger}=0.\end{array}$
(6.37)

Here

$\begin{array}[]{ll}\mathbf{A}^{\dagger}(t^{\dagger},\mathbf{U}^{\dagger})=\begin{bmatrix}\mathbf{A}^{\dagger}_{-}(t^{\dagger},\mathbf{U}^{\dagger}_{-})\\\
\mathbf{A}^{\dagger}_{+}(t^{\dagger},\mathbf{U}^{\dagger}_{-})\\\
\mathbf{G}^{\dagger}(\mathbf{V}^{\dagger*},\bm{\Psi}^{\dagger})\end{bmatrix},&\mathbf{U}=\begin{bmatrix}\mathbf{U}^{\dagger}_{-}\\\
\mathbf{U}^{\dagger}_{+}\\\ \bm{\Psi}^{\dagger}\end{bmatrix},\end{array}$
(6.38)
$\begin{array}[]{ll}\mathbf{A}^{\dagger}_{\pm}(t^{\dagger},\mathbf{U}^{\dagger}_{\pm})=\begin{bmatrix}\mathbf{v}^{\dagger}_{\pm}\\\
\bar{\bar{\bm{\rho}}}^{-1}_{\pm}(\mathbf{L}(\mathbf{u}^{\dagger}_{\pm},\mathbf{u}^{\dagger*}_{\pm},\bm{\tau}^{\dagger*}_{\pm})+\mathbf{Q}^{\dagger}_{\pm}(t^{\dagger}))\\\
\mathbf{L}^{*}(\mathbf{v}^{\dagger}_{\pm},\tilde{\bm{\tau}}^{\dagger}_{\pm})\\\
\end{bmatrix},&\mathbf{U}^{\dagger}_{\pm}=\begin{bmatrix}\mathbf{u}^{\dagger}_{\pm}\\\
\mathbf{v}^{\dagger}_{\pm}\\\ \mathbf{u}^{\dagger*}_{\pm}\\\
\end{bmatrix},\end{array}$ (6.39)

where for each $t$ (6.29) is solved for $\mathbf{V}^{\dagger*}$. Note that the
forward variables $\mathbf{V}^{*}$ and $\bm{\Psi}$ entering as variable
coefficients in $\mathbf{F}^{\dagger}$ and $\mathbf{G}^{\dagger}$, also are
evaluated in the reversed time.

To evaluate the integrals over time in (6.21) and (6.33), the Runge–Kutta
quadrature is utilized. That is, for time interval $n$ with step size $\Delta
t_{n}$, the Runge–Kutta quadrature weight at substage $s\in[1,4]$ is given by
$\mathbf{H}_{\mathcal{T},4n+s}=\Delta t_{n}b_{s}$, where $b=[1/6,1/3,1/3,1/6]$
Butcher (2016). Then for $N$ time steps, the $4N\times 4N$ diagonal matrix
$\mathbf{H}_{\mathcal{T}}$ is a quadrature on the temporal grid
$\bm{\mathcal{T}}$. The fully discrete optimization problem reads

$\min_{\mathbf{p}}\bm{\mathcal{F}}\quad\text{subject to
\eqref{eq:antiplane_shear_ode} solved using RK4},$ (6.40)

where

$\bm{\mathcal{F}}=\frac{1}{2}\sum_{k=1}^{N_{rec}}\sum_{n=0}^{N-1}\sum_{s=1}^{4}\mathbf{H}_{\mathcal{T},4n+s}|\widehat{\mathbf{r}}^{(k)}_{4n+s}|^{2}.$
(6.41)

with the discrete residual given by

$\widehat{\mathbf{r}}^{(k)}_{4n+s}=\left(\hat{\bm{\delta}}_{k},\widehat{\mathbf{m}}^{(k)}_{4n+s}\right)_{\mathbf{\Omega}}-m^{(k)}_{data}(t_{n}+c_{s}\Delta
t_{n}).$ (6.42)

Here $\widehat{\mathbf{m}}^{(k)}$ is the measured field at receiver $k$ for
all Runge–Kutta stages between $t=0$ and $t=T$, $c=[0,1/2,1/2,1]$ are the RK4
stage nodes, and $t_{0}=0$. Note that for each index $i$,
$\widehat{\mathbf{m}}^{(k)}_{i}$ is a grid function.

We now arrive at one of the major results of this work, namely the discrete
counterpart to Theorem 5.1 for dynamic rupture inversions in antiplane shear.

######  Theorem 6.2.

Let $\widehat{\mathbf{V}}^{\dagger*}$, $\widehat{\bm{\Psi}}^{\dagger}$, be the
spacetime solutions (including Runge–Kutta substages) to (6.37) solved using
RK4. Further, let $\mathbf{p}$ be the fault grid function of a parameter in
$F$ or $G$, and let $\bm{\Psi}_{0}$ be the initial state. If in each time step
$\mathbf{V}^{*}$ is computed exactly from (6.13), then the gradient of the
misfit in the discrete inverse problem (6.40) is given by

$\displaystyle\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}}$
$\displaystyle=-\sum_{n=0}^{N-1}\sum_{s=1}^{4}\mathbf{H}_{\mathcal{T},4n+s}\mathbf{H}_{\Gamma}\left(\bar{\bar{\mathbf{F}}}_{p,4n+s}\widehat{\mathbf{V}}^{\dagger*}_{4n+s}+\bar{\bar{\mathbf{G}}}_{p,4n+s}\widehat{\bm{\Psi}}^{\dagger}_{4n+s}\right),$
(6.43) $\displaystyle\frac{\partial\bm{\mathcal{F}}}{\partial\bm{\Psi}_{0}}$
$\displaystyle=\mathbf{H}_{\Gamma}\widehat{\bm{\Psi}}^{\dagger}_{0},$

where
$\mathbf{F}_{p}(\widehat{\mathbf{V}}^{*},\widehat{\bm{\Psi}})\approx\frac{\partial
F}{\partial p}$,
$\mathbf{G}_{p}(\widehat{\mathbf{V}}^{*},\widehat{\bm{\Psi}})\approx\frac{\partial
G}{\partial p}$ are evaluated point-wise using the Runge–Kutta stage
approximations.

###### Proof.

RK4 is self-adjoint with respect to $\mathbf{H}_{\mathcal{T}}$ under time
reversal, see Sanz-Serna (2016); Matsuda and Miyatake (2021). Therefore, by
Lemma 6.1, the result follows. ∎

Again, errors in the non-linear solve of $\mathbf{V}^{*}$ will result in an
error in $\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}}$, see the remark
at the end of Appendix C.

To summarize, when solving the optimization problem (6.40), the gradient of
(6.41) is obtained through the following procedure:

  1. 1.

Solve (6.34) using RK4, storing the Runge–Kutta stage approximations
$\widehat{\mathbf{V}}^{*}$, $\widehat{\bm{\Psi}}$ and $\widehat{\mathbf{r}}$.

  2. 2.

Solve (6.37) using RK4, with $\widehat{\mathbf{V}}^{*}$, $\widehat{\bm{\Psi}}$
and the residuals $\widehat{\mathbf{r}}^{(k)}$ in $\mathbf{Q}^{\dagger}$
reversed in time. Use the time steps from the forward solve in reverse. Store
the Runge–Kutta stage approximations $\widehat{\mathbf{V}}^{\dagger*}$,
$\widehat{\bm{\Psi}}^{\dagger}$.

  3. 3.

Reverse $\widehat{\mathbf{V}}^{\dagger*}$, $\widehat{\bm{\Psi}}^{\dagger}$ in
time and compute the gradient according to (6.43).

###  6.4 Lower-dimensional parameter representation

The misfit functional in the inverse problem (6.40) is typically non-convex
with many local minima. As a way to regularize (6.40) by reducing the discrete
solution space, a lower-dimensional parameter representation of $\mathbf{p}$
may therefore be beneficial. Examples of such lower-dimensional
representations are polynomials or splines Samareh (2001); Bader et al.
(2023). Here we make use of a coarse grid representation of the parameters and
move between the lower-dimensional parameter grid and the higher-dimensional
computational grid using interpolation operators.

Denote an inversion parameter on the coarse fault grid by $\mathbf{p}^{C}$ and
let $\mathbf{I}^{C2F}$ be an interpolation operator from the coarse grid to
the fine computational grid. Using superscript $F$ for fault grid functions on
the computational grid, $\mathbf{p}^{F}$ is given by

$\mathbf{p}^{F}=\mathbf{I}^{C2F}\mathbf{p}^{C}.$ (6.44)

The gradient expressions in (6.43) are derived with respect to
$\mathbf{p}^{F}$, and are of the form

$\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}^{F}}=\mathbf{H}^{F}_{\Gamma}\bm{\Phi}^{F},$
(6.45)

where $\bm{\Phi}^{F}$ is a combination of forward and adjoint fields. By the
chain rule, the gradient with respect to $\mathbf{p}^{C}$ is

$\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}^{C}}=\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}^{F}}\frac{\partial\mathbf{p}^{F}}{\partial\mathbf{p}^{C}}.$
(6.46)

Differentiating (6.44) stated in indical notation (where $i,j,k$ denote
components in the matrices and vectors) leads to

$\frac{\partial\mathbf{p}^{F}_{i}}{\partial\mathbf{p}^{C}_{k}}=\mathbf{I}_{ij}^{C2F}\delta_{jk}=\mathbf{I}_{ik}^{C2F}=(\mathbf{I}^{C2F})^{T}_{ki}.$
(6.47)

Hence,

$\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}^{C}}=(\mathbf{I}^{C2F})^{T}\mathbf{H}^{F}_{\Gamma}\bm{\Phi}^{F}.$
(6.48)

If $\mathbf{I}^{C2F}$ is any interpolation operator, it is not guaranteed that
$(\mathbf{I}^{C2F})^{T}$ is a consistent approximation of interpolation of
fields from the fine to the coarse grid. Thus, it is not guaranteed that
$\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}^{C}}$ is consistent with
the continuous gradient. However, by requiring that $\mathbf{I}^{C2F}$ is
constructed such that

$\mathbf{I}^{F2C}:=(\mathbf{H}_{\Gamma}^{C})^{-1}(\mathbf{I}^{C2F})^{T}\mathbf{H}^{F}_{\Gamma},$
(6.49)

is a consistent interpolation operator from the fine grid to the coarse, then
(6.48) reads

$\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}^{C}}=\mathbf{H}^{C}_{\Gamma}\mathbf{I}^{F2C}\bm{\Phi}^{F}.$
(6.50)

This corresponds to first interpolating $\bm{\Phi}^{F}$ to the coarse grid and
then applying the coarse fault quadrature. In Almquist et al. (2019)
interpolation operators satisfying (6.49) are shown to exist, provided
$\mathbf{H}_{\Gamma}^{C,F}$ are accurate quadratures. Here, we use
intermediate glue grids to construct interpolation operators for the boundary-
optimized SBP difference operators, see Kozdon and Wilcox (2016); Stiernström
and Almquist (2022).

######  Remark.

Replacing $(\mathbf{I}^{C2F})^{T}$ in (6.48) with any accurate interpolation
operator $\mathbf{I}^{F2C}$ would result in consistency with the continuous
gradient and an approximation of (6.45) to the order of accuracy of
$\mathbf{I}^{F2C}$. The resulting gradient computation would then correspond
to a continuous-adjoint approach Hicken and Zingg (2014). The additional
requirement (6.49) guarantees that the continuous-adjoint approach is
equivalent to the discrete-adjoint approach, meaning that we obtain the exact
gradient of (6.41) while remaining consistent with the continuous gradient.
The property (6.49) states that
$\mathbf{I}^{F2C}=(\mathbf{I}^{C2F})^{\dagger}$ in the $\mathbf{H}_{\Gamma}$
norm. It is termed SBP or inner product preserving and is commonly used to
couple SBP schemes across non-conforming interfaces Mattsson and Carpenter
(2010); Nissen et al. (2015); Lundquist et al. (2018); Kozdon and Wilcox
(2016); Almquist et al. (2019); Lundquist et al. (2020). Interestingly, the
property is useful also in the present setting.

###  6.5 Numerical studies

In this section, we corroborate our theoretical results and present inversions
for friction parameters in dynamic rupture simulations using synthetic data.
The experiments are performed on the domain in Figure 6.1, where the two
elastic blocks are separated by a rough self-similar fractal fault Dunham et
al. (2011); Fang and Dunham (2013). The parameter set used for synthetic data
is listed in Table 6.1. Here $\Gamma_{VW}=\\{\bar{x}\in\Gamma,x\in[-5,6]\text{
km}\\}$ denotes the velocity weakening (VW) region of the fault, i. e. , the
part of the fault where $b-a>0$. In contrast, the region where $b-a<0$ is
velocity strengthening (VS). The initial stress tensor is such that the normal
stress $\sigma_{n}$ is constant while shear traction varies along the fault
according to $\tau_{0}=\sigma_{yz}\hat{n}_{y}$, where $\hat{n}_{y}$ is the
(spatially variable) $y$-component of the fault normal. A one-dimensional
Gaussian distribution centered about $x_{c}=3$ km, with amplitude $A=25$ MPa
and standard deviation $d=2$ km is added as loading to the fault, i. e. ,
setting $\tau_{L}(\bar{x})=Ae^{-\frac{(x-x_{c})^{2}}{2d^{2}}}$ in (6.3). The
problem setup is such that the external loading will cause an earthquake to
nucleate at $x=x_{c}$. Left- and right-going ruptures will propagate along the
fault until they reach the VS region where they then will halt.

Using $m$ grid points in the $x$-direction, the domain is discretized such
that the grid spacing is approximately equal in both directions, resulting in
$m\times(m+1)$ grid points in total across the multiblock grid
$\mathbf{\Omega}=\mathbf{\Omega}_{-}\cup\mathbf{\Omega}_{+}$. The number of
grid points used for lower-dimensional parameter representations is denoted
$m_{p}$. 8th-order accurate boundary-optimized SBP operators are used for the
spatial discretization. In Stiernström et al. (2023) it was shown that these
operators provide an efficient alternative for wave propagation problems where
boundary or interface effects are significant. The non-linear equation (6.13)
for target slip velocity $\mathbf{V}^{*}$ is solved using bisection,
specifying an absolute tolerance $10^{-13}$ m/s. A forward simulation using
$m=1001$ and $\Delta t=0.0005$ s up to $T=6$ s is used to generate high-
resolution synthetic data. Snapshots of the velocity are presented in Figure
6.2, illustrating the nucleation, propagation, and arrest of the rupture. In
Figure 6.3 the slip $[\\![\mathbf{u}]\\!]$ on the fault is presented.

(a) $t=2.00$ s

(b) $t=3.00$ s

(c) $t=4.50$ s

(d) $t=5.00$ s

Figure 6.2: Velocity (m/s) at different points in time, showing rupture propagation and arrest. Figure 6.3: Slip $[\\![\mathbf{u}]\\!]$ (m) plotted every $t=0.25$ s. Table 6.1: Parameter values used in reference problem setup. Parameter | Symbol | Value
---|---|---
Density | $\rho_{\pm}$ | 2.6700 g/$\text{cm}^{3}$
Shear modulus | $\mu_{\pm}$ | 32.0381 GPa
Initial stress | $\sigma_{n}$ | 120 MPa
| $\sigma_{yz}$ | 72 MPa
Reference friction coefficient for steady sliding | $f_{0}$ | 0.6
Direct effect parameter | $a$ | 0.009 $\bar{x}\in\Gamma_{VW}$,
|  | 0.013 otherwise
State evolution effect parameter | $b$ | 0.011
State evolution distance | $D_{c}$ | 0.2 m $\bar{x}\in\Gamma_{VW}$,
|  | 1.0 m otherwise
Reference slip velocity | $V_{0}$ |  $10^{-6}$ m/s
Initial state | $\Psi_{0}$ | 0.7243
Initial displacement | $u_{0\pm}$ | 0 m
Initial velocity | $v_{0\pm}$ |  $\mp 5\cdot 10^{-13}$ m/s

####  6.5.1 Verification of discrete gradient

To verify Theorem 6.2 we compare with a first-order finite difference
approximation of the gradient. Consider the parameter-normalized norm

$\|\mathbf{u}\|_{\mathbf{p}}=\|\bar{\bar{\mathbf{p}}}^{-1}\mathbf{u}\|_{\infty},$
(6.51)

and the relative error

$e(\Delta
p)=\frac{\|\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}}-\mathbf{D}_{+}\bm{\mathcal{F}}(\Delta
p)\|_{\mathbf{p}}}{\|\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}}\|_{\mathbf{p}}},$
(6.52)

where $\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}}$ is the adjoint-
based gradient and $\mathbf{D}_{+}\bm{\mathcal{F}}(\Delta p)$ is the first-
order finite difference approximation, with the $i$th component given by

$\left(\mathbf{D}_{+}\bm{\mathcal{F}}(\Delta
p)\right)_{i}=\frac{\bm{\mathcal{F}}(\mathbf{p}+\Delta
p\mathbf{e}_{i})-\bm{\mathcal{F}}(\mathbf{p})}{\Delta p}.$ (6.53)

Here, $\mathbf{e}_{i}$ is a vector with a 1 at the $i$th entry and zeros
elsewhere. Synthetic data is generated up to $T=6$ s using the parameter
values in Table 6.1, and discretization parameters $m=101$, $m_{p}=11$ and
$\Delta t=0.005$ s. Note that the grid is too coarse for accurate dynamic
rupture simulations but the setup is still useful in verifying that the
correct discrete gradient is obtained since any approximation errors in the
gradient should be clearly visible. To interpolate between the computational
and the parameter grid we construct interpolation operators based on the 8th-
order accurate boundary-optimized norm, and the second-order accurate
traditional SBP norm on an equispaced grid. This means that the interpolation
is second-order accurate. The reason for using a lower-order accurate
interpolation is that it reduces oscillations when interpolating non-smooth
parameters. In total $N_{rec}=88$ receivers are placed in the cut-out
rectangle with outer boundary $-9\leq x\leq 9$, $-9\leq y\leq 9$, and inner
boundary $-7\leq x\leq 7$, $-3\leq y\leq 3$, spaced 2 km apart, see Figure
6.4.

Figure 6.4: Position of $N_{rec}=88$ receivers spaced 2 km apart ($\times$).
Velocity weakening part of fault (- -).

Setting $a_{0}=1.1a$ as the initial value, the error $e(\Delta a)$ is then
computed for $\Delta a$ decreasing from $10^{-5}$ to $10^{-12}$. The results
are presented in Figure 6.5 where it can be observed that $e(\Delta a)$
decreases with a first-order rate with decreasing $\Delta a$ until $e(\Delta
a)\approx 1.8\cdot 10^{-6}$, at $\Delta a=10^{-8.5}$ with a displacement
misfit and $e(\Delta a)\approx 1.4\cdot 10^{-6}$, at $\Delta a=10^{-9}$ with a
velocity misfit. At these points, cancellation errors in the first-order
approximation $\mathbf{D}_{+}\bm{\mathcal{F}}(\Delta a)$ take effect after
which $e(\Delta a)$ increases linearly (although oscillatory) as $\Delta a$ is
further decreased. Note that if
$\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}}$ included an
approximation error $\varepsilon>|\mathbf{D}_{+}\bm{\mathcal{F}}(\Delta a)|$
then the error curve would become horizontal at $e(\Delta a)=\varepsilon$.
Since no such piecewise constant parts of the error curves are observed, we
conclude that (6.43) indeed is the gradient of $\bm{\mathcal{F}}$, and that
errors in computing $\mathbf{V}^{*}$ are negligible.

Figure 6.5: Relative error $e(\Delta a)$ between the adjoint-based gradient
and the first order finite-difference gradient, as a function of $\Delta a$.

####  6.5.2 Inverse crimes

We now proceed to perform inversions in an inverse crime setting, i.e., we aim
to reconstruct parameters based on synthetic data generated from a forward
simulation using the same spatio-temporal resolution. The optimization problem
(6.40) is solved using the quasi-Newton limited-memory Broyden-Fletcher-
Goldfarb-Shanno (L-BFGS) algorithm provided by the Matlab function fmincon.
The forward and adjoint solves are performed using $m=251$ and $\Delta
t=0.003$ s, with $m_{p}$ varying depending on the parameter considered for
inversion. As before, synthetic data is generated using the parameters in
Table 6.1 as the true values.

Inverting for the direct effect parameter $a$ we set $m_{p}=26$, and start
from the initial value $a=0.0135$, which corresponds to a 50% error in the VW
region. Figure 6.6 presents $\mathbf{a}$ plotted for different iterations.
After 200 iterations, $\mathbf{a}$ matches the true value fairly well in the
VW region. Note that the difference in $\mathbf{a}$ in the first and final
iterations is negligible in the VS parts of the fault. This is likely due to
the lack of signal emitted, resulting in a very low sensitivity (small value
of $\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}})$ in the VS region.
See Figure 6.3 where it is clear that the fault has not experienced any
significant slip in the VS region.

(a) 50 iterations

(b) 100 iterations

(c) 150 iterations

(d) 200 iterations

Figure 6.6: Direct effect parameter $\mathbf{a}$ at different iterations for
$m_{p}=26$.

To illustrate the effect of the parameter resolution $m_{p}$ on the
optimization procedure, we perform the same inversion using $m_{p}=51$ and
$m_{p}=101$. The computed values of $\mathbf{a}$ after 200 iterations, and the
misfit histories are presented in Figures 6.7. Increasing the resolution on
the parameter grid, setting $m_{p}=51$, we see larger oscillations about the
hypocenter $x_{c}$. This shows that reducing the discrete parameter space may
be beneficial for faster convergence of the misfit, and for reducing the
number of local minima. The oscillations around the hypocenter observed for
higher $m_{p}$ could likely be remedied by introducing regularization to the
misfit, for instance by penalizing large gradients in $\mathbf{a}$. Finding a
suitable regularization is however out of the scope of this paper.

(a) 200 iterations

(b) Misfit history

Figure 6.7: Direct effect parameter $\mathbf{a}$ and misfit history for
different $m_{p}$

Next, we invert for the shear traction $\tau_{0}=\sigma_{yz}\hat{n}_{y}$,
setting $m_{p}=51$, and starting from $\tau_{0}$ given by $\sigma_{yz}=68$ MPa
corresponding to approximately 5% initial error. The considerably smaller
error in the initial guess is due to inversions of $\tau_{0}$ proving to be
significantly more difficult for this problem setup, with a higher likelihood
of incorrectly converging to a local minimum. Two inverse-crime inversions are
performed, one with the receiver distribution in Figure 6.4 and one where the
number of receivers is increased by placing them 1 km apart and decreasing the
inner bounding box to $-6\leq x\leq 7$, $-2\leq y\leq 2$, totaling
$N_{rec}=325$ (see Figure 6.9). In Figure 6.8 $\bm{\tau}_{0}$ obtained with
$N_{rec}=88$ is plotted for different iterations. The misfit history is shown
in Figure 10(b). The optimization stagnates at a local minimum after about 150
iterations. A comparison of the results at 200 iterations with $N_{rec}=88$
and $N_{rec}=325$ is shown in Figure 6.10. The result using $N_{rec}=325$
matches the true value well in the VW region.

(a) 50 iterations

(b) 100 iterations

(c) 150 iterations

(d) 200 iterations

Figure 6.8: Shear traction $\bm{\tau}_{0}$ at different iterations using
$N_{rec}=88$. Figure 6.9: Position of $N_{rec}=325$ receivers spaced 1 km
apart ($\times$). Velocity weakening part of fault (- -).

(a) $\bm{\tau}_{0}$ at 200 iterations.

(b) Normalized misfit $\bm{\mathcal{F}}/N_{rec}$ history.

Figure 6.10: Comparison of shear traction $\bm{\tau}_{0}$ at 200 iterations
for $N_{rec}=88$ and $N_{rec}=325$.

######  Remark.

Numerical experiments (not presented here) indicate that the optimization
procedure is more robust when approaching the true value from a VS simulation.
For instance, if initially $b-a>0$ in parts of the fault that truly are VS,
the optimization procedure is very likely to converge to an incorrect local
minimum.

####  6.5.3 Inversion using high-resolution synthetic data

Slightly increasing the difficulty of the optimization problem, we perform
inversions using the same setups as in Section 6.5.2, but use high-resolution
data generated from the forward simulation presented at the beginning of this
section (illustrated in Figures 6.2 \- 6.3). In this setting, we effectively
have no modeling errors (other than the perturbed parameters considered for
inversion), but the data contain features that are not resolvable on the
computational grid. When inverting for $a$ we use $m_{p}=26$ and $N_{rec}=88$
receivers, while for $\tau_{0}$, we set $m_{p}=51$ and $N_{rec}=325$. The
results after 200 iterations are presented in Figures 6.11 \- 6.12. In this
setting, $\mathbf{a}$ oscillates around the true value in the VW region, while
$\bm{\tau}_{0}$ seems to interpolate the true value fairly well.

In Figure 6.13 the velocity field is presented for high-resolution data, the
solution using the true parameter set on the computational grid, the computed
parameters, and the initial parameters. Comparing Figure 13(a) to 13(b) it is
clear that spurious oscillations due to lower spatial resolution are present
in the velocity fields on the computational grid. The oscillations observed in
$\mathbf{a}$ in Figure 6.11 are therefore likely a result of numerical errors
in the forward and adjoint fields being mapped into the parameter. Increasing
the number of receivers to $N_{rec}=325$ had no significant effect, and the
results are therefore omitted. To further improve the results, a combination
of regularization and filtering of the misfit residual is likely needed.

(a) $\mathbf{a}$ after 200 iterations

(b) Misfit history

Figure 6.11: Inversion for direct effect parameter $\mathbf{a}$ using high-
resolution data and $N_{rec}=88$.

(a) $\bm{\tau}_{0}$ after 200 iterations

(b) Misfit history

Figure 6.12: Inversion for shear traction $\bm{\tau}_{0}$ using high-
resolution data and $N_{rec}=325$.

(a) True parameters, computational grid

(b) True parameters, high-resolution grid

(c) $\mathbf{a}$ 200 iterations

(d) $\mathbf{a}$ initial parameters

(e) $\bm{\tau}_{0}$ 200 iterations

(f) $\bm{\tau}_{0}$ initial parameters

Figure 6.13: Comparison of velocity (m/s) at $t=5$ using the true parameter
sets on the high-resolution grid and the computational grid, $\mathbf{a}$ and
$\bm{\tau}_{0}$ obtained after 200 L-BFGS iterations and initial parameter
values. Note the difference in scales using the initial parameter values.

##  7 Conclusion

This work presents an adjoint-based optimization framework for earthquake
modeling with rate-and-state friction. The main contributions are two-fold.
Firstly, the continuous adjoint equations to linear elasticity with rate-state
frictional faults are derived, considering displacement and velocity misfit
functionals. The adjoint equations satisfy the equations of linear elasticity
in reversed time, together with an adjoint friction law and state evolution
equation. The adjoint friction law and state evolution equation resemble those
of linearized rate-and-state friction but with time-dependent variable
coefficients from the forward problem. By convolving the adjoint fields with
derivatives of the friction law and state evolution equation, the gradient of
the misfit with respect to the inversion parameter is obtained.

Secondly, we present an SBP-SAT finite difference discretizations of the
forward and adjoint equations in antiplane shear, utilizing the non-stiff
interface treatment in Erickson et al. (2022) and the high-order accurate
boundary-optimized SBP operators of Stiernström et al. (2023). SBP-preserving
interpolation allows the use of a lower-dimensional parameter representation
while maintaining dual consistency. In combination with RK4 time integration,
a fully discrete gradient expression is presented. Due to the self-adjointness
of RK4 Sanz-Serna (2016); Matsuda and Miyatake (2021), the space-time
discretization is dual consistent, meaning that the gradient expression is the
true gradient of the discrete misfit and consistent with the continuous
gradient. This claim is corroborated by numerical experiments in a dynamic
rupture setting, presented in Section 6.5.1.

We then proceed with inversions of model problems in increasingly complex
settings in Sections 6.5.2 \- 6.5.3. To start with, inverse-crime inversions
are performed. Here, we are able to reconstruct the direct effect parameter
and fault shear traction fairly accurately in the velocity-weakening region of
the fault with appreciable slip. In the velocity-strengthening region that has
negligible slip, the parameter sensitivity is close to zero, likely due to
weak signals. Next, inversions utilizing synthetic high-resolution data are
performed. For the direct effect parameter, a local minimum is found, which
exhibits oscillations in the parameter around the true value. Inversions for
shear traction interpolate the true value fairly well. We conclude that
regularization, likely in combination with filtering of the misfit and adjoint
source signals, is required to further improve the results.

The primary future research direction is thus to investigate suitable
regularizations of the misfit functional. In addition, we will extend the
adjoint framework to include other types of measurements, e. g. , strain rate
from fiber optics cables, pressure from hydrophones, and InSAR measurements of
surface displacement. Moreover, we will apply the inversion framework to
aseismic slow slip events.

## Acknowledgments

We sincerely thank Stéphanie Chaillat (ENSTA Paris, Stanford visiting
scholar), Laura Bagur and Alice Nassor (ENSTA Paris) for fruitful discussions
during the fall of 2022, on the treatment of displacement misfits in Section
5.2.

V. Stiernström was supported by the Swedish Research Council (grant no.
2017-04626). The computations were enabled by resources provided by the
National Academic Infrastructure for Supercomputing in Sweden (NAISS) at
UPPMAX partially funded by the Swedish Research Council through grant
agreement no. 2022-06725.

##  Appendix A Summation-by-parts properties

Integration by parts (IBP) is key in deriving energy balances and showing
conservation properties of the continuous equations. SBP operators allow us to
mimic the analysis in the discrete setting. It is shown in Almquist and Dunham
(2020, 2021) that $\mathbf{D}_{II}(\bm{\mu})$ satisfy SBP properties on
$\mathbf{\Omega}$, consistent with the continuous IBP property for
$\partial_{I}\mu\partial_{I}$ on $\Omega$. The results are summarized below.

For scalar functions $u$ and $v$, the variable coefficient Laplace operator
satisfies

$\left(v,\partial_{I}\mu\partial_{I}u\right)_{\Omega}=\left(v,\tau_{u}\right)_{\partial\Omega}-\left(\partial_{I}v,\mu\partial_{I}u\right)_{\Omega},$
(A.1)

where $\tau_{u}=\mu\partial_{\hat{n}}u$. Integrating once more yields

$\left(v,\partial_{I}\mu\partial_{I}u\right)_{\Omega}=\left(v,\tau_{u}\right)_{\partial\Omega}-\left(\tau_{v},u\right)_{\partial\Omega}+\left(\partial_{I}\mu\partial_{I}v,u\right)_{\Omega},$
(A.2)

In the discrete setting let $\mathbf{D}_{I}\approx\partial_{I}$ on
$\mathbf{\Omega}$ be a first-derivative SBP operator, based on the same SBP
norm $\mathbf{H}$ as $\mathbf{D}_{II}$. Then, $\mathbf{D}_{II}(\bm{\mu})$
satisfies the SBP property

$\left(\mathbf{v},\mathbf{D}_{II}(\bm{\mu})\mathbf{u}\right)_{\mathbf{\Omega}}=\left(\mathbf{v},\bm{\tau}_{u}\right)_{\bm{\partial\Omega}}-\left(\mathbf{D}_{I}\mathbf{v},\bar{\bar{\bm{\mu}}}\mathbf{D}_{I}\mathbf{u}\right)_{\mathbf{\Omega}}-\mathbf{v}^{T}\mathbf{R}_{\Delta}\mathbf{u},$
(A.3)

where $\bm{\tau}_{u}=\mathbf{T}\mathbf{u}$, with
$\mathbf{T}\approx\mu\partial_{\hat{n}}$ and
$\mathbf{R}_{\Delta}=\mathbf{R}_{\Delta}^{T}\geq 0$,
$\mathbf{R}_{\Delta}\approx 0$. Thus (A.3) is a discrete counterpart to the
IBP property (A.1). Applying the SBP property once more yields

$\left(\mathbf{v},\mathbf{D}_{II}(\bm{\mu})\mathbf{u}\right)_{\mathbf{\Omega}}=\left(\mathbf{v},\bm{\tau}_{u}\right)_{\bm{\partial\Omega}}-\left(\bm{\tau}_{v},\mathbf{u}\right)_{\bm{\partial\Omega}}+\left(\mathbf{D}_{II}(\bm{\mu})\mathbf{v},\mathbf{u}\right)_{\mathbf{\Omega}},$
(A.4)

For details on how $\mathbf{D}_{I}$ and $\mathbf{T}$ are constructed see
Almquist and Dunham (2020, 2021); Stiernström et al. (2023).

The SBP property (A.3) facilitates the derivation of a discrete energy rate
for scalar wave equations. This is achieved by replacing $\mathbf{v}$ by
$\dot{\mathbf{u}}$, in which case the volume terms
$\left(\cdot,\cdot\right)_{\mathbf{\Omega}}$, contribute to the discrete
potential energy rate. The discrete energy rates are then useful for guiding
the enforcement of boundary or interface conditions, such that an energy-
stable scheme may be obtained. See, e. g. , Duru et al. (2019); Wang and
Petersson (2019); Almquist and Dunham (2020, 2021); Erickson et al. (2022).
The property (A.4) states that $\mathbf{D}_{II}(\bm{\mu})$ augmented with
suitable boundary or interface treatment is self-adjoint with respect to
$\mathbf{H}_{\mathbf{\Omega}}$.

##  Appendix B Semi-discrete adjoint sources

To derive the semi-discrete gradient we seek an adjoint source
$\mathbf{Q}^{\dagger}$ satisfying

$\left(\mathbf{Q}^{\dagger},\frac{\partial\dot{\mathbf{u}}}{\partial\mathbf{p}}\right)_{\mathbf{\Omega}\times\mathcal{T}}=\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}}.$
(B.1)

Consider the gradient of the misfit (6.21) with a single receiver. For
convenience, we drop the superscripts. By the chain rule, it follows that

$\displaystyle\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}}$
$\displaystyle=\frac{\partial}{\partial\mathbf{p}}\left(\frac{1}{2}\int_{\mathcal{T}}\mathbf{r}^{2}dt\right)=\int_{\mathcal{T}}\mathbf{r}(t)\left(\hat{\bm{\delta}},\frac{\partial\mathbf{m}}{\partial\mathbf{p}}\right)_{\mathbf{\Omega}}dt=\left(\mathbf{r}\hat{\bm{\delta}},\frac{\partial\mathbf{m}}{\partial\mathbf{p}}\right)_{\mathbf{\Omega}\times\mathcal{T}}.$
(B.2)

If $\mathbf{m}=\dot{\mathbf{u}}$, we have
$\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}}=\left(\mathbf{r}\hat{\bm{\delta}},\frac{\partial\dot{\mathbf{u}}}{\partial\mathbf{p}}\right)_{\mathbf{\Omega}\times\mathcal{T}}$,
which is of the desired form (B.1). If $\mathbf{m}=\mathbf{u}$, introduce
$\hat{\mathbf{r}}(t)=\int_{0}^{t}\mathbf{r}(t^{\prime})dt^{\prime}+\hat{\mathbf{r}}_{0}$
with $\hat{\mathbf{r}}_{0}=-\int_{0}^{T}\mathbf{r}(t^{\prime})dt^{\prime}$ and
integrate
$\left(\mathbf{r}\hat{\bm{\delta}},\frac{\partial\mathbf{u}}{\partial\mathbf{p}}\right)_{\mathbf{\Omega}\times\mathcal{T}}$
by parts in time to obtain
$\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}}=\left(-\hat{\mathbf{r}}(t),\frac{\partial\dot{\mathbf{u}}}{\partial\mathbf{p}}\right)_{\mathbf{\Omega}\times\mathcal{T}}$.
Thus,

$\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}}=\left(\mathbf{S},\frac{\partial\dot{\mathbf{u}}}{\partial\mathbf{p}}\right)_{\mathbf{\Omega}\times\mathcal{T}},\quad\mathbf{S}(t)=\begin{cases}-\hat{\mathbf{r}}(t)\bm{\delta}(\bar{\mathbf{x}}-\bar{\mathbf{x}}_{r}),&\mathbf{m}=\mathbf{u},\\\
\mathbf{r}(t)\bm{\delta}(\bar{\mathbf{x}}-\bar{\mathbf{x}}_{r}),&\mathbf{m}=\dot{\mathbf{u}}.\end{cases}$
(B.3)

Summing $\mathbf{S}$ over all receivers results in $\mathbf{Q}^{\dagger}$ in
(6.31).

##  Appendix C Derivation of the semi-discrete gradient

To derive the gradient of (6.21), the procedure in Section 5 is followed.
Omitting external forcing (since it does not contribute to the gradient) the
Lagrangian misfit functional to (6.21) reads

$\displaystyle\bm{\mathcal{L}}=\bm{\mathcal{F}}$
$\displaystyle+\left(\dot{\mathbf{u}}^{\dagger},\bar{\bar{\bm{\rho}}}\ddot{\mathbf{u}}-\mathbf{L}(\mathbf{u},\mathbf{u}^{*},\bm{\tau}^{*})\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}$
(C.1)
$\displaystyle+\left(\dot{\mathbf{u}}^{\dagger},\bar{\bar{\bm{\rho}}}\ddot{\mathbf{u}}-\mathbf{L}(\mathbf{u},\mathbf{u}^{*},\bm{\tau}^{*})\right)_{\mathbf{\Omega}_{-}\times\mathcal{T}}$
$\displaystyle+\left(\bm{\Psi}^{\dagger},\dot{\bm{\Psi}}-\mathbf{G}\right)_{\bm{\Gamma}\times\mathcal{T}},$

where we omit subscripts $\pm$ of the grid functions since they are implied by
the inner products. Consider first
$\left(\dot{\mathbf{u}}^{\dagger},\mathbf{L}(\mathbf{u},\mathbf{u}^{*},\bm{\tau}^{*})\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}$.
By (6.6) it follows that

$\displaystyle\left(\dot{\mathbf{u}}^{\dagger},\mathbf{L}(\mathbf{u},\mathbf{u}^{*},\bm{\tau}^{*})\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}=$
(C.2)
$\displaystyle\left(\dot{\mathbf{u}}^{\dagger},\mathbf{D}_{II}(\bm{\mu})\mathbf{u}+\text{SAT}(\mathbf{u},\mathbf{u}^{*},\bm{\tau}^{*})\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}=$
$\displaystyle\left(\dot{\mathbf{u}}^{\dagger},\mathbf{D}_{II}(\bm{\mu})\mathbf{u}\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}+\left(\dot{\mathbf{u}}^{\dagger},\bm{\tau}^{*}-\bm{\tau}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}-\left(\dot{\bm{\tau}}^{\dagger},\mathbf{u}^{*}-\mathbf{u}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}.$

Swapping forward and adjoint variables yields

$\displaystyle\left(\dot{\mathbf{u}},\mathbf{L}(\mathbf{u}^{\dagger},\mathbf{u}^{\dagger*},\bm{\tau}^{\dagger*})\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}=$
(C.3)
$\displaystyle\left(\dot{\mathbf{u}},\mathbf{D}_{II}(\bm{\mu})\mathbf{u}^{\dagger}\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}+\left(\dot{\mathbf{u}},\bm{\tau}^{\dagger*}-\bm{\tau}^{\dagger}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}-\left(\dot{\bm{\tau}},\mathbf{u}^{\dagger*}-\mathbf{u}^{\dagger}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}.$

Next, using the SBP property (A.4) and IBP in time it follows that

$\displaystyle\left(\dot{\mathbf{u}}^{\dagger},\mathbf{D}_{II}(\bm{\mu})\mathbf{u}\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}=$
(C.4) $\displaystyle-$
$\displaystyle\left(\mathbf{D}_{II}(\bm{\mu})\mathbf{u}^{\dagger},\dot{\mathbf{u}}\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}+\left(\dot{\mathbf{u}}^{\dagger},\bm{\tau}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}-\left(\dot{\bm{\tau}}^{\dagger},\mathbf{u}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}},$

where the initial and terminal conditions in (6.19) and (6.31) were used to
cancel boundary terms in the integration. Due to the $p$-independence of data,
we assumed homogeneous initial conditions for the forward problem, since such
terms would vanish in the gradient expression regardless. Using (C.4) in (C.2)
together with (C.3) yields

$\displaystyle\left(\dot{\mathbf{u}}^{\dagger},\mathbf{L}(\mathbf{u},\mathbf{u}^{*},\bm{\tau}^{*})\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}=$
(C.5) $\displaystyle-$
$\displaystyle\left(\mathbf{D}_{II}(\bm{\mu})\mathbf{u}^{\dagger},\dot{\mathbf{u}}\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}+\left(\dot{\mathbf{u}}^{\dagger},\bm{\tau}^{*}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}-\left(\dot{\bm{\tau}}^{\dagger},\mathbf{u}^{*}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}=$
$\displaystyle-$
$\displaystyle\left(\dot{\mathbf{u}}^{\dagger},\mathbf{L}(\mathbf{u},\mathbf{u}^{*},\bm{\tau}^{*})\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}+BT$

where

$\displaystyle BT$
$\displaystyle=\left(\bm{\tau}^{\dagger*}-\bm{\tau}^{\dagger},\dot{\mathbf{u}}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}-\left(\mathbf{u}^{\dagger*}-\mathbf{u}^{\dagger},\dot{\bm{\tau}}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}$
(C.6)
$\displaystyle+\left(\dot{\mathbf{u}}^{\dagger},\bm{\tau}^{*}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}-\left(\dot{\bm{\tau}}^{\dagger},\mathbf{u}^{*}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}$
$\displaystyle=\left(\bm{\tau}^{\dagger*}-\bm{\tau}^{\dagger},\dot{\mathbf{u}}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}+\left(\dot{\mathbf{u}}^{\dagger*}-\dot{\mathbf{u}}^{\dagger},\bm{\tau}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}$
$\displaystyle+\left(\dot{\mathbf{u}}^{\dagger},\bm{\tau}^{*}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}+\left(\bm{\tau}^{\dagger},\dot{\mathbf{u}}^{*}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}$

with the last equality obtained by IBP in time. Next we perform a change of
variables in the term
$\left(\dot{\mathbf{u}}^{\dagger*}-\dot{\mathbf{u}}^{\dagger},\bm{\tau}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}$
from $\bm{\tau}$ to $\tilde{\bm{\tau}}$ using (6.8)

$\displaystyle\left(\dot{\mathbf{u}}^{\dagger*}-\dot{\mathbf{u}}^{\dagger},\bm{\tau}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}=\left(\dot{\mathbf{u}}^{\dagger*}-\dot{\mathbf{u}}^{\dagger},\tilde{\bm{\tau}}-\gamma(\mathbf{u}^{*}-\mathbf{u})\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}=$
(C.7)
$\displaystyle\left(\dot{\mathbf{u}}^{\dagger*}-\dot{\mathbf{u}}^{\dagger},\tilde{\bm{\tau}}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}+\left(\gamma(\mathbf{u}^{\dagger*}-\mathbf{u}^{\dagger}),\dot{\mathbf{u}}^{*}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}-\left(\gamma(\mathbf{u}^{\dagger*}-\mathbf{u}^{\dagger}),\dot{\mathbf{u}}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}},$

where IBP in time was used to obtain the last equality. Combining the result
with (C.6) and substituting $\tilde{\bm{\tau}}^{\dagger}$ results in

$\displaystyle BT$
$\displaystyle=\left(\bm{\tau}^{\dagger*}-\tilde{\bm{\tau}}^{\dagger},\dot{\mathbf{u}}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}+\left(\dot{\mathbf{u}}^{\dagger*}-\dot{\mathbf{u}}^{\dagger},\tilde{\bm{\tau}}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}$
(C.8)
$\displaystyle+\left(\dot{\mathbf{u}}^{\dagger},\bm{\tau}^{*}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}+\left(\tilde{\bm{\tau}}^{\dagger},\dot{\mathbf{u}}^{*}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}$

Next we use (6.7) and (6.24) to obtain

$\displaystyle\tilde{\bm{\tau}}$
$\displaystyle=\bm{\tau}^{*}+Z(\dot{\mathbf{u}}-\dot{\mathbf{u}}^{*})$ (C.9)
$\displaystyle\tilde{\bm{\tau}}^{\dagger}$
$\displaystyle=\bm{\tau}^{\dagger*}-Z(\dot{\mathbf{u}}^{\dagger}-\dot{\mathbf{u}}^{\dagger*})$

Substituting (C.9) into (C.8) and eliminating terms yields

$\displaystyle BT$
$\displaystyle=\left(\dot{\mathbf{u}}^{\dagger*},\bm{\tau}^{*}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}+\left(\bm{\tau}^{*\dagger},\dot{\mathbf{u}}^{*}\right)_{\bm{\partial\Omega}_{+}\times\mathcal{T}}$
(C.10)

Since the target values satisfy the boundary conditions, $BT$ vanishes on all
boundaries but the fault. This is straightforwardly seen on displacement or
traction boundaries, where $\mathbf{u}^{*}$ and $\mathbf{u}^{\dagger*}$, or
$\bm{\tau}^{*}$ and $\bm{\tau}^{\dagger*}$, respectively, are zero. On outflow
boundaries, substituting (6.9), (6.10), (6.25), (6.26) with
$\mathbf{R}=\mathbf{0}$ into (C.10) yields the result. Adding contributions of
(C.10) from $\mathbf{\Omega}_{-}$ yields

$IT=BT_{+}+BT_{-}=\left(\bm{\tau}^{\dagger*}_{+},\mathbf{V}^{*}\right)_{\bm{\Gamma}\times\mathcal{T}}-\left(\mathbf{V}^{\dagger*},\bm{\tau}^{*}_{+}\right)_{\bm{\Gamma}\times\mathcal{T}}.$
(C.11)

where it is assumed that $\mathbf{V}^{*}=[\\![\dot{\mathbf{u}}^{*}]\\!]$ and
$\mathbf{V}^{\dagger*}=-[\\![\dot{\mathbf{u}}^{\dagger*}]\\!]$.

By (C.1) and (C.5) it then follows that

$\displaystyle\bm{\mathcal{L}}=\bm{\mathcal{F}}$
$\displaystyle-\left(\bar{\bar{\bm{\rho}}}\ddot{\mathbf{u}}^{\dagger}-\mathbf{L}(\mathbf{u}^{\dagger},\mathbf{u}^{\dagger*},\bm{\tau}^{\dagger*}),\dot{\mathbf{u}}\right)_{\mathbf{\Omega}_{+}\times\mathcal{T}}$
(C.12)
$\displaystyle-\left(\bar{\bar{\bm{\rho}}}\ddot{\mathbf{u}}^{\dagger}-\mathbf{L}(\mathbf{u}^{\dagger},\mathbf{u}^{\dagger*},\bm{\tau}^{\dagger*}),\dot{\mathbf{u}}\right)_{\mathbf{\Omega}_{\times}\mathcal{T}-}$
$\displaystyle-ST-IT,$

where

$ST=\left(\dot{\bm{\Psi}}^{\dagger},\bm{\Psi}\right)_{\bm{\Gamma}\times\mathcal{T}}+\left(\bm{\Psi}^{\dagger},\mathbf{G}\right)_{\bm{\Gamma}\times\mathcal{T}}-\left(\bm{\Psi}_{0}^{\dagger},\bm{\Psi}_{0}\right)_{\bm{\Gamma}},$
(C.13)

is obtained from IBP in time, where the term at $t=T$ vanishes due to the
terminal condition in (6.31).

We will now derive the gradient
$\frac{\partial\bm{\mathcal{L}}}{\partial\mathbf{p}}$. Combining (B.1) with
the adjoint scheme (6.31), the gradient of the volume terms vanishes such that
$\frac{\partial\bm{\mathcal{L}}}{\partial\mathbf{p}}=-\frac{\partial
ST}{\partial\mathbf{p}}-\frac{\partial IT}{\partial\mathbf{p}}$. First,
consider the inversion for a parameter in either $F$ or $G$. Since

$\frac{\partial\bm{\tau}_{+}^{*}}{\partial\mathbf{p}}=-\left(\bar{\bar{\mathbf{F}}}_{V}\frac{\partial\mathbf{V}^{*}}{\partial\mathbf{p}}+\bar{\bar{\mathbf{F}}}_{\Psi}\frac{\partial\bm{\Psi}}{\partial\mathbf{p}}+\bar{\bar{\mathbf{F}}}_{p}\frac{\partial\mathbf{p}}{\partial\mathbf{p}}\right),$
(C.14)

and

$\frac{\partial
ST}{\partial\mathbf{p}}=\left(\dot{\bm{\Psi}}^{\dagger},\frac{\partial\bm{\Psi}}{\partial\mathbf{p}}\right)_{\bm{\Gamma}\times\mathcal{T}}+\left(\bm{\Psi}^{\dagger},\bar{\bar{\mathbf{G}}}_{V}\frac{\partial\mathbf{V}^{*}}{\partial\mathbf{p}}+\bar{\bar{\mathbf{G}}}_{\Psi}\frac{\partial\bm{\Psi}}{\partial\mathbf{p}}+\bar{\bar{\mathbf{G}}}_{p}\frac{\partial\mathbf{p}}{\partial\mathbf{p}}\right)_{\bm{\Gamma}\times\mathcal{T}},$
(C.15)

we obtain

$\displaystyle\frac{\partial IT}{\partial\mathbf{p}}$
$\displaystyle=\left(\bm{\tau}^{\dagger*}_{+},\frac{\partial\mathbf{V}^{*}}{\partial\mathbf{p}}\right)_{\bm{\Gamma}\times\mathcal{T}}+\left(\mathbf{V}^{\dagger*},\bar{\bar{\mathbf{F}}}_{V}\frac{\partial\mathbf{V}^{*}}{\partial\mathbf{p}}+\bar{\bar{\mathbf{F}}}_{\Psi}\frac{\partial\bm{\Psi}}{\partial\mathbf{p}}+\bar{\bar{\mathbf{F}}}_{p}\frac{\partial\mathbf{p}}{\partial\mathbf{p}}\right)_{\bm{\Gamma}\times\mathcal{T}}.$
(C.16)

and

$\displaystyle\frac{\partial ST}{\partial\mathbf{p}}+\frac{\partial
IT}{\partial\mathbf{p}}=$
$\displaystyle\left(\bm{\tau}^{\dagger*}_{+}+\bar{\bar{\mathbf{F}}}_{V}\mathbf{V}^{\dagger*}+\bar{\bar{\mathbf{G}}}_{V}\bm{\Psi}^{\dagger},\frac{\partial\mathbf{V}^{*}}{\partial\mathbf{p}}\right)_{\bm{\Gamma}\times\mathcal{T}}$
(C.17) $\displaystyle+$
$\displaystyle\left(\dot{\bm{\Psi}}^{\dagger}+\bar{\bar{\mathbf{F}}}_{\Psi}\mathbf{V}^{\dagger*}+\bar{\bar{\mathbf{G}}}_{\Psi}\bm{\Psi}^{\dagger},\frac{\partial\bm{\Psi}}{\partial\mathbf{p}}\right)_{\bm{\Gamma}\times\mathcal{T}}$
$\displaystyle+$
$\displaystyle\left(\bar{\bar{\mathbf{G}}}_{p}\bm{\Psi}^{\dagger}+\bar{\bar{\mathbf{F}}}_{p}\mathbf{V}^{\dagger*},\frac{\partial\mathbf{p}}{\partial\mathbf{p}}\right)_{\bm{\Gamma}\times\mathcal{T}}.$

Using the discrete friction law (6.27) and state evolution equation (6.30), we
arrive at

$\frac{\partial ST}{\partial\mathbf{p}}+\frac{\partial
IT}{\partial\mathbf{p}}=\left(\bar{\bar{\mathbf{G}}}_{p}\bm{\Psi}^{\dagger}+\bar{\bar{\mathbf{F}}}_{p}\mathbf{V}^{\dagger*},\frac{\partial\mathbf{p}}{\partial\mathbf{p}}\right)_{\bm{\Gamma}\times\mathcal{T}}.$
(C.18)

Similarly, if inverting for the initial state $\bm{\Psi}_{0}$

$\frac{\partial\bm{\tau}_{+}^{*}}{\partial\bm{\Psi}_{0}}=-\left(\bar{\bar{\mathbf{F}}}_{V}\frac{\partial\mathbf{V}^{*}}{\partial\bm{\Psi}_{0}}+\bar{\bar{\mathbf{F}}}_{\Psi}\frac{\partial\bm{\Psi}}{\partial\bm{\Psi}_{0}}\right),$
(C.19)

while

$\displaystyle\frac{\partial ST}{\partial\bm{\Psi}_{0}}$
$\displaystyle=\left(\dot{\bm{\Psi}}^{\dagger},\frac{\partial\bm{\Psi}}{\partial\bm{\Psi}_{0}}\right)_{\bm{\Gamma}\times\mathcal{T}}+\left(\bm{\Psi}^{\dagger},\bar{\bar{\mathbf{G}}}_{V}\frac{\partial\mathbf{V}^{*}}{\partial\bm{\Psi}_{0}}+\bar{\bar{\mathbf{G}}}_{\Psi}\frac{\partial\bm{\Psi}}{\partial\bm{\Psi}_{0}}\right)_{\bm{\Gamma}\times\mathcal{T}}$
(C.20)
$\displaystyle-\left(\bm{\Psi}_{0}^{\dagger},\frac{\partial\bm{\Psi}_{0}}{\partial\bm{\Psi}_{0}}\right)_{\bm{\Gamma}},$

such that

$\frac{\partial ST}{\partial\bm{\Psi}_{0}}+\frac{\partial
IT}{\partial\bm{\Psi}_{0}}=-\left(\bm{\Psi}_{0}^{\dagger},\frac{\partial\bm{\Psi}_{0}}{\partial\bm{\Psi}_{0}}\right)_{\bm{\Gamma}}$
(C.21)

Since (in indical notation)
$\frac{\partial\mathbf{p}_{j}}{\partial\mathbf{p}_{i}}=\delta_{ij}$ (and
similarly for $\bm{\Psi}_{0}$) (6.33) follows.

######  Remark.

To reach (C.18) and (C.21) we have assumed
$\mathbf{V}^{*}=[\\![\dot{\mathbf{u}}^{*}]\\!]$ and
$\mathbf{V}^{\dagger*}=-[\\![\dot{\mathbf{u}}^{\dagger*}]\\!]$. This is true
for $\mathbf{V}^{\dagger*}$, since the adjoint scheme satisfies (6.29)
exactly, due to $F^{\dagger}$ being linear in $V^{\dagger}$ such that the
solve is trivial. However, due to the non-linear solve required in (6.13),

$[\\![\dot{\mathbf{u}}^{*}]\\!]=\mathbf{V}^{*}+\varepsilon.$ (C.22)

When adding the contributions from the interface terms, this leads to

$IT=\left(\bm{\tau}^{\dagger*}_{+},\mathbf{V}^{*}+\varepsilon\right)_{\bm{\Gamma}\times\mathcal{T}}-\left(\mathbf{V}^{\dagger*},\bm{\tau}^{*}_{+}\right)_{\bm{\Gamma}\times\mathcal{T}}.$
(C.23)

such that

$\frac{\partial\bm{\mathcal{F}}^{(\varepsilon)}}{\partial\mathbf{p}_{i}}=\frac{\partial\bm{\mathcal{F}}}{\partial\mathbf{p}_{i}}+\left(\bm{\tau}^{\dagger*}_{+},\frac{\partial\varepsilon}{\partial\mathbf{p}_{i}}\right)_{\bm{\Gamma}\times\mathcal{T}}$
(C.24)

is the resulting gradient. Thus, the error in the non-linear solve results in
an error in the semi-discrete gradient. However, note that $F$ is invertible
and monotone, with a single root within the interval $[0,F^{-1}(\tau_{l})$]
Erickson et al. (2022). Therefore we can control $\varepsilon$ arbitrarily
through the tolerance in a bracketed root finding algorithm, e. g. bisection.

## References

  * Allison and Dunham [2018] K. L. Allison and E. M. Dunham.  Earthquake cycle simulations with rate-and-state friction and power-law viscoelasticity.  _Tectonophysics_ , 733:232–256, 2018.  doi:https://doi.org/10.1016/j.tecto.2017.10.021.  Physics of Earthquake Rupture Propagation.
  * Almquist and Dunham [2020] M. Almquist and E. M. Dunham.  Non-stiff boundary and interface penalties for narrow-stencil finite difference approximations of the Laplacian on curvilinear multiblock grids.  _J. Comput. Phys._ , 408:109294, 33, 2020.  doi:10.1016/j.jcp.2020.109294.
  * Almquist and Dunham [2021] M. Almquist and E. M. Dunham.  Elastic wave propagation in anisotropic solids using energy-stable finite differences with weakly enforced boundary and interface conditions.  _J. Comput. Phys._ , 424:109842, 33, 2021.  doi:10.1016/j.jcp.2020.109842.
  * Almquist et al. [2019] M. Almquist, S. Wang, and J. Werpers.  Order-preserving interpolation for summation-by-parts operators at nonconforming grid interfaces.  _SIAM J. Sci. Comput._ , 41(2):A1201–A1227, 2019\.  doi:10.1137/18M1191609.
  * Bader et al. [2023] M. Bader, M. Almquist, and E. M. Dunham.  Modeling and inversion in acoustic-elastic coupled media using energy-stable summation-by-parts operators.  _Geophysics_ , 88(3):T137–T150, 04 2023.  doi:10.1190/geo2022-0195.1.
  * Berg and Nordström [2012] J. Berg and J. Nordström.  Superconvergent functional output for time-dependent problems using finite differences on summation-by-parts form.  _J. Comput. Phys._ , 231(20):6846–6860, 2012\.  doi:10.1016/j.jcp.2012.06.032.
  * Butcher [2016] J. C. Butcher.  _Numerical methods for ordinary differential equations_.  John Wiley & Sons, Ltd., Chichester, third edition, 2016.  ISBN 978-1-119-12150-3.  doi:10.1002/9781119121534.  With a foreword by J. M. Sanz-Serna.
  * Dunham et al. [2011] E. M. Dunham, D. Belanger, L. Cong, and J. E. Kozdon.  Earthquake Ruptures with Strongly Rate-Weakening Friction and Off-Fault Plasticity, PartÂ 2: Nonplanar Faults.  _Bull. Seismol. Soc. Am._ , 101(5):2308–2322, 10 2011.  ISSN 0037-1106.  doi:10.1785/0120100076.
  * Duru and Dunham [2016] K. Duru and E. M. Dunham.  Dynamic earthquake rupture simulations on nonplanar faults embedded in 3D geometrically complex, heterogeneous elastic solids.  _J. Comput. Phys._ , 305:185–207, 2016.  doi:10.1016/j.jcp.2015.10.021.
  * Duru et al. [2019] K. Duru, K. L. Allison, M. Rivet, and E. M. Dunham.  Dynamic rupture and earthquake sequence simulations using the wave equation in second-order form.  _Geophys. J. Int._ , 219:796–815, 2019.  doi:10.1093/gji/ggz319.
  * Engquist and Majda [1977] B. Engquist and A. Majda.  Absorbing boundary conditions for the numerical simulation of waves.  _Math. Comp._ , 31(139):629–651, 1977.  doi:10.2307/2005997.
  * Erickson and Dunham [2014] B. A. Erickson and E. M. Dunham.  An efficient numerical method for earthquake cycles in heterogeneous media: Alternating subbasin and surface-rupturing events on faults crossing a sedimentary basin.  _J. Geophys. Res. Solid Earth_ , 119(4):3290–3316, 2014.  doi:https://doi.org/10.1002/2013JB010614.
  * Erickson et al. [2022] B. A. Erickson, J. E. Kozdon, and T. Harvey.  A non-stiff summation-by-parts finite difference method for the scalar wave equation in second order form: characteristic boundary conditions and nonlinear interfaces.  _J. Sci. Comput._ , 93(1):Paper No. 17, 39, 2022\.  doi:10.1007/s10915-022-01961-1.
  * Erickson et al. [2023] B. A. Erickson, J. Jiang, V. Lambert, S. D. Barbot, M. Abdelmeguid, M. Almquist, J. Ampuero, R. Ando, C. Cattania, A. Chen, L. Dal Zilio, S. Deng, E. M. Dunham, A. E. Elbanna, A. Gabriel, T. W. Harvey, Y. Huang, Y. Kaneko, J. E. Kozdon, N. Lapusta, D. Li, M. Li, C. Liang, Y. Liu, S. Ozawa, A. Perez-Silva, C. Pranger, P. Segall, Y. Sun, P. Thakur, C. Uphoff, Y. van Dinther, and Y. Yang.  Incorporating full elastodynamic effects and dipping fault geometries in community code verification exercises for simulations of earthquake sequences and aseismic slip (seas).  _Bull. Seismol. Soc. Am._ , 113(2):499–523, 01 2023.  doi:10.1785/0120220066.
  * Fang and Dunham [2013] Z. Fang and E. M. Dunham.  Additional shear resistance from fault roughness and stress levels on geometrically complex faults.  _J. Geophys. Res. Solid Earth_ , 118(7):3642–3654, 2013.  doi:https://doi.org/10.1002/jgrb.50262.
  * Fichtner [2011] A. Fichtner.  _Full Seismic Waveform Modelling and Inversion_.  Springer Berlin, Heidelberg, 2011.  doi:10.1007/978-3-642-15807-0.
  * Gallovič et al. [2019a] F. Gallovič, L. Valentová, J.-P. Ampuero, and A.-A. Gabriel.  Bayesian dynamic finite-fault inversion: 1. method and synthetic test.  _J. Geophys. Res.-Sol. Ea._ , 124(7):6949–6969, 2019a.  doi:10.1029/2019JB017510.
  * Gallovič et al. [2019b] F. Gallovič, L. Valentová, J.-P. Ampuero, and A.-A. Gabriel.  Bayesian dynamic finite-fault inversion: 2. application to the 2016 Mw 6.2 Amatrice, Italy, earthquake.  _J. Geophys. Res.-Sol. Ea._ , 124(7):6970–6988, 2019b.  doi:10.1029/2019JB017512.
  * Harris et al. [2009] R. A. Harris, M. Barall, R. Archuleta, E. Dunham, B. Aagaard, J. P. Ampuero, H. Bhat, V. Cruz-Atienza, L. Dalguer, P. Dawson, S. Day, B. Duan, G. Ely, Y. Kaneko, Y. Kase, N. Lapusta, Y. Liu, S. Ma, D. Oglesby, K. Olsen, A. Pitarka, S. Song, and E. Templeton.  The scec/usgs dynamic earthquake rupture code verification exercise.  _Seismol. Res. Lett._ , 80(1):119–126, 01 2009\.  doi:10.1785/gssrl.80.1.119.
  * Harris et al. [2018] R. A. Harris, M. Barall, B. Aagaard, S. Ma, D. Roten, K. Olsen, B. Duan, D. Liu, B. Luo, K. Bai, J.-P. Ampuero, Y. Kaneko, A. Gabriel, K. Duru, T. Ulrich, S. Wollherr, Z. Shi, E. Dunham, S. Bydlon, Z. Zhang, X. Chen, S. N. Somala, C. Pelties, J. Tago, V. M. Cruz-Atienza, J. Kozdon, E. Daub, K. Aslam, Y. Kase, K. Withers, and L. Dalguer.  A suite of exercises for verifying dynamic earthquake rupture codes.  _Seismol. Res. Lett._ , 89(3):1146–1162, 04 2018\.  doi:10.1785/0220170222.
  * Harvey et al. [2023] T. W. Harvey, B. A. Erickson, and J. E. Kozdon.  A high-order accurate summation-by-parts finite difference method for fully-dynamic earthquake sequence simulations within sedimentary basins.  _J. Geophys. Res. Solid Earth_ , 128(4):e2022JB025357, 2023.  doi:10.1029/2022JB025357.
  * Heimisson et al. [2019] E. R. Heimisson, E. M. Dunham, and M. Almquist.  Poroelastic effects destabilize mildly rate-strengthening friction to generate stable slow slip pulses.  _J. Mech. Phys. Solids_ , 130:262–279, 2019.  doi:10.1016/j.jmps.2019.06.007.
  * Hicken and Zingg [2011] J. E. Hicken and D. W. Zingg.  Superconvergent functional estimates from summation-by-parts finite-difference discretizations.  _SIAM J. Sci. Comput._ , 33(2):893–922, 2011\.  doi:10.1137/100790987.
  * Hicken and Zingg [2014] J. E. Hicken and D. W. Zingg.  Dual consistency and functional accuracy: a finite-difference perspective.  _J. Comput. Phys._ , 256:161–182, 2014.  doi:10.1016/j.jcp.2013.08.014.
  * Hirahara and Nishikiori [2019] K. Hirahara and K. Nishikiori.  Estimation of frictional properties and slip evolution on a long-term slow slip event fault with the ensemble Kalman filter: Numerical experiments.  _Geophys. J. Int._ , 219(3):2074–2096, 2019.  doi:10.1093/gji/ggz415.
  * Kano et al. [2013] M. Kano, S. Miyazaki, K. Ito, and K. Hirahara.  An adjoint data assimilation method for optimizing frictional parameters on the afterslip area.  _Earth, Planets and Space_ , 65(12):1575–1580, 2013.  doi:10.5047/eps.2013.08.002.
  * Kano et al. [2015] M. Kano, S. Miyazaki, Y. Ishikawa, Y. Hiyoshi, K. Ito, and K. Hirahara.  Real data assimilation for optimization of frictional parameters and prediction of afterslip in the 2003 Tokachi-oki earthquake inferred from slip velocity by an adjoint method.  _Geophys. J. Int._ , 203(1):646–663, 2015.  doi:10.1093/gji/ggv289.
  * Kano et al. [2020] M. Kano, S. Miyazaki, Y. Ishikawa, and K. Hirahara.  Adjoint-based direct data assimilation of GNSS time series for optimizing frictional parameters and predicting postseismic deformation following the 2003 Tokachi-oki earthquake.  _Earth, Planets and Space_ , 72(1), 2020.  doi:10.1186/s40623-020-01293-0.
  * Kozdon and Wilcox [2016] J. E. Kozdon and L. C. Wilcox.  Stable coupling of nonconforming, high-order finite difference methods.  _SIAM J. Sci. Comput._ , 38(2):A923–A952, 2016\.  doi:10.1137/15M1022823.
  * Kozdon et al. [2013] J. E. Kozdon, E. M. Dunham, and J. Nordström.  Simulation of dynamic earthquake ruptures in complex geometries using high-order finite difference methods.  _J. Sci. Comput._ , 55(1):92–124, 2013.  doi:10.1007/s10915-012-9624-5.
  * Lapusta et al. [2000] N. Lapusta, J. R. Rice, Y. Ben-Zion, and G. Zheng.  Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate- and state-dependent friction.  _J. Geophys. Res. Solid Earth_ , 105(B10):23765–23789, 2000.  doi:https://doi.org/10.1029/2000JB900250.
  * Lundquist et al. [2018] T. Lundquist, A. Malan, and J. Nordström.  A hybrid framework for coupling arbitrary summation-by-parts schemes on general meshes.  _J. Comput. Phys._ , 362:49–68, 2018.  doi:10.1016/j.jcp.2018.02.018.
  * Lundquist et al. [2020] T. Lundquist, A. G. Malan, and J. Nordström.  Efficient and error minimized coupling procedures for unstructured and moving meshes.  _J. Comput. Phys._ , 406:109158, 21, 2020.  doi:10.1016/j.jcp.2019.109158.
  * Marone [1998] C. Marone.  Laboratory-derived friction laws and their application to seismic faulting.  _Annu. Rev. Earth Planet. Sci._ , 26(1):643–696, 1998.  doi:10.1146/annurev.earth.26.1.643.
  * Matsuda and Miyatake [2021] T. Matsuda and Y. Miyatake.  Generalization of partitioned Runge-Kutta methods for adjoint systems.  _J. Comput. Appl. Math._ , 388:Paper No. 113308, 7, 2021\.  doi:10.1016/j.cam.2020.113308.
  * Mattsson and Carpenter [2010] K. Mattsson and M. H. Carpenter.  Stable and accurate interpolation operators for high-order multiblock finite difference methods.  _SIAM J. Sci. Comput._ , 32(4):2298–2320, 2010\.  doi:10.1137/090750068.
  * Nissen et al. [2015] A. Nissen, K. Kormann, M. Grandin, and K. Virta.  Stable difference methods for block-oriented adaptive grids.  _J. Sci. Comput._ , 65(2):486–511, 2015.  doi:10.1007/s10915-014-9969-z.
  * O’Reilly et al. [2015] O. O’Reilly, J. Nordström, J. E. Kozdon, and E. M. Dunham.  Simulation of earthquake rupture dynamics in complex geometries using coupled finite difference and finite volume methods.  _Commun. Comput. Phys._ , 17(2):337–370, 2015\.  doi:10.4208/cicp.111013.120914a.
  * Petersson and Sjögreen [2009] N. A. Petersson and B. Sjögreen.  An energy absorbing far-field boundary condition for the elastic wave equation.  _Commun. Comput. Phys._ , 6(3):483–508, 2009\.  doi:10.4208/cicp.2009.v6.p483.
  * Petersson et al. [2016] N. A. Petersson, O. O’Reilly, B. Sjögreen, and S. Bydlon.  Discretizing singular point sources in hyperbolic wave propagation problems.  _J. Comput. Phys._ , 321:532–555, 2016.  doi:10.1016/j.jcp.2016.05.060.
  * Peyrat and Olsen [2004] S. Peyrat and K. B. Olsen.  Nonlinear dynamic rupture inversion of the 2000 Western Tottori, Japan, earthquake.  _Geophys. Res. Let._ , 31(5), 2004.  doi:10.1029/2003GL019058.
  * Premus et al. [2022] J. Premus, F. Gallovič, and J.-P. Ampuero.  Bridging time scales of faulting: From coseismic to postseismic slip of the Mw 6.0 2014 South Napa, California earthquake.  _Sci. Adv._ , 8(38):eabq2536, 2022.  doi:10.1126/sciadv.abq2536.
  * Rice [1993] J. R. Rice.  Spatio-temporal complexity of slip on a fault.  _J. Geophys. Res. Solid Earth_ , 98(B6):9885–9907, 1993.  doi:https://doi.org/10.1029/93JB00191.
  * Rice et al. [2001] J. R. Rice, N. Lapusta, and K. Ranjith.  Rate and state dependent friction and the stability of sliding between elastically deformable solids.  _J Mech Phys Solids_ , 49(9):1865–1898, 2001\.  doi:10.1016/S0022-5096(01)00042-4.
  * Samareh [2001] J. A. Samareh.  Survey of shape parameterization techniques for high-fidelity multidisciplinary shape optimization.  _AIAA Journal_ , 39(5):877–884, 2001.  doi:10.2514/2.1391.
  * Sanz-Serna [2016] J. M. Sanz-Serna.  Symplectic Runge-Kutta schemes for adjoint equations, automatic differentiation, optimal control, and more.  _SIAM Rev._ , 58(1):3–33, 2016.  doi:10.1137/151002769.
  * Stiernström and Almquist [2022] V. Stiernström and M. Almquist.  Order-preserving non-conforming grid interfaces for boundary-optimized summation-by-parts operators.  In _The 15th International Conference on Mathematical and Numerical Aspects of Wave Propagation, Book of Abstracts, Version: July 17th_ , 2022.
  * Stiernström et al. [2023] V. Stiernström, M. Almquist, and K. Mattsson.  Boundary-optimized summation-by-parts operators for finite difference approximations of second derivatives with variable coefficients.  _J. Comput. Phys._ , 491:Paper No. 112376, 2023.  doi:10.1016/j.jcp.2023.112376.
  * Torberntsson et al. [2018] K. Torberntsson, V. Stiernström, K. Mattsson, and E. M. Dunham.  A finite difference method for earthquake sequences in poroelastic solids.  _Comput. Geosci._ , 22(5):1351–1370, 2018.  doi:10.1007/s10596-018-9757-1.
  * Tromp et al. [2005] J. Tromp, C. Tape, and Q. Liu.  Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels.  _Geophys. J. Int._ , 160(1):195–216, 2005.  doi:10.1111/j.1365-246X.2004.02453.x.
  * van Dinther et al. [2019] Y. van Dinther, H. R. Künsch, and A. Fichtner.  Ensemble data assimilation for earthquake sequences: probabilistic estimation and forecasting of fault stresses.  _Geophys. J. Int._ , 217(3):1453–1478, 2019.  doi:10.1093/gji/ggz063.
  * Wang and Petersson [2019] S. Wang and N. A. Petersson.  Fourth order finite difference methods for the wave equation with mesh refinement interfaces.  _SIAM J. Sci. Comput._ , 41(5):A3246–A3275, 2019\.  doi:10.1137/18M1211465.