# The Community Code Verification Exercise for Simulating Sequences of Earthquakes and Aseismic Slip (SEAS)

Brittany Erickson, Junle Jiang, Michael Barall, Nadia Lapusta, Eric Dunham, Ruth Harris, Lauren Abrahams, Kali Allison, Jean Paul Ampuero, Sylvain Barbot, et al.

► To cite this version:

Brittany Erickson, Junle Jiang, Michael Barall, Nadia Lapusta, Eric Dunham, et al.. The Community Code Verification Exercise for Simulating Sequences of Earthquakes and Aseismic Slip (SEAS). Seismological Research Letters, 2020, 91 (2A), pp.874-890. 10.1785/0220190248. hal-03578332

HAL Id: hal-03578332

https://hal.science/hal-03578332v1

Submitted on 9 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
**This manuscript is a preprint** and has been submitted for publication in *Seismological Research Letters*. Please note that the manuscript has yet to be peer-reviewed and has not been formally accepted for publication. Subsequent versions of this manuscript may have slightly different content. If accepted, the final version of this manuscript will be available via the '*Peer-reviewed Publication DOI*' link on the right-hand side of this webpage. Please feel free to contact any of the authors; we welcome feedback.
¹ The Community Code Verification Exercise for

² Simulating Sequences of Earthquakes and Aseismic Slip
(SEAS)

³ Brittany A. Erickson¹*, Junle Jiang²*, Michael Barall³, Nadia Lapusta⁴,
Eric M. Dunham⁵, Ruth Harris⁶, Lauren S. Abrahams⁵, Kali L. Allison⁷,
Jean-Paul Ampuero⁸, Sylvain Barbot⁹, Camilla Cattania⁵, Ahmed Elbanna¹⁰,
Yuri Fialko¹¹, Benjamin Idini⁴, Jeremy E. Kozdon¹², Valère Lambert⁴,
Yajing Liu¹³, Yingdi Luo⁴, Xiao Ma¹⁰, Maricela Best McKay¹⁴,
Paul Segall⁵, Pengcheng Shi¹⁵, Martijn van den Ende⁸, Meng Wei¹⁵

¹University of Oregon, 1477 E. 13th Ave., Eugene, OR 97403-1202

²Cornell University

³Invisible Software

⁴California Institute of Technology

⁵Stanford University

⁶United States Geological Survey

⁷University of Maryland

⁸Université Côte d'Azur, IRD, CNRS,
Observatoire de la Côte d'Azur, Géoazur, France

⁹University of Southern California

¹⁰University of Illinois Urbana-Champaign

¹¹University of California San Diego

¹²Naval Postgraduate School

¹³McGill University

¹⁴Portland State University

¹⁵University of Rhode Island

5

September 4, 2019

<footnote>*Correspondence to B.A.E (bae@uoregon.edu) and J.J. (jjiang@cornell.edu)</footnote>
⁶ Any use of trade, firm, or product names is for descriptive purposes only and does not

⁷ imply endorsement by the U.S. Government.
Numerical simulations of Sequences of Earthquakes and Aseismic Slip (SEAS) have made great progress over the past decades to address important questions in earthquake physics and fault mechanics. However, significant challenges in SEAS modeling remain in resolving multiscale interactions between aseismic fault slip, earthquake nucleation, and dynamic rupture; and understanding physical factors controlling observables such as seismicity and ground deformation. The increasing capability and complexity of SEAS modeling calls for extensive efforts to verify codes and advance these simulations with rigor, reproducibility, and broadened impact. In 2018, we initiated a community code-verification exercise for SEAS simulations, supported by the Southern California Earthquake Center (SCEC). Here we report the findings from our first two benchmark problems (BP1 and BP2), designed to test the capabilities of different computational methods in correctly solving a mathematically well-defined, basic problem in crustal faulting. These benchmarks are for a 2D antiplane problem, with a 1D planar vertical strike-slip fault obeying rate-and-state friction, embedded in a 2D homogeneous, linear elastic half-space. Sequences of quasi-dynamic earthquakes with periodic occurrences (BP1) or bimodal sizes (BP2) and their interactions with aseismic slip are simulated. The comparison of >70 simulation results from 11 groups using different numerical methods, uploaded to our online platform, show excellent agreements in long-term and coseismic evolution of fault properties. In BP1, we found that the truncated domain boundaries influence interseismic fault stressing, earthquake recurrence, and coseismic rupture process, and that agreement between models is only achieved with sufficiently large domain sizes. In BP2, we found that complexity of long-term fault behavior depends on how well important physical length scales related to spontaneous nucleation and rupture propagation are resolved. Poor numerical resolution can result in the generation of artificial complexity, impacting simulation results that are of potential interest for characterizing seismic hazard, such as earthquake size distributions, moment release, and earthquake recurrence times. These results inform the development of more advanced SEAS models, contributing to our further understanding of earthquake system
dynamics.
# Introduction and Motivation

When we develop models of physical systems, credible and reproducible results are essential to scientific progress. Robust predictive models of earthquake source processes have become important means for studying fundamental questions in earthquake science. Models of single earthquakes (known as *dynamic rupture simulations*) have emerged as powerful tools for understanding the influence of fault geometry, friction and prestress on rupture propagation, and for explaining observations of high-frequency ground motions and damage zones (Day, 1982; Olsen et al., 1997; Nielsen et al., 2000; Duan and Oglesby, 2006; Ripperger et al., 2007; Bhat et al., 2007; Dunham et al., 2011a,b; Lozos et al., 2011; Gabriel et al., 2012; Shi and Day, 2013; Kozdon and Dunham, 2013; Xu et al., 2015; Wollherr et al., 2018; Ma and Elbanna, 2019). Many of the codes used for these studies incorporate advanced features such as 3D domains and complex fault geometries, leading to very large problems for which rigorous convergence tests can be too computationally expensive. An alternative means for verifying model results are code comparisons made across the different modeling groups, using cell sizes at the limit of computational feasibility. Over the past decade, the SCEC/USGS S spontaneous Rupture Code Verification Project has made significant progress in using code comparison studies to provide confidence in model outcomes (Harris et al., 2009; Barall and Harris, 2015; Harris et al., 2018).

Although these dynamic rupture simulations have contributed greatly to our understanding of the physical factors that govern ground motion, they are limited to single-event scenarios with imposed artificial prestress conditions and *ad hoc* nucleation procedures. In order to understand earthquake source processes and how fault slip history influences subsequent events, it has been widely recognized that we need models that simulate behavior over multiple seismic events and the intervening periods of aseismic deformation. To address this need, models of Sequences of Earthquakes and Aseismic Slip (SEAS) have emerged that consider all phases of earthquake faulting, from slow tectonic loading to earthquake nucleation (under self-consistent prestress
conditions), propagation and termination. However, so far codes for SEAS simulations remain untested. Inspired by the success of the SCEC/USGS Spontaneous Rupture Code Verification Project, this paper describes the efforts of the SEAS initiative – a SCEC (Southern California Earthquake Center) funded working group who has initiated the first code-verification study for earthquake sequence simulations. In this paper, we present the initial benchmark problems and results from the code comparisons submitted to our online platform (http://scecdata.usc.edu/cvws/seas/). Through these exercises, we aim to provide confidence in SEAS model outcomes, determine best practices for improvement of accuracy and efficiency of SEAS simulations, and provide other scientists strategies for verification during code development.

In SEAS models the goal is to capture the interplay of interseismic periods and the associated aseismic fault slip that ultimately lead to earthquake nucleation and earthquakes (dynamic rupture events) themselves, in an effort to understand which physical factors control the full range of observables such as aseismic deformation, nucleation locations of earthquakes, ground shaking during dynamic rupture, recurrence times and magnitudes of major earthquakes, see Figure 1. These features distinguish SEAS models from both dynamic rupture models which only consider single events, and the so-called earthquake simulators (Tullis et al., 2012). Earthquake simulators are capable of simulating seismicity patterns over millennium time scales in complex fault network systems (Richards-Dinger and Dieterich, 2012) but are missing key physical features that could potentially dominate earthquake and fault interaction, such as stress transfer generated by dynamic waves, aseismic slip within fault segments, and inelastic responses.

SEAS modeling is not without significant challenges, due to the varying temporal and spatial scales that characterize earthquake source behavior. For computational efficiency the vast majority of SEAS models do not consider full dynamics during earthquake rupture, but rather take a "quasi-dynamic" approach, where inertia is only approximated (see section for further details). Computations are further complicated when material heterogeneities, bulk inelastic responses and fault nonplanarity are in-
cluded. However, accounting for such complexity is widely recognized as crucial for understanding the real Earth and predicting seismic hazards. Significant developments in SEAS models over the past decade have incorporated some of these complexities and connected model outcomes to geophysical observations. For example, seismological and geodetic observations have been combined with modeling of coseismic and quasi-static (aseismic) deformation to infer the spatial distribution of fault frictional properties (Johnson et al., 2006; Barbot et al., 2009; Mitsui and Hirahara, 2011; Dublanchet et al., 2013; Floyd et al., 2016; Jiang and Fialko, 2016), the decay rate of aftershocks (Perfettini and Avouac, 2004, 2007), the role of tremor and slow slip (Mele Veedu and Barbot, 2016; Dublanchet, 2017; Luo and Ampuero, 2017), and long-term models have been used to reproduce characteristics of multiple and/or repeating events (Chen and Lapusta, 2009; Barbot et al., 2012). The framework of earthquake cycle modeling is also adopted to explain geodetic and geologic data (Meade et al., 2013; Kaneko et al., 2011; Wei et al., 2013, 2018), study subduction zones (Hori et al., 2004; van Dinther et al., 2013; Noda and Lapusta, 2013; Liu and Rice, 2005, 2007; Li and Liu, 2016, 2017), collision zones (Qiu et al., 2016; Michel et al., 2017), and explore induced seismicity phenomena (McClure and Horne, 2011; Dieterich et al., 2015), among many applications.

While SEAS models are being used to explain, reproduce, and predict earthquake behavior and other geophysical phenomena, a critical step must be to ensure that these methodologies are accurate. The SEAS initiative is also taking the step to improve and promote a new generation of verified numerical SEAS models that can simulate much longer periods of earthquake activity than single-event dynamic rupture simulations but with the same level of computational rigor, while incorporating qualitatively different features such as (a) pre-, inter-, and post-seismic slip and the resulting stress redistribution, (b) spontaneous earthquake nucleation, and (c) physical processes relevant to long-term slip such as interseismic healing of the fault zone, viscoelasticity, and fluid flow. Such SEAS models can provide physics-based approximations for larger-scale and longer-term earthquake simulators. In addition they can inform the initial condi-
tions and nucleation procedures for dynamic rupture simulations, however our vision for SEAS models is to develop them all to include full dynamic ruptures, capturing the range of processes and heterogeneities known to be essential for realistic ground motion modeling.

# SEAS Modeling Challenges and Initial Benchmark Problems

Although the ultimate SEAS modeling framework would naturally include dynamic rupture modeling, current methods for simulating SEAS problems require computational codes that are fundamentally different from those used in single-event dynamic rupture simulations. The use of variable time stepping and possible switching between different computational schemes is required in order to resolve sub-seconds to year-long changes. The interaction between the highly nonlinear nature of the problems and round-off errors can lead to model divergence. The need to distinguish between legitimate solution differences due and improper choices of algorithm and modeling procedures necessitates new and more suitable comparison metrics.

SEAS models are unique in that they cover a wide range of numerical methodologies and applications in earthquake science. Methods based on spectral boundary integral formulations (BIEM) are efficient in solving for earthquake ruptures with quasi-dynamic or full inertial effects (Lapusta and Rice, 2003; Lapusta and Liu, 2009; Jiang and Lapusta, 2016). Methods based on the finite difference method (FDM) or a hybrid finite element/spectral BIEM have been used to simulate quasi-dynamic ruptures on faults with more complex bulk rheologies (Erickson and Dunham, 2014; Erickson et al., 2017; Allison and Dunham, 2018; Mckay et al., 2019; Abdelmeguid et al., 2019). Other SEAS modeling approaches include boundary element methods (BEM) for simulating slow slip and tremor (e.g., Tse and Rice, 1986; Rice and Tse, 1986; Ong et al., 2019; Goswami and Barbot, 2018; Luo and Ampuero, 2011; Nakata et al., 2012; Liu, 2013; Wei et al., 2019).
2013), coupling faulting with fluid/heat transport and inelastic dilatancy (Segall and Bradley, 2012a), effects of surface topography (Ohtani and Hirahara, 2015), frictional heterogeneities (Kato, 2016) and viscoelastic response (Kato, 2002; Lambert and Barbot, 2016; Barbot, 2018). A spectral element method (SEM) has also been developed for simulating fully dynamic earthquakes in a heterogeneous bulk (Kaneko et al., 2010).

To verify the accuracy of SEAS models based on these different computational methods, the SEAS group developed our first benchmark problem, BP1, to test the capabilities of different computational methods in correctly solving a mathematically well-defined problem in crustal faulting. The overall strategy of our benchmark exercises is to produce robust results and maximize participation, with the goal of obtaining agreements in resolving detailed fault slip history over a range of time scales. These efforts required us to better understand the dependence of fault slip history on initial conditions, model spin-up, fault properties, and friction laws. Given the complexity of this task, it was important to start from the most basic problem and gradually add model complexity. BP1 is a 2D antiplane problem, with a 1D planar vertical strike-slip fault embedded in a 2D homogeneous, linear elastic half-space with a free surface, see Figure 2. Full details of this benchmark (and subsequent benchmarks), including governing equations and initial and fault interface conditions, are available online on the SEAS platform (http://scecdata.usc.edu/cvws/seas/index.html). We include some of the details on the friction law here, for clarity of important concepts.

The fault is governed by rate- and state-dependent friction (Dieterich, 1979; Ruina, 1983; Marone, 1998) where shear stress on the fault $\tau$ is set equal to fault strength $F$, namely

$$ \tau = F(V, \theta), \qquad (1) $$

where $\tau = \tau^0 + \tau^{qs} - \eta V$ is the sum of the prestress $\tau^0$, the shear stress due to quasi-static deformation $\tau^{qs}$, and the radiation damping term $-\eta V$ as approximation to inertia (Rice, 1993). $\eta = \mu/2c_s$ is half the shear-wave impedance for shear wave speed $c_s = \sqrt{\mu/\rho}$, where $\mu$ is the elastic shear modulus and $\rho$ is the material density. The
fault strength $F = \sigma_n f(V, \theta)$, where $V$ is the slip rate and $\theta$ is a state variable. $\sigma_n$ is the effective normal stress on the fault. For this first benchmark problem we assume $\theta$ evolves according to the aging law

$$ \frac{d\theta}{dt} = 1 - \frac{V\theta}{L}, \qquad (2) $$

where $L$ is the critical slip distance. The friction coefficient $f$ is given by a regularized formulation (Lapusta et al., 2000)

$$ f(V, \theta) = a \sinh^{-1} \left[ \frac{V}{2V_0} \exp \left( \frac{f_0 + b \ln(V_0\theta/L)}{a} \right) \right] \qquad (3) $$

for reference friction coefficient $f_0$ and reference slip rate $V_0$. Depth-dependent frictional parameters $a$ and $b$ define a shallow seismogenic region with velocity-weakening (VW) friction and a deeper velocity-strengthening (VS) region, below which a relative plate motion rate is imposed. A periodic sequence of spontaneous, quasi-dynamic earthquakes and slow slip are simulated in the model, see Figure 3a, where results from the BICyclE code (Lapusta et al., 2000; Lapusta and Liu, 2009) show slip contours plotted against fault depth in blue every 5 yr during interseismic loading and in red every 1 s during the coseismic phase. Over a 1200 year simulation period, approximately 13 events take place, nucleating at a depth of ~12 km, rupturing to a depth of ~18 km, and accumulating ~3 m of slip at the Earth's surface. Model parameters used for the benchmark are given in Table 1.

A critical physical length scale present in this first benchmark problem, often referred to as the process zone or cohesive zone $\Lambda$, describes the spatial region near the rupture front under which breakdown of fault resistance occurs, and shrinks as ruptures propagate faster (Palmer and Rice, 1973). For fault models governed by rate-and-state friction, the quasi-static process zone at a rupture speed of $0^+$, $\Lambda_0$, can be estimated (Day et al., 2005; Ampuero and Rubin, 2008; Perfettini and Ampuero, 2008) as

$$ \Lambda_0 = C \frac{\mu L}{b\sigma_n}, \qquad (4) $$
where C is a constant of order 1. Another characteristic length scale which has been shown to control model behavior is the critical nucleation size $h^*$, which governs the minimum extent of the rate-weakening region under which spontaneous nucleation may occur, (see Andrews, 1976a,b; Rubin and Ampuero, 2005; Ampuero and Rubin, 2008). For 2D problems, the critical nucleation size can be estimated for the aging law (with $0.5 < a/b < 1$) as

$$h^* = \frac{2 \mu b L}{\pi (b-a)^2 \sigma_n} \qquad (5)$$

A cell size of 50 m was used for BP1, resolving $\Lambda_0$ with approximately 6 grid points and $h^*$ with approximately 40 grid points.

We developed the second benchmark BP2 that is similar to BP1 to explore the model resolution issues, which will be important in future benchmarks in 3D when computational efficiency demands a larger cell size. Complexity of event sizes and recurrence times is known to emerge through a reduction in the characteristic slip distance $L$ (Lapusta and Rice, 2003; Mitsui and Hirahara, 2011; Wu and Chen, 2014; Kato, 2014; Barbot, 2019; Viesca, 2016a,b; Cattania, 2019). Thus BP2 is exactly the same as BP1 except that $L$ is halved, resulting in bimodal sequences of full and partial ruptures of the velocity-weakening region (every large event is accompanied by a smaller event and the sequence repeats periodically). Besides aiming for agreements between different models, one main objective is to understand complexity in simulated events and how to deal with numerical resolution issues. A reduction in $L$ corresponds to a reduction in the quasi-static process zone size $\Lambda_0$. BP2 requests model outputs using a cell size of 25 m, 50 m, 100 m, 200 m, 300 m, 400 m and 800 m. The first three cases resolve $\Lambda_0$ with approximately 6, 3, and 1.7 grid points, and the other four cases do not resolve $\Lambda_0$. Figures 3b-d show results from the BICycle code using a cell size of 25 m, 100 m and 200 m respectively. Small cell sizes of 25 m and 50 m (the latter is not shown) show nearly indistinguishable, bimodal patterns of events nucleating at ~15 km depth, suggesting model convergence. A cell size of 100 m leads to a resolution issue where periodic behavior is observed, but the bimodal sequence of events is replaced by
an alternating sequence of large, small and medium sized events. A cell size of 200 m,
which does not resolve the process zone, reveals a loss of periodic behavior altogether
in favor of a broad range of event sizes and nucleation locations.

Modeling Groups and Working Platforms

For these benchmark exercises, we have used two SCEC-funded workshops (hosted in April and November 2018, http://scecdata.usc.edu/cvws/seas/workshop_presentations.html) as open platforms for modelers to share and follow recent scientific progress in the field, discuss details in benchmark design/results, and collectively decide the directions of our future efforts, with considerable inputs from students and early career scientists.

Over 10 modeling groups participated in these first two benchmarks; the details of the group members and different computational methods are summarized in Table. Note that the modeler name refers to the member of the modeling group who uploaded the data to the platform for simulations done by the group. It does not necessarily refer to the code author(s) - see the references in Table for authorship and code availability. For time-stepping schemes, the majority of groups used adaptive Runge-Kutta methods for both benchmark problems (the details of which can be found in the references listed in Table ), with the exception of QDYN, which applies a Bulirsch-Stoer method for BP1, and BICyclE, which incorporates adaptive time-stepping based on stability conditions derived from the choice of constitutive relationship.

To facilitate the submission and comparison of simulation results, we established an online platform that provides access to community resources and supports the submission, storage, visualization, and comparison of benchmark results, see Figure 4. For our first benchmarks, we adopted a platform with similar functionality developed for the SCEC dynamic rupture simulation group (http://scecdata.usc.edu/cvws/seas/). All modelers can upload and immediately plot time-series data to quickly assess the overall agreements between models for the time evolution of fault slip, slip rates and shear stress at representative locations on fault. We use the online platform for prelim-
inary model comparisons and analyze more detailed model observables to verify these
computational codes.

Model Comparisons and What We Learned

It is important to note that the problem descriptions for BP1 and BP2 consider a semi-
infinite half-space. Codes based on a volume discretization (FDM/FEM) therefore had
to make their own decisions regarding computational domain truncation and far-field
boundary conditions. The figures in the following sections contain labels generated
by the platform which state the model group name and correspond to results from a
particular model set-up. Some results are followed by the version corresponding to an
alternative set-up, e.g. abrahams.3 corresponds to results from the abrahams group
with an increased computational domain size of (Lₓ, Lz) = (400 km, 200 km) and a
remote displacement boundary condition, see the lower right of Figure 4. We discuss
in the next sections the implications that these choices had on model comparisons.

Results from BP1

For the first benchmark problem, BP1, we found qualitative agreements in nucleation
sites, depth extent of rupture, and slip with depth similar to those exemplified by the
slip contours in Figure 3a. In Figure 5 we plot time series of local shear stress and
slip rates at mid-seismogenic depth (z = 7.5 km) from BP1 over the first 700 years for
different model results. Results from several BEM codes as well as codes with volume
discretization (abrahams and kozdon modeling groups) and varying computational
domain sizes are compared in Figure 5a-b. The legends indicate the computational
domain size and boundary condition. For BEM codes, HS refers to a half-space, and
(Lz, boundary condition) refers to computational domain depth and boundary condi-
tion, where BC3 corresponds to a periodic boundary condition. For codes with a vol-
ume discretization, (Lx/Lz/boundary condition), provides the computational domain
size used and BC1 and BC2 refers to a far-field free surface or a far-field displacement
284

boundary condition, respectively.

285

Figure 5(a-b) show model results from a BEM simulation (liu, in black) along with four model results from volume discretization codes, revealing quantitative differences in interevent times and peak values. Interevent times for different models range from approximately 78.3 to 78.8 years over the whole 3000 year simulation period, leading to model divergence at a near-constant rate. We found that these discrepancies were caused by choices in domain truncation and boundary conditions. We were surprised to find that far-field boundary condition type leads to quantitative differences in long-term fault behavior for relatively small domains (revealed by the blue and orange curves). This in part is due to small differences in the physical problem being solved by implementations that use periodic or finite domain boundary conditions compared to the spatial domain BEM methods which represent a truly infinite domain, and therefore larger loading regions. The green and red curves show how the discrepancy in long-term behavior among computational methodologies decreases as the physical domain size is increased, suggesting convergence of results across the modeling groups. Figure 5(c-d) shows comparisons of all models with $L_z > 160$ km, further illustrating that excellent agreements between model results can be achieved with sufficiently large domain sizes.

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

While computational domain size and boundary conditions can lead to model divergence over the long term, the coseismic behavior of individual earthquake are qualitatively well reproduced by all models. In Figure 6 we show the time series of shear stress evolution near the nucleation depth (12.5 km) and slip rate (at a mid-seismogenic depth of 7.5 km) during the coseismic phase for the 8th event in the sequence from Figure 5. We chose these plotting depths as they best illustrate model discrepancies, with time series aligned relative to the rupture initiation time at the depth of 12.5 km. Peak values in slip rates at 7.5 km depth occur approximately 10 s later, and co-seismic surface reflection phases are marked for all four plots with black arrows. Figures 6(a-b) show results from models on relatively small computational domains, revealing discrepancies in pre-rupture stress levels near the locked-creeping transition due to differences in interseismic loading, and resultant coseismic rupture behavior, including peak shear

302

303

304

305

306

307

308

309

310

311

312
stress and rupture speeds as evidenced by rupture initiation times of the direct and surface-reflection phases at depth of 7.5 km. Figures 6(c-d) illustrate excellent agreements for model results on larger domains. The discrepancy of <1 MPa in prestress levels at transitional depths does not result in pronounced difference in fault slip rate evolution.

## Results from BP2

For BP2 we suggested submissions of multiple models with different spatial resolutions from each group, see Table. By design, models with a cell size/node spacing that does not resolve critical length scales – process zone size and nucleation zone size defined in (4) and (5) – would produce increased complexity in earthquake sequences, observed previously (Rice, 1993; Ben-Zion and Rice, 1997; Day et al., 2005; Lapusta and Liu, 2009), and illustrated in the cumulative slip profiles in Figure 3(b-d).

While drastic differences in small event patterns arise for large cell sizes, we found that with increasing resolution results converge to an alternating sequence of large and small events among most models. Figure 7a shows the long term evolution of slip rates at 9.6 km (near the bottom of the seismogenic zone and above the earthquake initiation depth) for the best model results (with a cell size of 25 m and large computational domain sizes). We found that even models with similar cell/domain sizes tend to produce results that are initially closely matching, but diverge over time, likely due to accumulation of numerical round-off errors and differences in computational techniques. However, if we zoom in on the tenth event in the sequence (gray bar in Figure 7a), the time series of fault slip rates, aligned with respect to the start time of seismic slip at the depth of 12 km within each model, show good agreements (Figure 7b). While small discrepancies exist in peak slip rates and early source complexity, partly due to differences in interevent times, the models with the highest resolution exhibit good agreements in their overall coseismic behavior despite their divergence in the long term.

Figure 8 illustrates how model agreement is gradually lost with decreased model resolution. For cell sizes of 25 m and 50 m, long-term stress evolution near the locked-
341

342

343

344

345

346

creeping transition is qualitatively similar for the three models shown and the offset in the timing of earthquakes does not significantly affect coseismic behavior of major events, as indicated by comparable coseismic stress drops. For large cell sizes of 100 m and 200 m, not only is the time offset more random, but also coseismic stress drops and event patterns vary between models. Numerical artifacts and different computational techniques likely contribute to the divergence of simulation results.

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

In Figure 9 we plot the distribution of earthquake sizes, seismic moment release and frequency-size relation for two groups of models (jiang and cattania) with increasing cell sizes. For the 2D problem, we define earthquake size as moment release per length for each event, $M = \int \mu s \, dz$, where shear modulus $\mu = c_s^2 \rho \approx 32$ GPa and $s$ is total coseismic slip over the cell. While better resolved models (cell sizes of 25 m and 50 m) show excellent agreements between the two groups, models produce dramatically different earthquake statistics when cell size increases to 400 m, with the most significant discrepancies in smaller earthquakes between the two models (Figure 9a). The distribution of total seismic moment release, $M_t$, calculated as the sum of moment release during all earthquakes within a certain magnitude range, also changes with cell sizes, though in a similar manner for the two model groups (Figure 9b). Overall, models with larger cell sizes tend to produce large earthquakes with reduced total moment; part of the moment deficit is accommodated through many more smaller earthquakes and the rest through additional aseismic slip. For example, the total moment release through largest earthquakes in 400-m models is only half of that in 25-m models. These results demonstrate that simulated small earthquakes are especially sensitive to model resolution and large earthquake behavior can also be affected. In addition, Figure 9c reveals how different simulations with poor resolution can produce similar power-law features in frequency-size distributions over certain ranges of earthquake sizes, as a result of numerical artifacts rather than well-resolved physics.

367

368

369

In Figure 10 we illustrate the effect of model resolution on the partition between seismic and aseismic slip. Normalized seismic moment release $R_s$ is plotted against depth for several modeling groups, in solid lines for total seismic moment release and
370

371

372

373

374

375

376

377

dashed lines for seismic moment due to surface-breaching events. $R_s = 0$ implies that all plate motion is accommodated by aseismic slip on the fault, while $R_s = 1$ means that all moment is released through earthquakes. A transitional zone in this partitioning around $z = H$ and down to $z = H+h$ ($H = 15$ and $h = 3$ in this exercise) is evidenced in the well-resolved models ($\Delta z = 25, 50, 100$ m). The poorly resolved models, however, illustrate model discrepancies in the seismic/aseismic partitioning, with the near-surface slip budget being increasingly accommodated by small earthquakes and aseismic slip with increased cell sizes.

378

379

380

381

382

383

384

385

386

In Figure 11 we show interevent times for large surface-breaching events for all models and cell sizes, showing a strong agreement of ~110 years for a cell size of 25 m, with an increasing variability and discrepancies among models with increased cell size. Although the range of earthquake recurrence intervals are highly dependent on cell sizes, the median values across models with larger cell sizes do not significantly deviate from the uniform recurrence intervals in well-resolved models. This suggests that at least some observables in these models retain information of the true behavior of physical models and the larger cell sizes can be viewed as a factor that leads to increased modeling errors.

387

# Conclusions and Perspectives

388

389

390

391

392

393

394

395

396

For the first two SEAS benchmarks we found that discrepancies among well-resolved models were significantly influenced by computational domain size, with larger domains yielding improvements in agreements, regardless of domain boundary conditions. Spin-up periods (time required for system to be independent of initial conditions) for well-resolved models was relatively short - approximately 2-3 events. Results on large domains agree well initially but still diverge over time, which was not unexpected due to accumulation of round-off errors and differences in computational techniques. For BP2 we investigated model resolution and observed qualitative similarities of bimodal events when the process zone was resolved by approximately 3 and 6 grid points, sug-
gesting model convergence. A failure to resolve this length scale however, can lead to substantial differences in long-term fault behavior as well as earthquake statistics relevant to seismic hazard, such as frequency-size distributions and interevent times.

Although our initial benchmarks have a simple setup, comparison of results for tens of models have yielded some unexpected and important insights, affirming the importance of starting simple in a community code verification exercise. The results and lessons from our initial benchmarks prepare us for future benchmark problems that incrementally incorporate additional, potentially dominating physical factors, including fully dynamic ruptures, coupling with fluids, multiple fault segments, nonplanar fault geometries, and inelastic bulk constitutive behavior (e.g., *Segall and Rice*, 1995; *Noda and Lapusta*, 2010; *Segall and Rice*, 2006; *Segall et al.*, 2010; *Erickson et al.*, 2017; *Lambert and Barbot*, 2016; *Qiu et al.*, 2016; *Barbot*, 2018; *Ong et al.*, 2019). For future verification exercises, we plan to address important issues in SEAS simulations, such as 3D effects, heterogeneous fault frictional properties, and full dynamics, which should advance the state-of-the-art computational capabilities in our field.

The goal of the SEAS initiative is to promote advanced models with robust physical features—a large spectrum of rupture styles and patterns, including slow-slip events, complex earthquake sequences, fluid effects, dynamic stress changes, and inelastic deformation—that are currently missing in the large-scale, long-term earthquake simulator frameworks such as ViscoSim, RSQSIM, Virtual California, and ALLCAL (*Pollitz*, 2012; *Richards-Dinger and Dieterich*, 2012; *Sachs et al.*, 2012; *Ward*, 2012). This new generation of verified SEAS models will help determine the controlling physical mechanisms of earthquake nucleation, propagation, and arrest. The community-wide initiative would also provide incentives and new ideas to characterize modeling uncertainty for the increasingly complex earthquake source models, an important step in using physics-based models for the assessment of seismic hazard. Future validation efforts comparing physics-based models with geophysical observations will bridge studies in paleoseismology, geodesy, and seismology to understand fault behavior over multiple temporal and spatial scales.
426

427

428

**Data and Resources:** Our online platform (http://scecdata.usc.edu/cvws/seas/) is being developed and maintained by M.B. The data for local fault properties are stored on the platform.

429

430

431

432

433

434

435

**Author Contributions:** B.A.E. and J.J. designed the benchmark problems, analyzed model results, co-organized the workshops and co-wrote this article. M.B. developed and maintains the online platform. M.B., N.L., E.M.D. and R.H. provided major support and advice in forming the working group, obtaining funding, and manuscript writing. Remaining co-authors provided feedback on benchmark design, participated in the benchmark exercises, helped revise the manuscript, and are listed alphabetically.

436

437

438

439

440

**Acknowledgments:** B.A.E., J.J. and M.B. were supported through the Southern California Earthquake Center, grant no. 18099 and 19109. Two SEAS-themed workshops were funded by SCEC awards no. 17151 and 18102. SCEC is funded by NSF Cooperative Agreement EAR-0529922 and USGS Cooperative Agreement 07HQAG0008. This is SCEC contribution no. 9066.
Abdelmeguid, M., X. Ma, and A. Elbanna, A novel hybrid finite element-spectral boundary integral scheme for modeling earthquake cycles: Application to rate and state faults with low-velocity zones, *J. Geophys. Res. (to appear)*, doi:10.31223/osf.io/xwhbs, 2019.

Allison, K., and E. M. Dunham, Earthquake cycle simulations with rate-and-state friction and power-law viscoelasticity, *Tectonophysics*, 733, 232 – 256, doi:https://doi.org/10.1016/j.tecto.2017.10.021, 2018.

Ampuero, J.-P., and A. M. Rubin, Earthquake nucleation on rate and state faults – Aging and slip laws, *J. Geophys. Res.-Sol. Ea.*, 113 (B1), doi:10.1029/2007JB005082, 2008.

Andrews, D. J., Rupture propagation with finite stress in antiplane strain, *J. Geophys. Res.*, 81, 3575–3582, 1976a.

Andrews, D. J., Rupture velocity of plane strain shear cracks, *J. Geophys. Res.*, 81, 5679–5687, 1976b.

Barall, M., and R. A. Harris, Metrics for comparing dynamic earthquake rupture simulations, *Seismol. Res. Lett.*, 86(1), doi:10.1785/0220140122, 2015.

Barbot, S., Asthenosphere flow modulated by megathrust earthquake cycles, *Geophys. Res. Lett.*, 45, 6018–6031, doi:10.1029/2018GL078197, 2018.

Barbot, S., Slow-slip, slow earthquakes, period-two cycles, full and partial ruptures, and deterministic chaos in a single asperity fault, *Tectonophysics*, doi:10.1016/j.tecto.2019.228171, 2019.

Barbot, S., Y. Fialko, and Y. Bock, Postseismic deformation due to the Mw 6.0 2004 Parkfield earthquake: Stress-driven creep on a fault with spatially variable rate-and-state friction parameters, *J. Geophys. Res.*, 114 (B07405), doi:10.1029/2008JB005748, 2009.
467

Barbot, S., N. Lapusta, and J.-P. Avouac, Under the hood of the earthquake machine: Toward predictive modeling of the seismic cycle, *Science*, **336**(6082), 707-710, doi: 10.1126/science.1218796, 2012.

468

469

470

Ben-Zion, Y., and J. R. Rice, Dynamic simulations of slip on a smooth fault in an elastic solid, *J. Geophys. Res.-Sol. Ea.*, **102**(B8), 17,771-17,784, doi:10.1029/97JB01341, 1997.

471

472

473

Bhat, H. S., M. Olives, R. Dmowska, and J. R. Rice, Role of fault branches in earthquake rupture dynamics, *J. Geophys. Res.*, **112**(B11309), 16, doi:10.1029/2007JB005027, 2007.

474

475

476

Bradley, A. M., Software for efficient static dislocation-traction calculations in fault simulators, *Seismol. Res. Lett.*, **85**(6), 1358-1365, doi:10.1785/0220140092, 2014.

477

478

Cattania, C., Complex earthquake behavior on simple faults, submitted to *Geophys. Res. Lett.*, doi:10.31223/osf.io/hgbjx, 2019.

479

480

Chen, T., and N. Lapusta, Scaling of small repeating earthquakes explained by interaction of seismic and aseismic slip in a rate and state fault model, *J. Geophys. Res.*, **114**(B01311), doi:10.1029/2008JB005749, 2009.

481

482

483

Day, S. M., Three-dimensional finite difference simulation of fault dynamics: rectangular faults with fixed rupture velocity, *Bull. Seismol. Soc. Am.*, **72**(3), 705-727, 1982.

484

485

486

Day, S. M., L. A. Dalguer, N. Lapusta, and Y. Liu, Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture, *J. Geophys. Res.-Sol. Ea.*, **110**(B12), doi:10.1029/2005JB003813, 2005.

487

488

489

Dieterich, J. H., Modeling of rock friction: 1. experimental results and constitutive equations, *J. Geophys. Res.-Sol. Ea.*, **84**(B5), 2161-2168, doi:10.1029/JB084iB05p02161, 1979.

490

491
492

Dieterich, J. H., K. B. Richards-Dinger, and K. A. Kroll, Modeling injection-induced seismicity with the physics-based earthquake simulator RSQSim, *Seismol. Res. Lett.*, 86(4), doi:10.1785/0220150057, 2015.

493

494

495

Duan, B., and D. D. Oglesby, Heterogeneous fault stresses from previous earthquakes and the effect on dynamics of parallel strike-slip faults, *J. Geophys. Res.*, 111 (B5), 2006.

496

497

498

Dublanchet, P., The dynamics of earthquake precursors controlled by effective friction, *Geophys. J. Int.*, 212(2), 853-871, 2017.

499

500

Dublanchet, P., P. Bernard, and P. Favreau, Interactions and triggering in a 3-d rate-and-state asperity model, *J. Geophys. Res.*, 118(5), 2225-2245, 2013.

501

502

Dunham, E. M., D. Belanger, L. Cong, and J. E. Kozdon, Earthquake ruptures with strongly rate-weakening friction and off-fault plasticity, Part:1: Planar faults, *Bull. Seismol. Soc. Am.*, 101(5), 2296-2307, doi:10.1785/0120100075, 2011a.

503

504

505

Dunham, E. M., 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, doi:10.1785/0120100076, 2011b.

506

507

508

Erickson, B. A., 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.-Sol. Ea.*, 119(4), 3290-3316, doi:10.1002/2013JB010614, 2014.

509

510

511

512

Erickson, B. A., E. M. Dunham, and A. Khosravifar, A finite difference method for off-fault plasticity throughout the earthquake cycle, *J. Mech. Phys. Solids*, 109, 50-77, doi:10.1016/j.jmps.2017.08.002, 2017.

513

514

515

Floyd, M. A., et al., Spatial variations in fault friction related to lithology from rupture and afterslip of the 2014 South Napa, California, earthquake, *Geophys. Res. Lett.*, 43(13), 6808-6816, doi:10.1002/2016GL069428, 2016.

516

517
518

Gabriel, A.-A., J.-P. Ampuero, L. A. Dalguer, and P. M. Mai, The transition of dynamic rupture styles in elastic media under velocity-weakening friction, *J. Geophys. Res.*, **117**(B9), 2012.

519

520

521

Goswami, A., and S. Barbot, Slow-slip events in semi-brittle serpentinite fault zones, *Scientific reports*, **8**(1), 6181, 2018.

522

523

Harris, R. A., et al., The SCEC/USGS dynamic earthquake rupture code verification exercise, *Seismol. Res. Lett.*, **80**, 119–126, doi:10.1785/gssrl.80.1.119, 2009.

524

525

Harris, R. A., et al., A suite of exercises for verifying dynamic earthquake rupture codes, *Seismol. Res. Lett.*, **89**(3), doi:https://doi.org/10.1785/0220170222, 2018.

526

527

Hori, T., N. Kato, K. Hirahara, T. Baba, and Y. Kaneda, A numerical simulation of earthquake cycles along the Nankai Trough in southwest Japan: Lateral variation in frictional property due to the slab geometry controls the nucleation position, *Earth Planet Sc. Lett.*, **228**(3-4), 215–226, 2004.

528

529

530

531

Jiang, J., and Y. Fialko, Reconciling seismicity and geodetic locking depths on the Anza section of the San Jacinto fault, *Geophys. Res. Lett.*, **43**(20), 10,663–10,671, doi:10.1002/2016GL071113, 2016.

532

533

534

Jiang, J., and N. Lapusta, Deeper penetration of large earthquakes on seismically quiescent faults, *Science*, **352**(6291), 1293–1297, doi:10.1126/science.aaf1496, 2016.

535

536

Johnson, K. M., R. Bürgmann, and K. Larson, Frictional properties of the San Andreas Fault near Parkfield, California, inferred from models of afterslip following the 2004 earthquake, *Bull. Seismol. Soc. Am.*, **96**(4B), doi:10.1785/0120050808, 2006.

537

538

539

Kaneko, Y., J.-P. Avouac, and N. Lapusta, Towards inferring earthquake patterns from geodetic observations of interseismic coupling, *Nat. Geosci.*, **3**, 2010.

540

541

Kaneko, Y., J.-P. Ampuero, and N. Lapusta, Spectral-element simulations of long-term fault slip: Effect of low-rigidity layers on earthquake-cycle dynamics, *J. Geophys. Res.-Sol. Ea.*, **116**(B10), doi:10.1029/2011JB008395, 2011.

542

543
544

Kato, N., Seismic cycle on a strike-slip fault with rate-and state-dependent strength in an elastic layer overlying a viscoelastic half-space, *Earth Planets Space*, 54(11), 1077-1083, 2002.

545

546

547

Kato, N., Deterministic chaos in a simulated sequence of slip events on a single isolated asperity, *Geophys. J. Int.*, 198(2), 727-736, 2014.

548

549

Kato, N., Earthquake cycles in a model of interacting fault patches: Complex behavior at transition from seismic to aseismic slip, *Bull. Seismol. Soc. Am.*, 106(4), 1772-1787, 2016.

550

551

552

Kozdon, J. E., and E. M. Dunham, Rupture to the trench: Dynamic rupture simulations of the 11 March 2011 Tohoku earthquake, *Bull. Seismol. Soc. Am.*, 103(2B), 1275-1289, doi:10.1785/0120120136, 2013.

553

554

555

Lambert, V., and S. Barbot, Contribution of viscoelastic flow in earthquake cycles within the lithosphere-asthenosphere system, *Geophys. Res. Lett.*, 43(19), 142-154, 2016.

556

557

558

Lapusta, N., and Y. Liu, Three-dimensional boundary integral modeling of spontaneous earthquake sequences and aseismic slip, *J. Geophys. Res.-Sol. Ea.*, 114 (B9), doi: 10.1029/2008JB005934, 2009.

559

560

561

Lapusta, N., and J. R. Rice, Nucleation and early seismic propagation of small and large events in a crustal earthquake model, *J. Geophys. Res.*, 108 (B4, 2205), 2003.

562

563

Lapusta, N., 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.-Sol. Ea.*, 105 (B10), 23,765-23,789, doi:10.1029/2000JB900250, 2000.

564

565

566

567

Li, D., and Y. Liu, Spatiotemporal evolution of slow slip events in a nonplanar fault model for northern Cascadia subduction zone, *J. Geophys. Res.*, 121(9), 6828-6845, 2016.

568

569
570

Li, D., and Y. Liu, Modeling slow-slip segmentation in Cascadia subduction zone con-
strained by tremor locations and gravity anomalies, *J. Geophys. Res.*, **122**(4), 3138-
3157, 2017.

571

572

573

Liu, Y., Numerical simulations on megathrust rupture stabilized under strong dilatancy strengthening in slow slip region, *Geophys. Res. Lett.*, **40**(7), 1311–1316, doi:10.1002/grl.50298, 2013.

574

575

576

Liu, Y., and J. R. Rice, Aseismic slip transients emerge spontaneously in three-
dimensional rate and state modeling of subduction earthquake sequences, *J. Geophys.
Res.*, **110**(B08307), doi:10.1029/2004JB003424, 2005.

577

578

579

Liu, Y., and J. R. Rice, Spontaneous and triggered aseismic deformation transients in a
subduction fault model, *J. Geophys. Res.*, **112**(B09404), doi:10.1029/2007JB004930,
2007.

580

581

582

Lozos, J. C., D. D. Oglesby, B. Duan, and S. G. Wesnousky, The effects of double fault bends on rupture propagation: A geometrical parameter study, *Bull. Seismol. Soc. Am.*, **101**(1), 385–398, 2011.

583

584

585

Luo, Y., and J. P. Ampuero, Numerical simulation of tremor migration triggered by
slow slip and rapid tremor eversals, *AGU Fall Meeting Abstracts*, S33C-02, 2011.

586

587

Luo, Y., and J. P. Ampuero, Preprint: Tremor migration patterns and the collective
behavior of deep asperities mediated by creep, doi:10.31223/osf.io/mbcav, 2017.

588

589

Ma, X., and A. Elbanna, Dynamic rupture propagation on fault planes with explicit
representation of short branches, *Earth Planet Sc. Lett.*, **523**, 115,702, doi:10.1016/j.
epsl.2019.07.005, 2019.

590

591

592

Marone, C., Laboratory-derived friction laws and their application to seismic faulting,
*Ann. Rev. Earth Pl. Sc.*, **26**(1), 643–696, doi:10.1146/annurev.earth.26.1.643, 1998.

593
594

McClure, M. W., and R. N. Horne, Investigation of injection-induced seismicity using a coupled fluid flow and rate/state friction model, *Geophysics*, 76(6), WC181-WC198, doi:10.1190/geo2011-0064.1, 2011.

595

596

597

Mckay, M. B., B. A. Erickson, and J. E. Kozdon, A computational method for earthquake cycles within anisotropic media, *Geophys. J. Int.*, doi:10.1093/gji/ggz320, 2019.

598

599

600

Meade, B. J., Y. Klinger, and E. A. Hetland, Inference of multiple earthquake-cycle relaxation timescales from irregular geodetic sampling of interseismic deformation, *Bull. Seismol. Soc. Am.*, 103(5), 2824-2835, doi:10.1785/0120130006, 2013.

601

602

603

Mele Veedu, M., and S. Barbot, The Parkfield tremors reveal slow and fast ruptures on the same asperity, *Nature*, 532(7599), 361-365, doi:10.1038/nature17190, 2016.

604

605

Michel, S., J.-P. Avouac, N. Lapusta, and J. Jiang, Pulse-like partial ruptures and high-frequency radiation at creeping-locked transition during megathrust earthquakes, *Geophys. Res. Lett.*, 44(16), 8345-8351, doi:10.1002/2017GL074725, 2017.

606

607

608

Mitsui, Y., and K. Hirahara, Fault instability on a finite and planar fault related to early phase of nucleation, *J. Geophys. Res.*, 116(B6), 2011.

609

610

Nakata, R., M. Hyodo, and T. Hori, Numerical simulation of afterslips and slow slip events that occurred in the same area in Hyuga-nada of southwest Japan, *Geophys. J. Int.*, 190(2), 1213-1220, doi:10.1111/j.1365-246X.2012.05552.x, 2012.

611

612

613

Nielsen, S. B., J. Carlson, and K. B. Olsen, Influence of friction and fault geometry on earthquake rupture, *J. Geophys. Res.*, 105(B3), 6069-6088, 2000.

614

615

Noda, H., and N. Lapusta, Three-dimensional earthquake sequence simulations with evolving temperature and pore pressure due to shear heating: Effect of heterogeneous hydraulic diffusivity, *J. Geophys. Res.*, 115(B12314), doi:10.1029/2010JB007780, 2010.

616

617

618
619

Noda, H., and N. Lapusta, Stable creeping fault segments can become destructive as a result of dynamic weakening, *Nature*, **493** (7433), 518–521, 2013.

620

621

Ohtani, M., and K. Hirahara, Effect of the Earth's surface topography on quasi-dynamic earthquake cycles, *Geophys. J. Int.*, **203**(1), 384–398, doi:10.1093/gji/ggv187, 2015.

622

623

Olsen, K. B., R. Madariaga, and R. J. Archuleta, Three-dimensional dynamic simulation of the 1992 Landers earthquake, *Science*, **278**(5339), 834–838, 1997.

624

625

Ong, M., Su Qing, S. Barbot, and J. Hubbard, Physics-based scenario of earthquake cycles on the Ventura thrust system, California: The effect of variable friction and fault geometry, *J. Pure Appl. Geophys.*, doi:10.1007/s00024-019-02111-9, 2019.

626

627

628

Palmer, A. C., and J. R. Rice, The growth of slip surfaces in the progressive failure of over-consolidated clay, *Proc. Roy. Soc. Lond. A*, **332**(1591), 527–548, doi:10.1098/rspa.1973.0040, 1973.

629

630

631

Perfettini, H., and J.-P. Ampuero, Dynamics of a velocity strengthening fault region: Implications for slow earthquakes and postseismic slip, *J. Geophys. Res.-Sol. Ea.*, **113**(B9), doi:10.1029/2007JB005398, 2008.

632

633

634

Perfettini, H., and J.-P. Avouac, Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the Chi-Chi earthquake, Taiwan, *J. Geophys. Res.*, **109**(B02304), 2004.

635

636

637

Perfettini, H., and J.-P. Avouac, Modeling afterslip and aftershocks following the 1992 Landers earthquake, *J. Geophys. Res.*, **112**(B07409), 2007.

638

639

Pollitz, F. F., ViscoSim Earthquake Simulator, *Seismol. Res. Lett.*, **83**(6), 979–982, doi:10.1785/0220120050, 2012.

640

641

Qiu, Q., et al., The mechanism of partial rupture of a locked megathrust: The role of fault morphology, *Geology*, **44**(10), 875–878, doi:10.1130/G38178.1, 2016.

642
643

Rice, J. R., Spatio-temporal complexity of slip on a fault, *J. Geophys. Res.*, **98**(B6), 9885–9907, 1993.

644

645

Rice, J. R., and S. T. Tse, Dynamic motion of a single degree of freedom system following a rate and state dependent friction law, *J. Geophys. Res.*, **91**(B1), 521–530, 1986.

646

647

Richards-Dinger, K., and J. H. Dieterich, RSQSim earthquake simulator, *Bull. Seismol. Soc. Am.*, **83**(6), 983–990, doi:https://doi.org/10.1785/0220120105, 2012.

648

649

Ripperger, J., J.-P. Ampuero, P. M. Mai, and D. Giardini, Earthquake source characteristics from dynamic rupture with constrained stochastic fault stress, *J. Geophys. Res.-Sol. Ea.*, **112**(B4), doi:10.1029/2006JB004515, 2007.

650

651

652

Rubin, A. M., and J.-P. Ampuero, Earthquake nucleation on (aging) rate and state faults, *J. Geophys. Res.-Sol. Ea.*, **110**(B11), doi:10.1029/2005JB003686, 2005.

653

654

Ruina, A., Slip instability and state variable friction laws, *J. Geophys. Res.-Sol. Ea.*, **88**(B12), 10,359–10,370, doi:10.1029/JB088iB12p10359, 1983.

655

656

Sachs, M. K., E. M. Heien, D. L. Turcotte, M. B. Yikilmaz, J. B. Rundle, and L. H. Kellogg, Virtual California Earthquake Simulator, *Seismol. Res. Lett.*, **83**(6), 973–978, doi:10.1785/0220120052, 2012.

657

658

659

Segall, P., and A. M. Bradley, The role of thermal pressurization and dilatancy in controlling the rate of fault slip, *J. Appl. Mech.*, **79**(3), doi:10.1115/1.4005896, 2012a.

660

661

Segall, P., and A. M. Bradley, Slow-slip evolves into megathrust earthquakes in 2d numerical simulations, *Geophys. Res. Lett.*, **39**(18), doi:10.1029/2012GL052811, 2012b.

662

663

Segall, P., and J. R. Rice, Dilatancy, compaction, and slip instability of a fluid infiltrated fault, *J. Geophys. Res.*, **100**, 22,155–22,171, 1995.

664

665

Segall, P., and J. R. Rice, Does shear heating of pore fluid contribute to earthquake nucleation?, *J. Geophys. Res.*, **111**(B09316), 17, 2006.

666
667

Segall, P., A. M. Rubin, A. M. Bradley, and J. R. Rice, Dilatant strengthening as a mechanism for slow slip events, *J. Geophys. Res.*, **115** (B12305), 2010.

668

669

Shi, Z., and S. M. Day, Rupture dynamics and ground motion from 3-d rough-fault simulations, *J. Geophys. Res.-Sol. Ea.*, **118**(3), 1122–1141, doi:10.1002/jgrb.50094, 2013.

670

671

672

Tse, S. T., and J. R. Rice, Crustal earthquake instability in relation to the depth variation of frictional slip properties, *J. Geophys. Res.*, **91**(B9), 9452–9472, 1986.

673

674

Tullis, T. E., et al., Generic earthquake simulator, *Seismol. Res. Lett.*, **83**, 959–963, doi:10.1785/0220120093, 2012.

675

676

van Dinther, Y., T. V. Gerya, L. A. Dalguer, F. Corbi, F. Funiciello, and P. M. Mai, The seismic cycle at subduction thrusts: 2. dynamic implications of geodynamic simulations validated with laboratory models, *J. Geophys. Res.-Sol. Ea.*, **118**(4), 1502–1525, doi:10.1029/2012JB009479, 2013.

677

678

679

680

Viesca, R. C., Stable and unstable development of an interfacial sliding instability, *Phys. Rev. E*, **93**, 060202, doi:10.1103/PhysRevE.93.060202, 2016a.

681

682

Viesca, R. C., Self-similar slip instability on interfaces with rate- and state-dependent friction, *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, **472**(2192), 20160,254, doi:10.1098/rspa.2016.0254, 2016b.

683

684

685

Ward, S. N., ALLCAL Earthquake Simulator, *Seismol. Res. Lett.*, **83**(6), 964–972, doi:10.1785/0220120056, 2012.

686

687

Wei, M., Y. Kaneko, Y. Liu, and J. J. McGuire, Episodic fault creep events in California controlled by shallow frictional heterogeneity, *Nat. Geosci.*, **6**, 566–570, doi:10.1038/ ngeo1835, 2013.

688

689

690

Wei, M., Y. Kaneko, P. Shi, and Y. Liu, Numerical modeling of dynamically triggered shallow slow slip events in New Zealand by the 2016 Mw 7.8 Kaikoura earthquake, *Geophys. Res. Lett.*, **45**(10), 4764–4772, 2018.

691

692
693

Wollherr, S., A.-A. Gabriel, and P. M. Mai, Landers 1992 "reloaded": Integrative dynamic earthquake rupture modeling, *J. Geophys. Res.-Sol. Ea.*, 2018.

694

695

Wu, Y., and X. Chen, The scale-dependent slip pattern for a uniform fault model obeying the rate-and state-dependent friction law, *J. Geophys. Res.*, 119(6), 4890-4906, 2014.

696

697

698

Xu, J., H. Zhang, and X. Chen, Rupture phase diagrams for a planar fault in 3-d full-space and half-space, *Geophys. J. Int.*, 202(3), 2194-2206, 2015.

699
<table><caption>Table 1: Parameter values used in the benchmark problem</caption><thead><tr><th>Parameter</th><th>Definition</th><th>Value, Units</th></tr></thead><tbody><tr><td>ρ</td><td>density</td><td>2670 kg/m³</td></tr><tr><td>c<sub>s</sub></td><td>shear wave speed</td><td>3.464 km/s</td></tr><tr><td>σ<sub>n</sub></td><td>effective normal stress on fault</td><td>50 MPa</td></tr><tr><td>a</td><td>rate-and-state parameter</td><td>variable (see Fig. 1)</td></tr><tr><td>b</td><td>rate-and-state parameter</td><td>variable (see Fig. 1)</td></tr><tr><td>L</td><td>critical slip distance</td><td>BP1: 0.008 m<br>BP2: 0.004 m</td></tr><tr><td>V<sub>p</sub></td><td>plate rate</td><td>10<sup>-9</sup> m/s</td></tr><tr><td>V<sub>init</sub></td><td>initial slip rate</td><td>10<sup>-9</sup> m/s</td></tr><tr><td>V<sub>0</sub></td><td>reference slip rate</td><td>10<sup>-6</sup> m/s</td></tr><tr><td>f<sub>0</sub></td><td>reference friction coefficient</td><td>0.6</td></tr><tr><td>H</td><td>depth extent of uniform VW region</td><td>15 km</td></tr><tr><td>h</td><td>width of VW-VS transition zone</td><td>3 km</td></tr><tr><td>W<sub>f</sub></td><td>width of rate-and-state fault</td><td>40 km</td></tr><tr><td>Δz</td><td>suggested cell sizes</td><td>BP1: 25 m<br>BP2: 25 m, 50 m, 100 m, 200 m, 300 m, 400 m, 800 m</td></tr><tr><td>t<sub>f</sub></td><td>final simulation time</td><td>BP1: 3000 years<br>BP2: 1200 years</td></tr><tr><td>L<sub>z</sub></td><td>depth of computational domain</td><td>not specified</td></tr><tr><td>L<sub>x</sub></td><td>off-fault distance of computational domain</td><td>not specified</td></tr></tbody></table>
<table><caption>Table 2: Details of participating SEAS codes and modeling groups.</caption><thead><tr><td>Code Name</td><td>Type</td><td>Modeler Name &amp; Group Members</td><td>References</td></tr></thead><tbody><tr><td>SCycle</td><td>FDM</td><td>abrahams (Abrahams/ Allison/Dunham)</td><td>Erickson and Dunham (2014)<br/>Allison and Dunham (2018)<br/>https://github.com/kali-allison/SCycle</td></tr><tr><td>FDCycle</td><td>FDM</td><td>erickson (Erickson/Mckay)</td><td>Erickson and Dunham (2014)<br/>https://github.com/brittany-erickson/FDCycle</td></tr><tr><td>QDESDG</td><td>DG-FEM</td><td>kozdon (Kozdon)</td><td>https://github.com/jkozdon/QDESDG</td></tr><tr><td>Unicycle</td><td>BEM</td><td>barbot (Barbot)</td><td>Barbot (2019)<br/>http://bitbucket.org/sbarbot</td></tr><tr><td>FDRA</td><td>BEM</td><td>cattania (Cattania/Segall)</td><td>Segall and Bradley (2012b); Bradley (2014)</td></tr><tr><td>BICyclE</td><td>BEM</td><td>jiang (Jiang)<br/>lambert (Lambert/Lapusta)<br/>xma (Ma/Elbanna)</td><td>Lapusta et al. (2000); Lapusta and Liu (2009)</td></tr><tr><td>QDYN</td><td>BEM</td><td>luo (Luo/Idini/ van den Ende/Ampuero)</td><td>Luo and Ampuero (2017)<br/>https://github.com/ydluo/qdyn</td></tr><tr><td>ESAM</td><td>BEM</td><td>liu (Liu)<br/>wei (Wei/Shi)</td><td>Liu and Rice (2007)</td></tr></tbody></table>
<figure><img src="image_2.png" /><figcaption>Figure 1: Ingredients and observables for SEAS (sequences of earthquakes and aseismic slip) models. In a conceptual fault-zone model, earthquakes initiate at seismogenic depths (red star) and rupture through the interseismically locked regions (gray), while aseismic slip occurs in deeper and sometimes shallower regions (yellow). For numerical models, given fault zone properties, computational simulations can reproduce long-term fault locking and creep over years to decades, punctuated by dynamic earthquake ruptures over seconds to minutes. Seismic shaking and aseismic deformation are typical observables from the surface.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Our first SEAS benchmark is based on the model in *Rice* (1993), where a planar fault is embedded in a homogeneous, linear elastic half-space with a free surface. A vertical cross-section of the 3D setting is taken so that slip varies only with depth and deformation is 2D antiplane strain. The fault is governed by rate-and-state friction with depth-dependent frictional parameters *a* and *b* above the depth *W*<sub>f</sub>, below which a steady slow loading rate *V*<sub>p</sub> is assumed. The friction-controlled fault is seismogenic due to velocity-weakening properties ((*a* − *b*) < 0) down to depth *H* and accommodates aseismic creep at greater depths due to velocity-strengthening properties ((*a* − *b*) > 0). Earthquakes nucleate spontaneously, with inertia approximated with radiation damping.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Cumulative slip profiles plotted over a 1,200 year period in blue every 5 years during interseismic loading and in red every second during quasi-dynamic rupture. Results were obtained using the BICyclE code for (a) BP1 with a cell size of 50 m, (b) BP2 with a cell size of 25 m, (c) BP2 with a cell size of 100 m and (d) BP2 with a cell size of 200 m. Number of events also listed, where we define a seismic event to be one with a local slip rate > 0.01 m/s separated by aseismic periods of at least 15 s.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Online platform for the SEAS working group. (Left) Home page for our website. (Top right) Currently available benchmarks. (Bottom right) Examples of BP1 model submissions.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Long-term behavior of BP1 models. (a) Shear stress and (b) slip rates at the depth of 7.5 km in models with different outer boundary conditions (BC) and computational domain sizes. (c) Shear stress and (d) slip rates at depth of 7.5 km in models with sufficiently large computational domain sizes. Legend labels indicate model names followed by information on BC and domain size, namely, ($L_x/L_z$/BC) for FDM/FEM, and ($L_z$/BC) or (HS, half-space) for BEM. BC1 and BC2 refer to the far-field free surface or displacement BC and BC3 refers to the peridic BC.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

Figure 6: Coseismic behavior of BP1 models. Coseismic phase during the 8th event in Figure 5 is shown. Models with smaller computational domain sizes show discrepancies in (a) shear stresses at 12.5 km depth and (b) slip rates at 7.5 km depth. Models with sufficiently large computational domain sizes are compared for (c) shear stresses at 12.5 km depth and (d) slip rates at 7.5 km depth. Time series are aligned relative to the rupture initiation time at the depth of 12.5 km in each model. Note that the half-space solution **luo** is the same in (b) and (d) and serves as a reference. The surface reflection phase is marked by a black arrow.

<figure><img src="image_8.png" /></figure>

Figure 7: Comparison of best-resolved BP2 models (cell size of ~25 m). (a) Long-term evolution of slip rates at depth of 9.6 km; (b) coseismic evolution of slip rates at the depth of 9.6 km for the 10th large events in the sequence (marked in gray in (a)). Time series are aligned relative to the rupture initiation time at the depth of 12 km in each model.
<figure><img src="image_9.png" /><figcaption>Figure 8: Increasing discrepancy in BP2 models due to an increased cell size of (a) 25 m, (b) 50 m, (c) 100 m, and (d) 200 m. Time evolution of shear stress at the depth of 9.6 km during the first 600 years is shown for models from three groups (abrahams, barbot, and liu).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Effect of model resolution on earthquake patterns. Distribution of (a, top row) earthquake sizes and (b, middle row) of total seismic moment release per unit length, M (in unit of N) and (c, bottom row) frequency-size relation. Models from two groups (jiang and cattania) are compared. The corresponding cell size ($\Delta z$) and total seismic event numbers ($N_s$) are marked in the titles. Seismic moment $M$ refers to the seismic moment of each earthquake; total seismic moment $M_t$ refers to the sum of moment release for all earthquakes within each magnitude bin. $N_s$ in (c) refers to the number of seismic events with moment above the corresponding $M$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Effect of model resolution on seismic-aseismic slip partitioning over depth. Depth distribution of the ratio of total seismic moment release to total moment release, $R_s$, is shown by solid lines. The ratio between seismic moment due to surface-breaching earthquakes (with surface slip greater than 0.1 m) to total moment release is indicated by dashed lines. Simulations with different resolutions are shown, with the same color for each modeling group. Note that not all groups have simulation results for all resolutions.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Effect of model resolution on recurrence intervals of large surface-breaching events. The vertical lines indicate the range of recurrence interval values, with the median value marked as dots.</figcaption></figure>