Köhne
et al. Earth, Planets and Space (2025) 77:3
https://doi.org/10.1186/s40623-024-02121-5
FULL PAPER
Open Access
© The Author(s) 2024.
Open Access
This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the
original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or
other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this
licence, visit
http:// creat iveco mmons. org/ licen ses/ by/4. 0/
.
Earth, Planets and Space
Probabilistic estimation of rheological
properties in subduction zones using sequences
of earthquakes and aseismic slip
Tobias Köhne
1*
, Rishav Mallick
2
and Mark Simons
1
Abstract
Constraining the effective rheology of major faults contributes to improving our understanding of the physics
of plate boundary deformation. Geodetic observations over the earthquake cycle are often used to estimate key
rheological parameters, assuming specific laboratory-informed classes of viscous or frictional rheological models.
However, differentiating between various rheological model classes using only observations of a single earthquake
(coseismic and postseismic deformation) is difficult—especially in the presence of coarse spatiotemporal sampling,
inherent observational noise, and non-uniqueness of the inverted properties. In this study, we present a framework
to estimate key rheological parameters of a subduction zone plate interface using simulations of sequences of earth-
quakes and aseismic slip, constrained by pre- and postseismic surface displacement timeseries. Our simplified forward
model consists of a two-dimensional subduction zone, represented by a discretized planar fault or narrow shear zone,
divided into a locked, shallow region (“asperity”) experiencing periodically imposed coseismic events, and a stress-
driven creeping section governed by power-law viscoelasticity or rate-dependent friction. Our inverse model fits
the rheological parameters of the interface to surface displacement timeseries in a Bayesian probabilistic way. We
validate that our proposed framework can successfully recover depth-dependent profiles of effective viscosity
using a synthetic dataset of pre- and postseismic observations. Our first set of numerical experiments show that our
framework is only mildly sensitive to uncertainties in the rupture history or assumed coseismic slip, making it robust
enough to be applied to real observations of subduction zones. Our second set of tests considers the similarities
of surface displacement timeseries between synthetic models that model the plate interface either as a shear zone
described by power-law viscosity, or a surface described by rate-dependent friction. Here, we find that the ability to fit
surface observations using functional or mechanical models assuming frictional behavior does not constitute suf-
ficient evidence to actually infer frictional behavior at depth, as the surface expressions are virtually indistinguishable
from deformation generated from models with depth-variable power-law viscous behavior. Based on our numerical
experiments, we conclude that studies that aim to infer the mechanical behavior and rheological properties at depth
in subduction zones should consider the surface expression from time periods representative of the entire seismic
cycle.
Keywords
Rheology and friction of fault zones, Subduction zone processes, Bayesian inference, Plate motions,
Satellite geodesy
*Correspondence:
Tobias Köhne
tkoehne@caltech.edu
Full list of author information is available at the end of the article
Page 2 of 19
Köhne
et al. Earth, Planets and Space (2025) 77:3
Graphical Abstract
1 Introduction
Constraining the effective rheology of subduction zone
megathrusts is crucial to our understanding of the phys
-
ics of convergent plate boundary deformation (e.g.,
Bürgmann and Dresen
2008
). Commonly adopted rheo
-
logical descriptions of the megathrust come from labo
-
ratory experiments conducted on a limited number of
rock types, at length scales varying from sub-microns to
meters (e.g., Blanpied et al.
1995
; Scholz
1998
; Marone
1998
; Hirth
2002
; Hirth and Kohlstedt
2004
). However, it
is not obvious how to correctly upscale laboratory-cali
-
brated rheologies to the geological scale, given that the
role of material heterogeneities in rocks and interactions
between physical (and chemical) processes at various
scales remain poorly understood (Handy
1994
; Yamash
-
ita et al.
2015
; Fagereng and Beall
2021
; Bercovici et al.
2023
). Surface geodetic measurements are sensitive to a
spatially and process averaged
effective
rheology, with the
averaging scale comparable to the distance between the
deforming source and the surface, and the specific form
of averaging ambiguous. Regardless of the source of any
proposed rheology, postseismic displacement timeseries
observations of GNSS-based networks near plate inter
-
faces can be used to estimate ranges of parameters for
rheological models. However, it is unclear if geodetic
observations can distinguish between different models at
megathrust scales.
In the case of Northeastern Japan, several constitu
-
tive laws have been used to model surface postseismic
deformation in the wake of the 2011
M
w
9.1
Tohoku-oki
earthquake. For example, Sun et al. (
2014
) modeled the
postseismic displacements using only flow in the viscoelas
-
tic mantle described by a Burgers rheology (Müller
1986
;
Hetland and Hager
2005
). Subsequent improvements in
the methodology of the inverse problem are highlighted in
the works of Freed et al. (
2017
), who estimated the after
-
slip on the fault kinematically after modeling the viscoe
-
lastic mantle using a linear viscous rheology, and Fukuda
and Johnson (
2021
), who jointly solve for stress-driven
frictional afterslip and viscoelastic deformation as a Burg
-
ers body. Both Agata et al. (
2019
) and Muto et al. (
2019
)
invoke non-steady-state or transient viscous rheologies
in the mantle, in the form of power-law Burgers bodies
(Chopra
1997
; Masuti et al.
2016
), combined with fric
-
tional afterslip on the plate interface. Notably, none of the
aforementioned studies attempted to verify whether their
model setup and parameters would also sufficiently repro
-
duce the preseismic deformation phase.
Based on the wide range of plausible models matching
postseismic observations, it is therefore useful to study
classes of rheological models that are consistent with
the
entire
interseismic period (including the post- and
preseismic time periods) using a
single
set of parameter
values (Hearn and Thatcher
2015
; Mallick et al.
2022
).
Such models may provide additional constraints on
which structural features or what ranges of rheological
parameters are indeed required by the data, as opposed
to being unconstrained free parameters that may lead
to overfitting of the observations. To this end, we build
on the formulation of Hetland et al. (
2010
); Hetland and
Simons (
2010
); Kanda et al. (
2013
) to simulate interseis
-
mic creep in an idealized subduction zone. In this model,
we create sequences of earthquakes and aseismic slip
(SEAS) by applying coseismic slip periodically on an oth
-
erwise locked zone (asperity) and letting the rest of the
system mechanically evolve following a predefined rheo
-
logical model. Conveniently, earthquake cycle models
become independent of the initial state once fully spun-
up (i.e., the simulation becomes cycle-invariant), avoid
-
ing the large influence these initial conditions impart
(e.g., Agata et al.
2019
; Muto et al.
2019
). They are, of
course, still sensitive to varying extent on other assump
-
tions, such as the fault geometry or the rupture history.
The purpose of this study is twofold. First, using our
proposed earthquake cycle framework, we investigate
the degree to which postseismic relaxation models are
Page 3 of 19
Köhne
et al. Earth, Planets and Space (2025) 77:3
able to differentiate between different constitutive laws
(i.e., frictional versus viscous behavior) or infer spatial
heterogeneity, given surface displacement timeseries.
Second, we explore the sensitivity of our framework with
respect to various assumptions that will become relevant
when applied to real surface observations. In Sect.
2
, we
provide a short, technical overview of our method, leav
-
ing the numerical details for Appendices
A
–
D
. Within
this framework, we then create a simplified model of a
two-dimensional subduction zone which we introduce
in Sect.
3
. We then explore our model through various
numerical experiments, and show inversion results in
Sect.
4
. Finally, we discuss the implications of our find
-
ings in terms of the recoverability of rheological param
-
eters and fault structure in Sect.
5
.
2 Proposed framework
Our framework is a continuation of the earthquake cycle
simulations in subduction zones proposed by Hetland
et al. (
2010
); Hetland and Simons (
2010
); Kanda et al.
(
2013
). Figure
1
summarizes the general workflow of our
proposed framework, which is broadly divided into a for
-
ward and an inverse part. The appendices describe the
individual components in detail.
The forward model starts from the Elastic Subduct
-
ing Plate Model (ESPM, Kanda and Simons
2010
), a
natural extension of the popular Backslip Model (BSM,
Savage
1983
). The ESPM approximates the subducting
plate as two parallel dislocation surfaces, each defining
an interface to the surrounding bulk material (see Fig.
2
and Appendix
A
). The two dislocation surfaces allow the
Fig. 1
Workflow schematic, separated into the forward and inverse parts as shown by the dashed gray line, and described in detail in Sect.
2
and Appendices
A
–
D
Fig. 2
Fault setup for our simulations, based on the Elastic Subducting Plate Model (ESPM, Kanda and Simons
2010
). The downgoing slab
is modeled by two interfaces bounding an effective elastic plate of constant thickness
H
. The upper and lower interfaces experience reverse
and normal motion, respectively. The upper interface is furthermore assumed fully coupled (i.e., locked) from the trench down to a depth of
D
shallow
.
Along the non-coupled upper interface, shear resistance is described using rheological models (see Appendix
B
) with depth-varying parameters
with knickpoints at
D
middle
and
D
deep
. Displacement timeseries are evaluated at discrete locations along the free surface and labeled S1–5. The
length scales, observer coverage, and dipping angle are loosely based on the Northern Japan subduction zone
Page 4 of 19
Köhne
et al. Earth, Planets and Space (2025) 77:3
ESPM to maintain physical consistency over multiple
earthquake cycles, i.e., the incoming plate subducts with
-
out net strain accumulation after a completed rupture
cycle. We assume all deformation on the plate interfaces
and in the bulk can effectively be collapsed onto the fault
planes, based on numerical simulations of the extent of
shear localization with power-law rheological param
-
eters (Moore and Parsons
2015
; Mallick et al.
2022
).
The upper fault interface can readily be separated into
a shallow zone that only slip coseismically (an asperity),
and the rest of the surface which creeps steadily in the
interseismic period but has negligible slip during earth
-
quake ruptures. This separation has been shown to be a
good approximation of the kinematics in the Northern
Japan region (e.g., Iinuma et al.
2016
; Philibosian and
Meltzner
2020
). The shallow asperity is forced to exhibit
(quasi-)periodic ruptures, a behavior predicted by elastic
rebound theory and substantiated for subduction zones
across the world (e.g., Griffin et al.
2020
; Philibosian and
Meltzner
2020
). To avoid numerical singularities and to
acknowledge that asperities may exhibit some transition
regions, the coseismic slip is tapered into the surround
-
ing area.
Between ruptures, the creeping surface evolves natu
-
rally based on a prescribed constitutive law, the coseis
-
mic stress changes, and the imposition of far-field plate
convergence. The first constitutive law explored here
is rate-dependent friction (Hetland et al.
2010
; Kanda
et al.
2013
) (also see Appendix
B.2
). Here, the shear
stress on the frictional interface can be modeled as
τ
∝
(
a
−
b
)σ
E
log
(
v
/
v
0
)
, where
a
−
b
≥
0
is the rate-
strengthening friction parameter,
σ
E
is the effective nor
-
mal stress acting on the interface,
v
is the sliding velocity,
and
v
0
is an arbitrary reference velocity. The second law
we explore is power-law viscosity (Hetland et al.
2010
)
(also see Appendix
B.1
). In this model, the shear stress
can be written as
τ
n
=
α
n
v
, where
α
n
and
n
are the two
parameters describing the material behavior. Regardless
of the constitutive law, we allow the rheological param
-
eters to vary with depth. Assuming a linearly elastic bulk
(for computational reasons), we then derive kernels using
Boundary Element Methods which relate slip and stress
at discretized fault patches to stress at all other patches
as well as to displacements at the free surface (Davis
2017
). Given the constitutive laws and the forcing of the
system, the system can be integrated in time using stand
-
ard numerical integration routines (e.g., SciPy, Virtanen
et al.
2020
). In order to obtain independence from the
initial conditions, the system is integrated over multiple
earthquake cycles, and only the output from the last cycle
is used. This concludes the forward model (left side in
Fig.
1
, also see Appendix
C
).
In the inverse model, we aim to derive the rheological
parameters of the constitutive laws considered that best
fit given (synthetic) surface observations. To this end, we
compute the residual timeseries between the predicted
displacement computed from the forward model using
candidate rheological properties and the observed sur
-
face displacement. The residuals are then passed on to
a Bayesian inversion scheme, enabling the probabilistic
estimation of the rheological parameters (e.g., Minson
et al.
2013
) (right side in Fig.
1
, also see Appendix
D
).
3 Model setup
In this section, we apply our framework to the modeling
of a simplified, two-dimensional subduction zone, exhib
-
iting power-law viscous or rate-dependent frictional
behavior.
3.1 Fault geometry
Our 2D fault geometry directly follows the ESPM (Fig.
2
,
also see Appendix
A
). The shallow part of the upper plate
interface is locked and only slips coseismically (the asper
-
ity). The rest of the upper plate interface creeps inter
-
seismically as a response to the stress induced by the
imposed coseismic slip in the asperity and the secular
loading rate, with the character of the response a function
of the chosen rheological model. The lower plate inter
-
face is assumed to be creeping at the plate convergence
rate, which is a convenient way to ensure the long-term
vertical surface motion is consistent with a steady-state
earthquake cycle, and that the incoming plate actually
subducts on geologic time scales (Kanda and Simons
2010
). The entire fault is assumed to be embedded in a
linear-elastic half-space, and the creeping part of the
upper plate interface is simulated through a collection of
line dislocation elements (“patches”). For the rest of this
study, we will omit the
upper
and
lower
qualifiers for the
plate interface for conciseness, and
plate interface
will
always refer to the upper plate interface unless otherwise
noted.
3.2
Synthetic target observations
We start our exploration of the framework with a refer
-
ence case assuming power-law viscosity (
τ
n
=
α
n
v
) as
the rheological model, which can be used to approxi
-
mate a wide variety of behavior, including friction and
linear viscosity, depending on the exponent
n
(Montési
and Hirth
2003
; Montési
2004
; Mallick et al.
2022
). The
rheological properties we (somewhat arbitrarily) choose
for our target model are based on the common assump
-
tion that the plate interface in subduction zones can be
approximated by a frictional (
n
→∞
), relatively weak
shallow interface, strengthening towards intermediate
Page 5 of 19
Köhne
et al. Earth, Planets and Space (2025) 77:3
depths, and then further transitioning into a weaker,
viscous interface at larger depths (
n
≤
4
) (see Table
1
for values of the rheological parameters). The resulting
profile of the
effective
viscosity, derived using Eq. (
B.2
),
is shown in Fig.
5
a. For the rupture history, we choose
an earthquake slip of 9.5 m with a recurrence time
interval of 100 years (approximately the convergence
rate between the Pacific Plate and Northern Japan) up
to a depth of 35 km. The generated surface displace
-
ments
d
target
are shown in Fig.
3
, covering a time period
of approximately 15 years before the earthquake and
10 years after the earthquake. For the inversions, we
add normally distributed noise of a standard deviation
of 10 mm to create the synthetic, noisy observations.
The goal of the inversion is to recover the rheological
parameters (including their depth-dependence) used to
generate the synthetic surface observations.
3.3 Fitted parameters
In the inverse model for the reference case, we aim to fit the
power-law viscosity parameters
α
eff
,
n
(see Appendix
B.1
for a mathematical description of the parameters) at three
different depths, as well as the depths themselves, fully
defining the viscosity profile that was used to generate the
target displacement timeseries. The three different depths
are denoted “shallow”, “middle”, and “deep”, leading to the
following parameter set:
We do not estimate
D
shallow
, corresponding to the lock
-
ing depth, which we keep constant to not have to recre
-
ate the fault mesh and kernels at every sampling step.
In later inversions, we estimate the rate-dependent
frictional strength
α
h
=
(
a
−
b
)σ
E
at different depths
instead of
α
eff
,
n
(see Appendix
B.2
for a full mathemati
-
cal description).
4 Inversion results
The output of the Bayesian inversion routine is the pos
-
terior sample distribution of each parameter in the set
θ
, which approximates the posterior probability density
function of each parameter. Derived from
θ
is the result
-
ing depth-varying profile of the effective viscosity. We
present the results of several experiments, in the form
of posterior probability densities (revealing parameter
covariances), as well as plots of the posterior error on
effective viscosity with depth. Since we know truth in
these synthetic experiments, we can quantify the quality
of fit to both the target (synthetic) surface displacement
timeseries, and the target underlying viscosity profiles, by
defining two error metrics,
δ
α
(
t
i
)
and
δ
d
, the average vis
-
cosity and surface displacement errors, respectively:
where
t
is time,
i
is the index of the
N
t
observations,
j
is
the index of the
N
o
observer stations at the surface,
c
is
one of the horizontal or vertical data components
H,V
,
respectively,
k
is the index of the
N
s
MCMC samples,
m
is
the index of the
N
p
simulated patches,
l
m
is the patch
(1)
θ
=
{
α
eff,shallow
,
α
eff,middle
,
α
eff,deep
,
n
shallow
,
n
middle
,
n
deep
,
D
middle
,
D
deep
}
.
(2)
δ
α
(
t
i
)
=
1
L
N
p
∑
m
=
1
l
m
(
1
N
s
N
s
∑
k
=
1
(
log
10
α
pred
eff,
m
,
k
(
t
i
)
−
log
10
α
target
eff,
m
(
t
i
)
)
2
)
1
/
2
,
(3)
δ
d
=
1
2
N
t
N
o
N
t
∑
i
=
1
N
o
∑
j
=
1
∑
c
=
H,V
(
1
N
s
N
s
∑
k
=
1
(
d
pred
i
,
j
,
c
,
k
−
d
target
i
,
j
,
c
)
2
)
1
/
2
,
Table 1
Rheological parameters
α
eff
,
n
used to create the
synthetic target dataset of surface displacement timeseries, and
their respective knickpoint depths
D
. Note that
α
eff
is defined in
(
B.2
) and depends on a reference slip velocity
Shallow
Middle
Deep
D
[
km
]
35
50
100
α
eff
[
Pas
/
m
]
10
12
10
15
10
13
n
[
−
]
10
6
1
Fig. 3
Simulated surface displacement timeseries used as input
for the reference case. The solid lines are the timeseries
d
target
,
with colors corresponding to observer location. The labels refer
to the names of the stations, with the range of trenchward distances
given as well. The black dots are the synthetic observations which
include a 10 mm standard deviation Gaussian noise. At approx.
14.2 years (vertical gray line), an earthquake occurs and starts
a postseismic transient process. The coseismic offset (including
the effect of tapered slip) are both removed in this plot, and in the
observations used in the inversion