Numerical Modeling of Earthquake Cycles Based On Navier-Stokes Equations With Viscoelastic-Plasticity Rheology

[YangHaibin1102]

-> submitter ORCID (or name)

0000-0002-8628-3704

-> slug

yang-2023-eqcycles

-> license

CC-BY-4.0

-> alternative license URL

No response

-> model category

model published in study, community benchmark, forward model

-> model status

completed

-> associated publication DOI

https://doi. org/10.1029/2023GC010872

-> model creators

Haibin Yang, 0000-0002-8628-3704
Louis Moresi, 0000-0003-3685-174X
Huihui Weng, 0000-0002-2936-2342
Julian Giordani, 0000-0003-4515-9296

-> title

Numerical Modeling of Earthquake Cycles Based On Navier-Stokes Equations With Viscoelastic-Plasticity Rheology

-> description

The numerical modeling method for long-term tectonic deformations
averages out the co-seismic fault displacement into thousands to tens of thousands of years, and neglects
near-fault damages of earthquakes; therefore, it may not be able to decipher fault activities in detail. Software
simulating earthquake rupture dynamics may not have a good estimation of background stress due to longterm
tectonic deformations. In this study, we develop a numerical framework that embeds earthquake rupture
dynamics into a long-term tectonic deformation model by adding inertial terms and using highly adaptive
time-stepping that can capture deformation at plate-motion rates as well as individual earthquakes. The inertia
term, which is neglected in long-term large-scale modeling methods, is considered to simulate the dynamic
rupture processes. The rate-and-state frictional relationship for co-seismic fault slip is implemented in
viscoelastic-plastic earth. Benchmarks of viscous flow, viscoelastic wave propagation and earthquake cycle
simulations are tested. Based on these benchmarks, we undertake a generic study of a thrust fault in crust.
We find that lower crustal rheology affects the periodic time of characteristic large earthquake cycles and the
inter-seismic free surface movement. Cratons with a relatively strong lower crust due to lower temperature
remain two peaks in surface uplift profiles around the fault zone for thousands of years after one characteristic
earthquake, which help identify active faults in cratons.

-> abstract

Visco-elastic-plastic modeling approaches for long-term tectonic deformation assume that
co-seismic fault displacement can be integrated over 1000s–10,000s years (tens of earthquake cycles) with the
appropriate failure law, and that short-timescale fluctuations in the stress field due to individual earthquakes
have no effect on long-term behavior. Models of the earthquake rupture process generally assume that the
tectonic (long-range) stress field or kinematic boundary conditions are steady over the course of multiple
earthquake cycles. This study is aimed to fill the gap between long-term and short-term deformations by
modeling earthquake cycles with the rate-and-state frictional (RSF) relationship in Navier-Stokes equations.
We reproduce benchmarks at the earthquake timescale to demonstrate the effectiveness of our approach. We
then discuss how these high-resolution models degrade if the time-step cannot capture the rupture process
accurately and, from this, infer when it is important to consider coupling of the two timescales and the level of
accuracy required. To build upon these benchmarks, we undertake a generic study of a thrust fault in the crust
with a prescribed geometry. It is found that lower crustal rheology affects the periodic time of characteristic
earthquake cycles and the inter-seismic, free-surface deformation rate. In particular, the relaxation of the
surface of a cratonic region (with a relatively strong lower crust) has a characteristic double-peaked uplift
profile that persists for thousands of years after a major slip event. This pattern might be diagnostic of active
faults in cratonic regions.

-> scientific keywords

Earthquake Cycles, Navier-Stokes, Viscoelastic-Plasticity, Cratonic Earthquakes

-> funder

https://ror.org/01h0zpd94, #42030306
https://ror.org/04yx6dh41

-> model embargo?

No response

-> include model code ?

• yes

-> model code/inputs DOI

https://doi.org/10.5281/ zenodo.8251145

-> model code/inputs notes

No response

-> include model output data?

• yes

-> data creators

Haibin Yang, 0000-0002-8628-3704
Louis Moresi, 0000-0003-3685-174X
Huihui Weng, 0000-0002-2936-2342
Julian Giordani, 0000-0003-4515-9296

-> model output data DOI

No response

-> model output data notes

No response

-> model output data size

2.8Gb

-> software framework DOI/URI

https://doi.org/10.5281/zenodo.3975252

-> software framework source repository

No response

-> name of primary software framework (e.g. Underworld, ASPECT, Badlands, OpenFOAM)

Underworld

-> software framework authors

No response

-> software & algorithm keywords

Python, Finite element, Particle-in-cell

-> computer URI/DOI

No response

-> add landing page image and caption

Figure 5. Evolution of the stress at the reference point (0, 0, −10 km) (a), maximum slip rate along the entire fault zone (b)
and the adaptive time step used in simulation (c) for the reference model.

-> add an animation (if relevant)

No response

-> add a graphic abstract figure (if relevant)

No response

-> add a model setup figure (if relevant)

Figure 4. The benchmark model BP5 for 3D sequence of earthquakes and
aseismic slip modeling. (a) A vertical planar fault is embedded in the middle
of a homogenous, isotropic half-space with a free surface at z = 0. Fault
behavior is controlled by the rate-and-state friction law. A periodic boundary
condition is applied in y direction. (b) The velocity-weakening (VW) region
(dark and light blue) is located within a transition zone (white), outside of
which is the velocity-strengthening (VS) region (gray). In y and z directions,
the frictional domain and VW region are (L2, L3) and (l, w), respectively. An
initial nucleation zone (dark blue square with a width of w) is designed at the
left end of the VW region.

-> add a description of your model setup

The BP5 benchmark is first simulated with a reference model size of 96 km (L1) × 100 km (L2) × 30 km
(L3) by 128 × 64 × 64 quadrilateral bilinear elements. Results from the SEAS code comparison platform (https://
strike.scec.org/cvws/seas/) with the mesh resolution of 1000 m are selected for comparison in this study (Table 2).
The stress at a depth of 10 km in the middle point along the fault strike (y = 0) and the maximum slip rate along
the entire fault are tracked for comparison (Figure 5). 0.1 m/s is taken as a threshold of fault slip rate to mark
the earthquake initiation.

Please provide any feedback on the model submission process?

No response

[m-te-bot]

Model Report

Thank you for submitting.

Using Github actions, we have regenerated a report summarising information about your model

• Please check the report below, including the Errors and Warnings section

• You can update any information, by editing the markdown file at the top of the issue

• these edits will trigger the report will be regenerated

• once you are satisfied with the results, please add a review requested Review requested label

Parsed data

Section 1: Summary of your model

Model Submitter:

Haibin Yang (0000-0002-8628-3704)

Model Creator(s):

• Haibin Yang (0000-0002-8628-3704)
• Louis-Noel Moresi (0000-0003-3685-174X)
• Huihui Weng (0000-0002-2936-2342)
• Julian Giordani (0000-0003-4515-9296)

Model slug:

yang-2023-eqcycles

(this will be the name of the model repository when created)

Model name:

Numerical Modeling of Earthquake Cycles Based On Navier-Stokes Equations With Viscoelastic-Plasticity Rheology

License:

Creative Commons Attribution 4.0 International

Model Category:

• model published in study
• community benchmark
• forward model

Model Status:

• completed

Associated Publication title:

Numerical Modeling of Earthquake Cycles Based On Navier‐Stokes Equations With Viscoelastic‐Plasticity Rheology

Short description:

The numerical modeling method for long-term tectonic deformations
averages out the co-seismic fault displacement into thousands to tens of thousands of years, and neglects
near-fault damages of earthquakes; therefore, it may not be able to decipher fault activities in detail. Software
simulating earthquake rupture dynamics may not have a good estimation of background stress due to longterm
tectonic deformations. In this study, we develop a numerical framework that embeds earthquake rupture
dynamics into a long-term tectonic deformation model by adding inertial terms and using highly adaptive
time-stepping that can capture deformation at plate-motion rates as well as individual earthquakes. The inertia
term, which is neglected in long-term large-scale modeling methods, is considered to simulate the dynamic
rupture processes. The rate-and-state frictional relationship for co-seismic fault slip is implemented in
viscoelastic-plastic earth. Benchmarks of viscous flow, viscoelastic wave propagation and earthquake cycle
simulations are tested. Based on these benchmarks, we undertake a generic study of a thrust fault in crust.
We find that lower crustal rheology affects the periodic time of characteristic large earthquake cycles and the
inter-seismic free surface movement. Cratons with a relatively strong lower crust due to lower temperature
remain two peaks in surface uplift profiles around the fault zone for thousands of years after one characteristic
earthquake, which help identify active faults in cratons.

Abstract:

Visco-elastic-plastic modeling approaches for long-term tectonic deformation assume that
co-seismic fault displacement can be integrated over 1000s–10,000s years (tens of earthquake cycles) with the
appropriate failure law, and that short-timescale fluctuations in the stress field due to individual earthquakes
have no effect on long-term behavior. Models of the earthquake rupture process generally assume that the
tectonic (long-range) stress field or kinematic boundary conditions are steady over the course of multiple
earthquake cycles. This study is aimed to fill the gap between long-term and short-term deformations by
modeling earthquake cycles with the rate-and-state frictional (RSF) relationship in Navier-Stokes equations.
We reproduce benchmarks at the earthquake timescale to demonstrate the effectiveness of our approach. We
then discuss how these high-resolution models degrade if the time-step cannot capture the rupture process
accurately and, from this, infer when it is important to consider coupling of the two timescales and the level of
accuracy required. To build upon these benchmarks, we undertake a generic study of a thrust fault in the crust
with a prescribed geometry. It is found that lower crustal rheology affects the periodic time of characteristic
earthquake cycles and the inter-seismic, free-surface deformation rate. In particular, the relaxation of the
surface of a cratonic region (with a relatively strong lower crust) has a characteristic double-peaked uplift
profile that persists for thousands of years after a major slip event. This pattern might be diagnostic of active
faults in cratonic regions.

Scientific Keywords:

• Earthquake Cycles
• Navier-Stokes
• Viscoelastic-Plasticity
• Cratonic Earthquakes

Funder(s):

• National Natural Science Foundation of China (https://ror.org/01h0zpd94)
• National Computational Infrastructure (https://ror.org/04yx6dh41)

Section 2: your model code, output data

No embargo on model contents requested

Include model code:

True

Include model output data:

True

Section 3: software framework and compute details

Software Framework DOI/URL:

Found software: Underworld

Name of primary software framework:

Underworld

Software & algorithm keywords:

• Python
• Finite element
• Particle-in-cell

Section 4: web material (for mate.science)

Landing page image:

Caption:
Figure 5. Evolution of the stress at the reference point (0, 0, −10 km) (a), maximum slip rate along the entire fault zone (b)
and the adaptive time step used in simulation (c) for the reference model.

Animation:

Filename:

Graphic abstract:

Filename:

Model setup figure:

Caption:
Figure 4. The benchmark model BP5 for 3D sequence of earthquakes and
aseismic slip modeling. (a) A vertical planar fault is embedded in the middle
of a homogenous, isotropic half-space with a free surface at z = 0. Fault
behavior is controlled by the rate-and-state friction law. A periodic boundary
condition is applied in y direction. (b) The velocity-weakening (VW) region
(dark and light blue) is located within a transition zone (white), outside of
which is the velocity-strengthening (VS) region (gray). In y and z directions,
the frictional domain and VW region are (L2, L3) and (l, w), respectively. An
initial nucleation zone (dark blue square with a width of w) is designed at the
left end of the VW region.
Description: The BP5 benchmark is first simulated with a reference model size of 96 km (L1) × 100 km (L2) × 30 km
(L3) by 128 × 64 × 64 quadrilateral bilinear elements. Results from the SEAS code comparison platform (https://
strike.scec.org/cvws/seas/) with the mesh resolution of 1000 m are selected for comparison in this study (Table 2).
The stress at a depth of 10 km in the middle point along the fault strike (y = 0) and the maximum slip rate along
the entire fault are tracked for comparison (Figure 5). 0.1 m/s is taken as a threshold of fault slip rate to mark
the earthquake initiation.

Errors and Warnings

Associated Publication
Error fetching metadata with application/ld+json from https://api.crossref.org/works/https://doi. org/10.1029/2023GC010872: 406 Client Error: Not Acceptable for url: https://api.crossref.org/works/https://doi.%20org/10.1029/2023GC010872
Software Framework DOI/URI
doi.org metadata record succesfully extracted in json-ld format
Software Repository
Warning: no repository URL provided.
Submitter
ORCID metadata record succesfully extracted in json-ld format

Model creators
ORCID metadata record succesfully extracted in json-ld format
ORCID metadata record succesfully extracted in json-ld format
ORCID metadata record succesfully extracted in json-ld format
ORCID metadata record succesfully extracted in json-ld format
Could not parse Embargo date. Check format is
Model code/inputs DOI
404 Client Error: Not Found for url: https://doi.org/10.5281/%20zenodo.8251145
Model code/inputs notes
Warning: No notes provided.
Data creators
ORCID metadata record succesfully extracted in json-ld format
ORCID metadata record succesfully extracted in json-ld format
ORCID metadata record succesfully extracted in json-ld format
ORCID metadata record succesfully extracted in json-ld format
Model output DOI
Warning: No DOI/URI provided.
Model data notes
Warning: No notes provided.
Computer URI/DOI
Warning: No URI/DOI provided.
Animation
Warning: No animation uploaded.

Graphic abstract
Warning: No image uploaded.

Next steps

• once the model_reviewers team has approved the model, we will create a repository for your model

[m-te-bot]

Review Requested

A review of this submission has been requested from @ModelAtlasofTheEarth/model_reviewers

[m-te-bot]

Model repository created at https://github.com/ModelAtlasofTheEarth/yang-2023-eqcycles

[m-te-bot]

Model repository created at https://github.com/ModelAtlasofTheEarth/yang-2023-eqcycles