Enter Zoom Meeting


Physics-based earthquake modeling and engineering

Numerical modeling of earthquakes provides new approaches to apprehend the physics of earthquake rupture and the seismic cycle, seismic wave propagation, fault zone evolution and seismic hazard assessment.
Recent advances in numerical algorithms and increasing computational power enable unforeseen precision and multi-physics components in physics-based earthquake simulation but also pose challenges in terms of fully exploiting modern supercomputing infrastructure, realistic parameterization of simulation ingredients and the analysis of large synthetic datasets while advances in laboratory experiments link earthquake source processes to rock mechanics.
This session aims to bring together modelers and data analysts interested in the physics and computational aspects of earthquake phenomena and earthquake engineering. We welcome studies focusing on all aspects of seismic hazard assessment and the physics of earthquakes - from slow slip events, fault mechanics and rupture dynamics, to wave propagation and ground motion analysis, to the seismic cycle and inter seismic deformation - and studies which further the state-of-the art in the related computational and numerical aspects.

Co-organized by NH4/TS3
Convener: Alice-Agnes GabrielECSECS | Co-conveners: Jean-Paul Ampuero, Hideo Aochi
| Mon, 23 May, 14:05–14:50 (CEST), 15:10–18:30 (CEST)
Room D3

Mon, 23 May, 13:20–14:50

Chairpersons: Alice-Agnes Gabriel, Hideo Aochi, James Biemiller

Mikihiro Kawabata et al.

Dehydration embrittlement was proposed to account for intermediate or deep earthquakes (e.g., Raleigh and Paterson, 1965). Many researchers have investigated the frictional instability induced by dehydration of hydrous minerals, such as gypsum (e.g., Milsch and Scholz, 2005; Brantut et al., 2011; Leclère et al., 2016). In addition, time dependence of dehydration of hydrous minerals has been studied based on reaction kinetics (e.g., Sawai et al., 2013). Since kinetics controls the dehydration rate, the effect of dehydration-derived pore fluid pressure on the mechanical strength of rocks can also be represented by kinetics. However, there is no experimental study to quantitatively investigate how pore fluid pressure builds up and controls the mechanical strength of fault gouges in terms of kinetics. Here, we derived time function of pore fluid pressure based on dehydration kinetics of simulated gypsum (bassanite) gouges. First, we conducted friction experiments of simulated gypsum gouges using gas apparatus under eight different conditions of pressures from 10 MPa to 200 MPa and temperatures from room temperature to 180 °C, spanning dehydration condition of gypsum. Each stress-strain curve showed stick-slip behaviors with almost constant stress drops and recurrence intervals depending on the effective pressures under the conditions of room temperature (RT): larger stress drops and longer intervals for higher effective stresses. On the other hand, stress drops and recurrence intervals gradually decrease with time under 200 MPa and 110 °C, close to the dehydration boundary. These results suggested that the elevated pore fluid pressure by dehydration decreases effective pressure and reduces the stress drops and the intervals. We tested this hypothesis as follows. Microstructural observations illuminated marked development of Riedel shears (R1 shear) in samples deformed under the stability field of gypsums (RT and 70 °C), while scarce development of Riedel shears in the sample deformed under 110 °C, being consistent with Leclère et al. (2016)’s observations on that the elevated pore pressure suppress the development of Riedel shears. Based on the equation of state for water (He and Zoller, 1991), we calculated the porosity of the sample deformed under 110 °C. Although the estimated value was smaller than that obtained from dehydration under hydrostatic conditions (Bedford et al., 2017), this result indicates that shear compaction may have occurred due to deformation caused by higher differential stress. Considering that the decrease in effective pressure modulates the amount of stress drops and recurrence intervals, we analyzed frictional coefficients with Mohr’s circle assuming pore fluid pressure. The estimated value of about 0.6 is consistent with Byerlee (1978)’s law. Based on the results, we created a time function for evolution of pore fluid pressure controlled by Avrami-type dehydration kinetics (Avrami, 1940). The estimated Avrami exponent, the important parameter for crystallization, of 3.121 indicated that the dehydration proceeded with nucleation and three-dimensional growth. This function enables more accurate prediction of pore fluid pressure evolution controlled by dehydration kinetics and may contribute to better understanding the effect of hydrous minerals on frequency of intermediate and deep earthquakes.

How to cite: Kawabata, M., Sasaki, Y., Iwasaki, M., Shiraishi, R., Muto, J., and Nagahama, H.: Time-varying stick-slip behaviors described by dehydration kinetics of gypsum, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-11092, https://doi.org/10.5194/egusphere-egu22-11092, 2022.


James Moore et al.

Advances in modelling and access to InSAR and GNSS observations have highlighted the role that rheological heterogeneities play in postseismic deformation. Here we discuss three recent studies (Muto et al. 2019, Sambuddha et al. 2022, and Takada et al. in prep) following the 2011 Tohoku-Oki and 2008 Iwate-Miyagi earthquakes, which reveal both localised and along-strike rheological heterogeneities. We construct a self-consistent physical model of the postseismic deformation for these two events using the Unicycle code (Moore et al. 2019, Barbot, Moore, and Lambert 2017), with which we consider coupled fault slip and viscoelastic flow utilising laboratory-derived constitutive laws to simulate the time series of geodetic observations. All three studies illuminate a crustal low viscosity rheological heterogeneity in the vicinity of Mt Kurikoma / Mt Naruko. This is perhaps to be expected, given the proximity to known active volcanic centres, and is commensurate with observations following the 2016 Kumamoto earthquake (Moore et al. 2017) where we found low-viscosity anomalies beneath Mt Aso and Mt Kuju. However, the heterogeneities the data reveal are not restricted to known volcanic regions, because our results also suggest along-arc heterogeneity in the forearc mantle rheology of north-eastern Japan; specifically we find a narrower cold nose in the Miyagi region and wider for the Fukushima forearc. We also find evidence of interaction between the localized crustal heterogeneity and afterslip in both events, highlighting the importance of addressing mechanical coupling for long-term studies of postseismic relaxation. Variations in rheological properties in the lithosphere are not restricted to viscous and thermal effects, and observations of the Iwate-Miyagi earthquake suggest elastic heterogeneities may also play a role. We therefore conclude by presenting expressions for computing displacements and stress due to localised (faulting) and distributed inelastic deformation in heterogeneous elastic spaces with piece-wise constant homogeneous elastic subregions (Sato & Moore 2022), and their application in the context of the seismic cycle.


Muto J, Moore J D P, Barbot S, Iinuma T, Ohta Y, Horiuchi S, Hikaru I, 2019. Coupled afterslip and transient mantle flow after the 2011 Tohoku earthquake. Science Advances

Dhar S, Muto J, Ito Y, Muira S, Moore J D P, Ohta Y, Iinuma T, 2022. Along-Arc Heterogeneous Rheology Inferred from Postseismic deformation of the 2011 Tohoku-oki Earthquake.

Moore J D P, Barbot S, Feng L, Hang Y, Lambert V, Lindsey E, Masuti S, Matsuzawa T, Muto J, Nanjundiah P, Salman R, Sathiakumar S, & Sethi H, 2019. jdpmoore/unicycle: Unicycle. In Coupled afterslip and transient mantle flow after the 2011 Tohoku earthquake, Science Advances 2019. Zenodo. https://doi.org/10.5281/zenodo.5688288

Barbot S, Moore J D P, Lambert V, 2017. Displacements and stress associated with distributed anelastic deformation in a half-space. BSSA

Moore J D P, Yu H, Tang C, Wang T, Barbot S, Peng D, Masuti S, Dauwels J, Hsu Y, Lambert V, Nanjundiah P, Wei S, Lindsey E, Feng L, Shibazaki B, 2017. Imaging the distribution of transient viscosity after the 2016 Mw7.1 Kumamoto earthquake. Science

Sato D, Moore J D P, 2022. Displacements and stress associated with localised and distributed inelastic deformation with piecewise-constant elastic variations.

How to cite: Moore, J., Dhar, S., Muto, J., Sato, D., and Takada, Y.: The role of rheological heterogeneities in postseismic deformation, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-7627, https://doi.org/10.5194/egusphere-egu22-7627, 2022.


Eyup Sopaci and Atilla Arda Özacar

The triggering mechanism of earthquakes and their synchronization in time and space can be considered the two sides of the same coin. Our previous studies on earthquake triggering reveal sensitive parameters affecting the triggering mechanism using simple spring slider systems. We pursue our previous analyses by considering a simulation set-up for synchronizing three strong asperity patches on a vertically oriented strike-slip fault with initial slip heterogeneity separated by barriers and strong creeping regions at the edges. This analogy intends to explore earthquake synchronization in time and mimic observed sequences of large earthquakes that ruptured most of the North Anatolian Fault within short time intervals. Using the quasi-dynamic and full-dynamic pseudo-spectral Fast Fourier Transform (FFT) method, we apply a periodic fault model governed with Rate-and-State Friction (RSF) law embedded in a 2.5D continuum. Simulation results so far using the quasi-dynamic approach revealed that the earthquake synchronization is mainly affected by direct velocity effect parameters, barrier dimension/properties, and RSF law (aging and slip law), particularly the weakening terms. Lower direct velocity effect parameters, state evolutions with a stronger weakening term such as slip law, and shorter barrier lengths promote better synchronization. In this respect, we observed fast, slow, or no synchronization depending on the parameter sets. It is also worth noting that slip localizes in the continuum at small critical slip distances, which cannot be inferred from simple 1D models, suggesting the size dependence. In order to minimize inherent non-uniqueness and uncertainties, the same set-up will also be simulated with the full-dynamic approach in which wave-mediated stress transfer is taken into account, and the long-term earthquake histories will be correlated with case-specific simulations.

How to cite: Sopaci, E. and Özacar, A. A.: Do Large Earthquakes along Major Faults Synchronize in Time?, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-508, https://doi.org/10.5194/egusphere-egu22-508, 2022.


Amir Sagy et al.

Geological and geophysical observations indicate that fault geometry is nonplanar, includes irregularities in all directions at many scales. The geometrical heterogeneity of faults is particularly critical during the interseismic stage of the earthquake cycle because it perturbs the stress field and thus affects the rupture nucleation along the fault zone and around it. We present a new analytical solution for the static stress field around a rough interlocked interface obtained under compressional stresses, and discuss its applications to faulting and seismic hazards. The model outputs are the local stress field and the Failure-Ratio, defined here as the susceptibility to failure of the bulk material around the interface. The calculations are then obtained by the following steps: First, the interface geometry is represented by a Fourier series. Then, the stress components around the irregular interface are calculated analytically using perturbation theory for any two dimensional far-field stress tensor. Finally, the Failure Ratio at any location near the interface is estimated by adopting a Coulomb failure criterion for the bulk material.

The model results can be applied to faulting mechanics because they demonstrate how the elastic stress field around rough fault is controlled by the geometry and by the tectonic stresses. We find that under a given tectonic stress state, stress heterogeneity increases with roughness. Therefore, some zones near rough faults are expected to yield at lower tectonic shear stress comparing to zones nearby smooth ones. However, the magnitudes of these events are expected to be relatively small, as they nucleate under relatively low tectonic stresses and fail as they propagating immediately to a stress shadow. This stress distribution promotes small seismic events near rough fault and therefore we suggest that increasing heterogeneity of the surface, contributes to increasing of the b-value in Gutenberg-Richter earthquakes distribution.

We compare the model predictions with results of experiments performed on rough rock surfaces and find good agreement between the locations of off-fault deformation zones and the calculated high Failure-Ratio values. We further test the model implications for stresses and failure around a natural fault system – the San Andreas Fault and find a first-order agreement between Failure-Ratio values and earthquake distribution around this fault system. We conclude that the proposed analytical approach is a useful and practical tool for evaluating the contribution of fault geometry to the seismic hazard potential around it.


How to cite: Sagy, A., Morad, D., H. Hatzor, Y., and Lyakhovsky, V.: The onset of faulting around geometrically irregular faults, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-5309, https://doi.org/10.5194/egusphere-egu22-5309, 2022.


Yuval Tal

Earthquakes occur by sudden slippage along pre-existing faults via a frictional instability. Laboratory-derived rate and state friction laws have emerged as powerful tools for investigating the mechanics of earthquakes. Two types of state‐variable evolution laws are commonly used to fit the experimental data, the aging and slip laws. The aging evolution law has been used extensively to model the earthquake cycle, including the nucleation, dynamic rupture propagation and arrest, and interseismic period. The slip law, which generally provides a better fit to rock friction experiments, has rarely been used in simulations of the whole seismic cycle. In addition, faults are zones with complex internal structure and non-planar geometry, which also affect the rupture process during the seismic cycle.

In this study, I examine the effects of fault geometry, state evolution law, and friction parameters on the earthquake source process with fully dynamic 2-D simulations of earthquake sequences on planar and non-planar faults. The numerical approach accounts for all stages in the seismic cycle and enables modeling slip that is comparable to the minimum wavelength of roughness. I test the statistics of the events in terms of static source parameters and analyze in detail the rupture process during the nucleation and dynamic propagation stages. For the same friction parameters and fault geometry, the slip law results in a more rapid weakening of the friction coefficient than the aging law. That leads to ruptures with smaller nucleation sizes, larger slip rates, and larger rupture speeds for the slip law, including transition to supershear. With the aging law, a small level of fault roughness is enough to introduce considerable complexity into the rupture process, with larger amount of aseismic slip and larger variability in earthquake sizes.  For the same level of roughness, those effects are significantly smaller in the case of the slip law.

How to cite: Tal, Y.: The Seismic Cycle on Rate and State Faults with Different Evolution Laws and Fault Geometries , EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-2368, https://doi.org/10.5194/egusphere-egu22-2368, 2022.


Ludovico Manna et al.

We present a study on the dependence of the frictional properties of a fault rock on its degree of damage. The purpose is therefore to gain insight into frictional sliding, the governing force that controls earthquake nucleation, propagation and arrest. The focus on this topic is to try to find a reason for the experimental evidence that the friction coefficient seems to be almost independent on lithology. A possible explanation to investigate through the numerical modelling could be that the frictional properties of a realistic fault rock depend mostly on the concentration of micro- to macroscopic cracks and/or of lamellar phyllosilicates in the host rock, rather than on the composition of its bulk materials. The formalism of the Linear Elastic Fracture Mechanics (LEFM) can quantitatively reproduce the stresses and the strains on the interface propagating frictional rupture. The purpose is to use a Finite Element Method (FEM) numerical code in order to simulate the plane strain elastic deformation of a two-dimensional medium crossed by elliptical fractures and weak anisotropic inclusions. The analysis of the distribution and orientation of the stresses resulting from the interaction of a system of randomly oriented elliptical fractures under different loading conditions could provide information on the onset and propagation of frictional ruptures, such as real contact area reduction, slip velocity, number and length of global sliding precursors. The magnitude and orientation of the principal stresses around the tips of elliptical voids are crucial for the understanding of fracture coalescence and frictional reactivation of shear cracks in an elastic rock, which in turn is one of the main factors that govern the seismic cycle of natural faults. 

How to cite: Manna, L., Dabrowski, M., Maino, M., Casini, L., Reali, A., and Toscani, G.: On the relation between the coefficient of friction of a fault and the variation of damage degree on the host rock: a numerical approach, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-4709, https://doi.org/10.5194/egusphere-egu22-4709, 2022.

Mon, 23 May, 15:10–16:40

Chairpersons: Hideo Aochi, Casper Pranger, Jorge N. Hayek

Ranjith Kunnath

The spectral boundary integral equation (SBIE) method is widely used for numerical modeling of earthquake ruptures at a planar interface between two elastic half-spaces. It was originally proposed by Geubelle and Rice (1995) based on the boundary integral formulation of Budiansky and Rice (1979). The distinguishing feature of the formulation is that it involves performing elastodynamic space-time convolution of the displacement discontinuities at the interface between the two solids. The method was extended to bi-material interfaces by Geubelle and Breitenfeld (1997) and Breitenfeld and Geubelle (1998). An alternative boundary integral formulation to that of Budiansky and Rice (1979) is that of Kostrov (1966), where the elastodynamic space-time convolution is done of the tractions at the interface between the two solids. A SBIE method based on the latter formulation was proposed by Ranjith (2015) for plane strain. In the present work, the SBIE method for antiplane strain based on the formulation of Kostrov (1966) is proposed and compared with other approaches. Illustrations of the use of the method for simulating dynamic antiplane ruptures at bi-material interfaces are given.

How to cite: Kunnath, R.: A new spectral boundary integral equation method for antiplane problems, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-9, https://doi.org/10.5194/egusphere-egu22-9, 2022.


Duo Li et al.

Earthquake fault zones are more complex, both geometrically and rheologically, than an idealized infinitely thin plane embedded in linear elastic material. Field and laboratory measurements have revealed intense fault weakening induced by flash heating and melting on natural fault (Di Toro et al., 2006; Goldsby & Tullis, 2011) and complex fault zone structure involving both tensile and shear fractures spanning a wide spectrum of length scales (e.g., Mitchell & Faulkner, 2009). Previous 2D numerical models explicitly accounting for off-fault fractures have demonstrated important feedback with rupture dynamics and ground motions (e.g., Thomas & Bhat 2018, Okubo et al., 2019). However, numerical studies of thermal-related weakening mechanisms usually avoid frictional melting due to the lack of the solid-fluid phase transition. 

In the work of Gabriel et al. (2021), we have presented our first-order hyperbolic and thermodynamically compatible mathematical model, namely the GPR model (Godunov & Romenski, 1972; Romenski, 1988), combined with a diffuse crack representation to incorporate finite strain nonlinear material behavior, natural complexities and multi-physics coupling within and outside of fault zones into dynamic earthquake rupture modeling. We compare our novel diffuse interface fault models of kinematic cracks, spontaneous dynamic rupture, and dynamically generated off-fault shear cracks to sharp interface reference models. Pre-damaged faults, as well as dynamically induced secondary cracks are therein described via a scalar function indicating the local level of material damage (Tavelli et al., 2020); arbitrarily complex geometries are represented via a diffuse interface approach based on a solid volume fraction function (Tavelli et al., 2019). 

Here we further extend the diffuse crack representation to more complicated scenarios including severe dynamic fault zone weakening as activated by flash heating, the effect of locally melting rocks, and off-fault cracks with complex topology in 3D materials, by taking advantage of adaptive Cartesian meshes (AMR) embedded in the extreme-scale hyperbolic PDE solver ExaHyPE (Reinarz et al., 2019). We intend to compare our thermally-weakened rupture in diffused fault zone with the semi-analytical thermal pressurization weakening implemented in the linear elastodynamic rupture on an infinitely-thin fault surface, using SeisSol (https://github.com/SeisSol). We will further qualitatively verify our model using the up-to-date observations in the 2020 M8.2 Chignik, Alaska, to illustrate the importance of thermal weakening on relatively deeper faults.

Our approach is part of the TEAR ERC project (www.tear-erc.eu) and will potentially allow to fully model volumetric fault zone shearing during earthquake rupture, which includes spontaneous partition of fault slip into intensely localized shear deformation within weaker (possibly cohesionless/ultracataclastic) fault-core gouge and more distributed damage within fault rocks and foliated gouges.

How to cite: Li, D., Gabriel, A.-A., Chiocchetti, S., Tavelli, M., Peshkov, I., Romenski, E., and Dumbser, M.: A unified model for thermally-activated fault weakening during nonlinear dynamic earthquake rupture and off-fault fracturing in 3D diffuse fault zones, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-13551, https://doi.org/10.5194/egusphere-egu22-13551, 2022.


Kyungjae Im and Jean-Philippe Avouac

Earthquakes often come in clusters formed of foreshock-mainshock-aftershock sequences. This clustering is generally thought to result from a cascading process which is commonly modeled using either the phenomenological ETAS model or a stress-based model assuming an earthquake nucleation process governed by Coulomb stress changes and Rate-and-State friction (CRS). In this work, we numerically investigated the foreshock and aftershock sequence in a discrete fault network with rate and state friction law and compared the result with the ETAS model. We set a fault zone consisting of dense fault segments and an off-fault area consisting of sparsely distributed smaller faults in the simulation domain. The CRS simulations are conducted 100 times with 1000 discrete faults with randomly generated fault location, initial velocity, and fault length within the weighted distribution, yielding a Gutenberg-Richter law. The simulations produce realistic foreshocks and aftershocks sequences. Aftershocks occur in the area of increased Coulomb stress and decay following Omori law as observed in nature. Individual foreshock sequences do not show a clear trend, but once stacked, they show an apparent inverse-Omori law acceleration. The prediction from our CRS model can be fitted with the ETAS model. This is not surprising since ETAS incorporates the Omori and Gutenberg-Richter laws. However, our CRS model predicts significantly more foreshocks than would be expected from the ETAS model. This results from the fact that the triggering productivity is lower in the aftershock sequence than in the foreshocks due to the depletion of critically stressed faults in our simulations. In other words, the ETAS is not compatible with the CRS model because Coulomb stress changes result in a time advance (if positive) or delay (if negative). This clustering process is fundamentally different from the additive process assumed in ETAS. As a result, the claim made that foreshocks more frequent than expected based on ETAS imply pre-seismic slip might be incorrect. It could alternatively be a manifestation of the nucleation process.

How to cite: Im, K. and Avouac, J.-P.: Numerical Modeling of Cascading Foreshocks and Aftershocks in Discrete Fault Network, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-3110, https://doi.org/10.5194/egusphere-egu22-3110, 2022.


Dave May et al.

Simulations of sequences of earthquakes and aseismic slip (SEAS) including more than one fault, complex geometries and elastic heterogeneities are challenging. We present a symmetric interior penalty discontinuous Galerkin (SIPG) method accounting for the complex geometries and heterogeneity of the subsurface. The method accommodates two- and three-dimensional domains, is of arbitrary order, handles sub-element variations in material properties and supports isoparametric elements, i.e. high-order representations of the exterior and interior boundaries and interfaces including intersecting faults.

We provide an open-source reference implementation, Tandem, that utilises highly efficient kernels, is inherently parallel and well suited to perform high resolution simulations on large scale distributed memory architectures. Further flexibility is provided by optionally defining the displacement evaluation via a discrete Green's function, using algorithmically optimal and scalable sparse parallel solvers and preconditioners. We highlight the characteristics of the SIPG formulation via an extensive suite of verification problems (analytic, manufactured and code comparison) for elasto-static and seismic cycle problems. We demonstrate that high-order convergence of the discrete solution can be achieved in space and time for elasto-static and SEAS problems.

Lastly, we apply the method to realistic demonstration models consisting of a 2D SEAS multi-fault scenario on a shallowly-dipping normal fault with four curved splay faults, and a 3D multi-fault scenario of instantaneous displacement due to the 2019 Ridgecrest, CA, earthquake sequence. We exploit the curvilinear geometry representation in both application examples and elucidate the importance of accurate stresses (or displacement gradients) representations on-fault. Our results exploit advantages of both the boundary integral and volumetric methods and is an interesting avenue to pursue in the future for extreme scale 3D SEAS simulations.

How to cite: May, D., Uphoff, C., and Gabriel, A.-A.: A discontinuous Galerkin method for sequences of earthquakes and aseismic slip on multiple faults using unstructured curvilinear grids, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-12166, https://doi.org/10.5194/egusphere-egu22-12166, 2022.


Luca Dal Zilio et al.

There is a growing interest in understanding how geologic faults respond to transient sources of fluid. However, the spatio-temporal evolution of sequences of seismic and aseismic slip in response to pore-fluid evolution is still poorly constrained. In this study, we present H-MEC (Hydro-Mechanical Earthquake Cycles), a newly-developed two-phase flow numerical code — which couples solid rock deformation and pervasive fluid flow — to simulate how crustal stress and fluid pressure evolve during the earthquake cycle on a fluid-bearing fault structure. This unified 2D numerical framework accounts for full inertial (wave) effects and fluid flow in a finite difference method and poro-visco-elasto-plastic compressible medium with rate-dependent strength. An adaptive time stepping allows the correct resolution of both long- and short-time scales, ranging from years to milliseconds during the dynamic propagation of dynamic rupture. We present a comprehensive plane strain strike-slip setup in which we test analytical benchmarks of pore-fluid pressure diffusion from an injection point. We then investigate how pore-fluid pressure evolution and solid–fluid compressibility control sequences of seismic and aseismic slip on a finite fault width. While the onset of fluid-driven shear cracks is controlled by localized collapse of pores and dynamic self-pressurization of fluids inside the undrained fault zone, subsequent dynamic ruptures are driven by solitary pulse-like fluid pressure wave propagating at seismic speed. Furthermore, shear strength weakening associated with rapid self-pressurization of pore-fluid can account for the slip–fracture energy scaling observed in large earthquakes. This numerical framework provides a viable tool to better understand fluid-driven dynamic ruptures — either as a natural process or induced by human activities — and highlight the importance of considering the realistic hydro-mechanical structure of faults to investigate sequences of seismic and aseismic slip.

How to cite: Dal Zilio, L., Hegyi, B., Behr, W., and Gerya, T.: Fluid-driven earthquake sequences and aseismic slip in a poro-visco-elasto-plastic fluid-bearing fault structure, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-6897, https://doi.org/10.5194/egusphere-egu22-6897, 2022.


Olaf Zielke et al.

It is well established in the seismology community that geometric complexity plays an important role for a fault’s seismotectonic behavior. It affects the initiation, propagation and termination of an earthquake as well as influencing the stress-slip relationship, the size of fault segments, and the probability of multi-segment rupture. Consequently, fault geometric complexity is studied intensively and increasingly incorporated into computational earthquake rupture simulations. These efforts reveal a problem: While we may be able to constrain a natural fault’s geometry with a high level of detail at the surface (i.e., the fault trace), we cannot do the same for the buried portion of the fault -where most of the rupture takes place. How much does a fault’s seismotectonic behavior vary as a result of this epistemic uncertainty?

We address this question computationally with a physics-based multi-cycle earthquake rupture simulator (MCQsim), enabling us to investigate how (for example) earthquake recurrence, slip accumulation, magnitude-frequency distribution, and fault segmentation vary (looking at the entire fault as well as individual locations on the fault) as function of our insufficient knowledge about the fault’s geometric complexity. To simulate fault geometric complexity, we generate 2-D random fields, using the “random midpoint displacement” method (RMD), representing the fault’s non-planar, self-similar geometry. The advantage of using RMD is that it allows us to create a 2-D random field while also keeping one or more of the field’s edges at a prescribed value. Hence, this approach allows us to generate a random field to represent fault roughness while also allowing us to incorporate what is known about the fault geometry (i.e., the fault surface trace, representing one of the random field’s edges). In doing so, we can investigate how the aforementioned seismo-tectonic parameters vary as a function of fault roughness uncertainty.

For this purpose, we create 5000-year long earthquake catalogs for a 150x18km large strike slip fault that is parameterized by more than 40k fault cells (average cell size 0.07km^2), containing earthquakes with 3.5 < M < 7.8. We create these catalogs for 100 roughness realizations while keeping the simulated fault’s surface trace constant for all realizations. The results of these simulations will be presented in our presentation.

How to cite: Zielke, O., Aspiotis, T., and Mai, P. M.: Epistemic uncertainty in fault geometry effects earthquake rupture behavior, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-2414, https://doi.org/10.5194/egusphere-egu22-2414, 2022.


Jorge Nicolas Hayek Valencia et al.

Natural fault system observations feature complexity that includes damage variation from the outer damage zone to the fault core and associated rheological degradation (e.g. variation in the frictional strength and spatio-temporal slip localisation). In earthquake dynamic rupture simulations, faults are typically treated as infinitesimally thin interfaces with distinct on- versus off-fault rheologies. Commonly, such faults are explicitly represented in the discretisation of the computational domain.

Here we present a diffuse interface approach for dynamic rupture modelling. We introduce a 2D spectral element method (SEM) with an embedded smeared discontinuity representing volumetric fault slip. Our diffuse fault SEM is inspired by the stress-glut method of Andrews, 1999. In our approach, a subdomain in which the tangential stresses are limited by a critical shear strength and an empirical friction law is embedded in a purely elastic domain, resembling classical discrete fault representations. Our approach is implemented on a structured quadrilateral mesh within an SEM framework for elastic wave propagation, with PETSc (Balay et al. 2019) as a linear algebra back-end.

Our method collapses volumetric complexities onto a distribution within a compact support instead of the traditional interface approach, making it a flexible inelastic zone alternative for mesh-independent fault representation in dynamic rupture simulations. We conduct 2D numerical experiments, including a kinematically driven Kostrov-like crack and spontaneous dynamic rupture as defined in SCEC community benchmarks (Harris et al., 2018) of increasing complexity. We extract the spectral response from seismograms at different receivers normal and along the fault. We also analyse the capacity of flexible fault representation by including mesh-independent fault geometries. 

Our approach will allow us to incorporate volumetric failure rheologies in SEM dynamic rupture simulations and is part of the TEAR ERC project (www.tear-erc.eu).

How to cite: Hayek Valencia, J. N., May, D., Pranger, C., and Gabriel, A.-A.: Diffuse thick fault representation in 2D SEM for earthquake dynamic rupture simulations, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-12539, https://doi.org/10.5194/egusphere-egu22-12539, 2022.


Louise Jeandet Ribes et al.

Understanding the mechanical properties of the off-fault medium and its interactions with earthquake rupture is essential for a better understanding of the behavior of fault zones. In this framework, two-dimensional, plane strain models are often used to investigate the interplay between seismic rupture propagation and inelastic deformation in the damage zone. The role of pre-stress conditions for faulting and damage has been studied, in particular the influence of Y, the angle between the largest principal stress and the fault strike. However, in plane strain conditions, the out-of-plane stress is often ignored when setting up the initial stress field, and its influence on dynamic rupture and stress evolution has not been inferred. In this study, we explore the role of the out-of-plane pre-stress for a 2D in-plane model in plane strain conditions. We model a 1D right-lateral, strike-slip vertical fault featuring slip-weakening friction law. We first demonstrate theoretically that if the out- of-plane stress is not considered properly in the initial stress field, pre-stress conditions may not correspond to actual strike-slip faulting. We then investigate how changing the initial stress field can influence the rupture and the stress evolution in the off-fault medium. Our results show that if it does not influence significantly the rupture dynamics, the out-of-plane stress is essential in controlling the evolution of the off-fault medium, especially the localization and extend of areas affected by plastic yielding. Therefore, our results demonstrate the importance of considering properly the initial out-of-plane stress to infer the extend, magnitude and distribution of damage in 2D plane strain simulations with off fault plastic deformation.

How to cite: Jeandet Ribes, L., Thomas, M., and Bhat, H.: Influence of pre-stress conditions in 2D plane strain simulations of a dynamic rupture with off fault damage, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-13390, https://doi.org/10.5194/egusphere-egu22-13390, 2022.


Zihua Niu et al.

Under dynamic perturbations, it has been observed that materials like sedimentary rocks show complex mechanical behaviors. They include the simultaneous dependence of the elastic moduli and attenuation on strain at the same time scale of the perturbations, as well as the conditioning and recovery of the elastic moduli that may happen at time scales that are much larger. The latter cases were recently referred to as non-classical nonlinearity. Aside from laboratory experiments, comparable observations of the non-classical nonlinearity have been made in the field over the past two decades with the development of long-term continuous monitoring of the velocity field inside the Earth using methods such as ambient noise interferometry.


A variety of mathematical models that can potentially quantify the non-classical nonlinearity have already been proposed, e.g., the Damage–Breakage Rheology Model, the Internal Variable Model and the Godunov–Peshkov–Romenski model. However, implementing them in numerical schemes suitable to reproduce nonlinear effects in wave propagation on the local, regional, or global scale is challenging. This can be of interest for constraining a more realistic dynamic rheology for the Earth with the field observations.


In this work, wave propagation in different non-classical nonlinear models is implemented in FEniCS using the discontinuous Galerkin (DG) method in 1D. Behaviors of the different models are systematically studied and quantitatively compared against measurements. This work lays the foundation for an extension to the simulation of 2D/3D wave propagation in the Earth on the large-scale DG simulation frameworks, e.g., SeisSol and ExaHyPE.

How to cite: Niu, Z., Gabriel, A.-A., May, D., Sens-Schönfelder, C., and Igel, H.: A Discontinuous-Galerkin approach to model non-classical nonlinearity observed from lab to global scales, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-10523, https://doi.org/10.5194/egusphere-egu22-10523, 2022.


Thomas Ulrich et al.

The 2018 Mw 7.5 Palu earthquake struck the Sulawesi island, Indonesia, in 2018 and was followed by an unexpected tsunami. Using a physics-based, coupled earthquake-tsunami model, Ulrich et al. (2019) showed that direct earthquake-induced uplift could have sourced the tsunami. The 3D dynamic rupture model of the earthquake captures key observations, including the supershear rupture speed and the deformation pattern derived from satellite data. Stress state and fault conditions were tightly constrained by observations combined with simple static analyses based on Mohr-Coulomb theory of frictional failure and a few trial models. The earthquake scenario predicts a combination of up to 6 m of left-lateral slip and 2 m of normal slip on a straight fault segment dipping 65 degrees beneath Palu Bay.

While most studies (e.g. Bai et al., 2018, Ulrich et al., 2019, Oral et al., 2019) suggest a very early supershear transition, the exact timing of the onset of supershear rupture and the driving mechanism of the supershear transition are elusive. Here we revisit the earthquake dynamic rupture modeling based on new high-resolution near-fault deformation maps derived from correlation of optical satellite data. We vary nucleation radius, fault geometry, and off-fault plasticity parametrization to obtain alternative dynamic rupture scenarios. Specific inputs allow delayed transition to supershear. The obtained scenarios are evaluated based on near-fault damage inference.

Additionally, we revisit the tsunami model, adopting advanced strategies for earthquake-tsunami linking and tsunami modeling. In Ulrich et al. (2019), a one-way linking approach with a shallow water equations solver allowed translating the time-dependent seafloor displacements into a tsunami model with wave amplitudes and periods matching those measured at the Pantoloan wave gauge and inundation that is consistent with field survey data. Such modeling workflow yet neglects tsunami generation complexity, acoustic waves, and dispersion, and only approximates horizontal momentum transfer.  We present a 3D fully coupled earthquake-tsunami model (Krenz et al., 2021), that releases these limitations. This allows us to assess how the standard earthquake-tsunami workflow affects our results, and to revisit our conclusions.

How to cite: Ulrich, T., Marconato, L., Gabriel, A.-A., and Klinger, Y.: Revisiting earthquake-tsunami models of the 2018 Palu events using near-fault high-resolution imaging and 3D fully-coupled earthquake-tsunami modeling, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-5488, https://doi.org/10.5194/egusphere-egu22-5488, 2022.


Yuk Po Bowie Chan et al.

In the West of Northern America, the Cascadia subduction zone that extends over one thousand kilometers has well-documented geological records of megathrust earthquakes. The most recent one occurred in 1700 AD with a moment magnitude of 9. Hence, it has been more than 300 years since the last earthquake, suggesting that Southern Cascadia is mature for the next large earthquake. Estimating future rupture scenarios is therefore crucial for earthquake hazard assessment in the region. Multiple interseismic locking distributions have been proposed for Cascadia. Since each locking model differs from another, it remains unclear how to estimate future rupture extents from interseismic locking distributions. Here, we use 3-D dynamic rupture simulations to investigate the potential rupture segmentation in Cascadia and test the dependency of rupture propagation on hypocenter, especially for the Southern Cascadia. We process the slip deficit distributions from locking models by interpolation and smoothening with a gaussian filter. We then calculate the corresponding stress changes with the assumption that all slip deficits would be released during a coseismic event and derive different initial stress distributions by prescribing constant dynamic stress. For the northern segment, the stress-shadowing (Lindsey et al. 2021) and the viscoelastic (Li et al. 2018) interseismic locking models based respectively on elastic and visco-elastic deformation have similar stress levels, lower than those derived from the Gamma model (Schmalzle et al. 2014). In addition, the Gamma model displays a distinct low-stress gap in the central segment but the stress-shadowing and viscoelastic models show smooth transition stress changes. Since the stress-shadowing and the viscoelastic locking models bear a resemblance, dynamic simulations are then developed based on the initial stress conditions derived from the viscoelastic and the Gamma models by prescribing artificial nucleation zones on the fault plane with varied hypocentre locations. Preliminary results demonstrate three major rupture scenario types - self-arrested, segmented, and full-margin ruptures for both stress models. Given the same conditions, both models indicate that Southern Cascadia with a shorter recurrence interval has a lower potential of growing into a margin-wide rupture compared to the central segment. The southern segment mainly hosts self-arrested and segmented ruptures with Mw ranging from ~6.7 to >7.3. Another finding is strong along-strike variations in stress distribution flavor segmented ruptures while homogeneous stress field promotes margin-wide ruptures. For ruptures initiating from the central segment, several segmented ruptures with Mw 8.14 to larger than 8.25 are observed from the Gamma model but such features are absent in the viscoelastic model. Apart from segmented ruptures, the margin-wide ruptures have amplitudes of ground surface vertical displacement comparable to the subsidence record in the A.D. 1700 megathrust earthquake, particularly for the along-strike variation in the Gamma model.

How to cite: Chan, Y. P. B., Yang, H., and Yao, S.: Estimation of Rupture Scenarios along the Cascadia Megathrust from Interseismic Locking Models, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-12624, https://doi.org/10.5194/egusphere-egu22-12624, 2022.


Mon, 23 May, 17:00–18:30

Chairpersons: Alice-Agnes Gabriel, Jorge N. Hayek, Casper Pranger

Prerna Singh et al.

The response of gravity retaining walls during ground motion is still a challenging field. Recent developments in computational methods have opened the possibility of enhancing the understanding of the non-linear nature of soil-structure systems, e.g., earth pressure thrust acting on the retaining wall, translational and rotational movements, propagation of waves in the soil more realistically and quickly. Till today, Mononobe Okabe (MO) method (pseudo-static) is the most used analytical method because of its simplicity. However, there are many limitations and gives over-conservative results in terms of earth pressure thrust, and many literatures have already justified such a response. Several improved studies are already available, but very few have considered proper soil-structure interaction, real-time input earthquake data (not sinusoidal), and a sufficient number of earthquakes to evaluate the response acting on the wall during dynamic loading.

We seek contribution by analyzing the problem numerically using FE software Plaxis 2D and studying the behavior of retaining wall during seismic loading (range of amax = 0.053g to 1.2g) in terms of acceleration, displacement, rotation, and earth pressure thrust of retaining wall. The main contribution observed is the acceleration was not uniform throughout the medium instead gets amplified up to around 0.6g and later gets attenuated with maximum amplification occurring at the top of the retaining wall followed by the top of backfill soil and base of the wall. The residual displacement and rotation showed an incremental trend with an increase in horizontal seismic coefficient (kh). The earth pressure thrust obtained using numerical analysis was comparatively less than predicted by the MO method.

Keywords: Gravity retaining wall; Acceleration amplification response; Earth pressure thrust; Finite element method; 

How to cite: Singh, P., Bhartiya, P., Chakraborty, T., and Basu, D.: Numerical Advances in Understanding the Behavior of Gravity Retaining Wall during Seismic Motions, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-349, https://doi.org/10.5194/egusphere-egu22-349, 2022.


Yi Zhang et al.

Acoustic wave equations are widely employed in wavefield extrapolation and inversion due to their simplicity compared to the elastic wave equations. In anisotropic media, qP- and qSV-waves are coupled. Multiple acoustic approximations in the vertical transversely isotropic (VTI) media have been proposed during the last decades. A classic way is to set the vertical S-wave velocity zero. As such, the S-wave artefacts still exist, whose amplitude increases with anisotropy. Setting S-wave velocity zero in all propagating directions tackles the issue. However, the higher-order spatial derivatives in the pure qP-wave equation make it hard to solve in the space domain. The spatial derivatives in the denominator of the pure qP-wave equation make the solution by the spatial-domain finite-difference unstable.  In this study, we employed the time-domain pseudospectral method to solve both the classic acoustic wave equation and the pure qP-wave equation in VTI media. Hybrid absorbing boundary conditions are used. Both equations are applied to reverse time migration (RTM) for the anisotropic Marmousi model. The new qP-wave equation outperformed the classic qP-wave equation regarding the computational time. Further work can be extended to waveform inversion with the pure qP-wave equation.

How to cite: Zhang, Y., De Siena, L., and Kaus, B.: Simulation of pure qP-wave in vertical transversely isotropic media, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-5432, https://doi.org/10.5194/egusphere-egu22-5432, 2022.


Hideo Aochi and Kenichi Tsuda

It is widely accepted that the rupture area of earthquake is controlled by fault geometry and the interaction between segments. Besides, many earthquakes do not rupture the whole seismogenic depth but only some limited depth zone. It is not so often to observe that a moderate earthquake such as the 2019 Mw4.9 Le Teil, France, earthquake shows a clear surface rupture and the very shallow rupture area limited at the first 1-2 km depth. Aochi and Tsuda (EGU, 2021) propose the concept that the fault is not uniformly loaded along dip due to the 1D layered structure. Namely, the stress is loaded mainly on the stiff layers, while the soft layers play a role of barrier. We use a boundary element method and a spectral element method for simulating the dynamic rupture propagation and wave propagation. We then demonstrate that the shallowest soft layer can be slipped if the rupture at deeper portion is sufficiently developed. On the other hand, a depth soft layer is difficult to be ruptured, mainly because the absolute stress level is high. In our synthetic scenarios, we compare the ground motions around the fault. In the usual model where the stress is uniformly loaded on all the depths, we observe a strong coherent pulse as the rupture progresses fast to the ground surface. However, we observe more than one pulse in our setting. Such heterogeneous condition along dip should be important to investigate the causality of the seismic asperity and the influence on the resultant near-field ground motion.

How to cite: Aochi, H. and Tsuda, K.: Numerical simulation of dynamic rupture and ground motion on a fault non-uniformly loaded along dip, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-1672, https://doi.org/10.5194/egusphere-egu22-1672, 2022.


Grzegorz Kwiatek and Yehuda Ben-Zion

We investigate theoretical limits to detection of fast and slow seismic events and discuss spatial variations of ground motion expected from an synthetic family of M6 earthquakes at short epicentral distances. The performed analyses are based on synthetic velocity seismograms calculated with the discrete wavenumber method assuming seismic velocities and attenuation properties of the crust in Southern California. The examined source properties include  magnitudes ranging from M -1.0 to M 6.0, static stress drops (0.1-10 MPa), and slow and fast ruptures (0.1-0.9 of shear wave velocity). For the M 6.0 events we also consider variations in rise times producing crack- and pulse-type events and different rupture directivities. We found slow events produce ground motions with considerably lower amplitude than corresponding regular fast earthquakes with the same magnitude, and hence are significantly more difficult to detect. The static stress drop and slip rise time also affect the maximum radiated seismic motion, and thus event detectability. Apart from geometrical factors, the saturation and depletion of seismic ground motion at short epicentral distances stem from radiation pattern, earthquake size (magnitude, stress drop), and rupture directivity. The rupture velocity, rise time and directivity affect significantly the spatial pattern of the ground motions. The results can help optimizing detection of slow and fast dynamic small earthquakes and understand the spatial distribution of ground motion generated by large events.

How to cite: Kwiatek, G. and Ben-Zion, Y.: Detection Limits and Near-Field Ground Motions of Fast and Slow Earthquakes, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-3555, https://doi.org/10.5194/egusphere-egu22-3555, 2022.


Fanny Lehmann et al.

The seismic risk in France, a region of low to moderate seismicity, is of paramount importance given the large number of industrial and nuclear installations. However, the large uncertainties on the geology and the poor knowledge of active faults make the seismic hazard estimation a challenging task. Despite being a promising tool to explore the underlying uncertainties, numerical simulations must be duly calibrated by reproducing specific events.

In this work, we considered the 2019 Mw4.9 earthquake that occurred at Le Teil village in southern France. This event was recorded by 17 stations of 3-component accelerometers, within an area of 50 km around the epicenter (French Accelerometric Network). We used these records to calibrate the numerical simulation. The seismological P- and S-wave speed profiles used result from a 3D weighted average model for Metropolitan France. In addition, the topography was included in the spatial discretization. The uncertainties on dip, strike, and rake angles were explored in order to calibrate the far-field synthetic ground motion model by determining the eigenquakes that efficiently span a large diversity of sources.

A good agreement between synthetic and recorded time histories was found, despite the simplicity of the geological and source model.

How to cite: Lehmann, F., Gatti, F., Bertin, M., and Clouteau, D.: First calibration of the physics-based ground motion model of the 2019 Mw4.9 Le Teil earthquake (France), EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-7558, https://doi.org/10.5194/egusphere-egu22-7558, 2022.


Giulia Sgattoni et al.

Ground motion prediction is one of the main goals in seismic hazard assessment. Empirical ground motion prediction equations may fail to reproduce the complexity of ground shaking in complex 3D media and therefore the use of full waveform modelling is increasingly adopted to model ground shaking. The knowledge of the 3D crustal structure in terms of geometries of the main discontinuities and velocities is fundamental to model wave propagation. However, we often lack detailed geological and geophysical information to build reliable models.

We exploit here a large set composed of high-resolution 2D and 3D seismic data and of about 40 wells with stratigraphic and velocity information, both onshore and offshore, to constrain a 3D crustal velocity model in a sector of the Calabrian accretionary prism (southern Italy). We interpret the main reflection discontinuities and constrain their depth at all available wells in the study area and we use well’s check-shots and velocity data to estimate interval-velocities of the main stratigraphic units. We then combine all depth and velocity information into a regional 3D crustal velocity model of the first 8-10 km. This is subsequently extended to a depth of ~50 km using available regional crustal models to obtain the final model used for ground motion simulation.

We implement our crustal model in the spectral-element code SPECFEM3D_Cartesian to simulate wave propagation in the 3D velocity model honoring surface topography. This allows reconstructing the low-frequency part of the waveforms (up to ~1 Hz), which is then combined with high-frequency seismograms obtained with a stochastic method following the hybrid broadband simulation approach by Graves and Pitarka (2010).

We evaluate the goodness of our model by simulating real earthquakes and comparing simulated and recorded waveforms at the available seismic stations in the area. We compare the results from our 3D model with the ones obtained using a local tomography model and the European crust model EPcrust. The maps of ground motion obtained from the simulated broadband waveforms are then compared with empirical ShakeMaps. These results will also be useful for earthquake scenario calculations, by simulating potential seismic sources identified from structural analysis of geological and seismic data.

How to cite: Sgattoni, G., Molinari, I., Lipparini, L., Faenza, L., and Argnani, A.: Ground-motion simulation in the Calabrian accretionary prism (Southern Italy) using a 3D geologic-based velocity model, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-2676, https://doi.org/10.5194/egusphere-egu22-2676, 2022.


Helena Latečki et al.

The south-eastern part of Adriatic Sea is seismically highly active region where numerous strong events have occurred in historic times. Among these, the most significant is the infamous Great Dubrovnik earthquake of 1667. This event, whose magnitude was estimated to be in the vicinity of Mw 7.0, caused widespread devastation in the whole region. More recently, a large Mw 7.1 event happening in 1979 in Montenegro caused extensive damage along 100 km of coastline, including the area around Dubrovnik. From this it is obvious that the city of Dubrovnik is seismically highly vulnerable and that there is an acute need to better understand possible consequences if an event of such a magnitude would happen today.  
One of the major steps in reducing the seismic risk in any region is to simulate seismic shaking and to evaluate expected ground motion for plausible earthquake scenarios. Therefore, our aim in this work is to create several earthquake scenarios for the city of Dubrovnik and estimate seismically most endangered parts of the region. For that purpose, we first assemble a detailed 3D crustal model which includes information on physical parameters of interest (velocity and density) and which reflects all the important geological features of the studied area. Then, we test whether the model is suitable for simulation by computing and comparing broadband seismograms against the recorded data of several moderate events. We validate the results by assessing the goodness of fit for different metrics describing ground-motion. Next, by combining seismic and geophysical data, we define the geometry of the main active faults and parameters required for the rupture model used in the simulation. We calculate synthetic waveforms on a dense grid and then extract intensity measures to determine the expected ground-motion features of a strong seismic event such was The Great Dubrovnik earthquake.

How to cite: Latečki, H., Sečanj, M., Dasović, I., and Stipčević, J.: Seismic shaking scenarios for city of Dubrovnik, Croatia, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-8291, https://doi.org/10.5194/egusphere-egu22-8291, 2022.


František Gallovič and Ľubica Valentová

Dynamic rupture modeling represents a preferable physics-based approach to strong ground motion simulations. However, its application in a broad frequency range (0-10Hz), interesting for engineering studies, is challenging. The main reason is that relatively simple models with smooth distributions of initial stress and frictional parameters on planar faults result in ground motions with depleted high-frequency content. Several studies suggested that nonplanar rupture surfaces can solve this issue. Nevertheless, fully accounting for rough ruptures typically requires supercomputers, preventing widespread use.

Here we test an efficient approach for the linear slip-weakening friction model on planar fault, based on the Ide and Aochi (2005) multiscale model, with a small-scale fractal distribution of the slip-weakening distance Dc. To intensify the incoherence of the rupture propagation, we also include a variation of the strength and initial stress correlated with Dc. We propose a way to combine the fractal variations of the dynamic parameters with a large-scale dynamic model. The planar fault assumption permits the use of the computationally very fast code FD3D_TSN (Premus et al., 2020). 

We illustrate the approach on a canonical elliptical model with linearly increasing fracture energy (i.e., constant rupture velocity) and the 2016 Mw6.2 Amatrice earthquake smooth rupture model from the dynamic source inversion by Gallovič et al. (2019). We demonstrate that the addition of the small-scale fractal properties results in sustained high-frequency radiation during the rupture propagation and omega-square (apparent) source time functions. The model improves the fit of the recordings of the Amatrice earthquake in the frequency range of 0-10Hz and generates synthetics agreeing with ground motion prediction equations up to 5Hz.

Our FD3D_TSN takes about 5 minutes to simulate the Mw6.2 rupture propagation on a single GPU. Nevertheless, the fractal dynamic model can be easily implemented in any dynamic rupture propagation code. This makes the proposed approach readily applicable in physics-based ground motion predictions for scenario earthquakes in seismic hazard assessment.

How to cite: Gallovič, F. and Valentová, Ľ.: Broadband strong ground motion modeling using planar dynamic rupture model with fractal parameters, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-8698, https://doi.org/10.5194/egusphere-egu22-8698, 2022.


Di Yin et al.

    The Sichuan-Yunnan region is located in the southern part of Chinese north-south seismic belt and has strong seismic activity. The prediction of future strong earthquake activity in this region has always been a research hotspot. In this study, firstly, we established a quasi-three-dimensional finite element elastic model, combined with the regional geological background and GPS observation data. Then, based on the information of 30 M>6.7 historical earthquakes that occurred in the region over the past 100 years, and constrained by the Coulomb-Mohr rupture criterion, we inverted a possible reasonable initial stress field at a specific time. Secondly, we simulated the development process of each historical earthquake and reproduced the 30 events orderly, by comprehensively considering the tectonic stress loading in the seismogenic stage and the stress change in the co-seismic adjustment stage. However, it is worth noting that there were some uncertainties in the numerical simulation process. We used Monte Carlo random experiments to obtain 5000 kinds of different possible initial values, which all can reproduce the development process of historical events. Then we got different current reginal stress values and calculated earthquake risk coefficient. Finally, we used mathematical methods to investigate the current seismic hazard of the different models, and assembled them into a probability distribution map of possible seismic risk in the region. The preliminary result shows that the seismic risk in the rupture zone of historical earthquakes is greatly reduced, which means relatively safe. Mainly due to the stress change caused by the 2008 Wenchuan Ms8.0 earthquake, the seismic probability in the northeastern segment of the Longmenshan fault is as high as 30%. At the junction of the southwestern section of the Longmenshan fault and the Xianshuihe fault zone, the seismic probability is about 15-20%. In addition, near the Longling Ruili fault and the Lancangjiang fault in southwestern Yunnan, the value is about 10-15%. In recent years, small earthquakes have occurred frequently in southwestern Yunnan, and the seismic risk in this area is also worth noting.

How to cite: Yin, D., Dong, P., and Shi, Y.: Numerical analysis of the seismic hazard in Sichuan-Yunnan region, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-10798, https://doi.org/10.5194/egusphere-egu22-10798, 2022.


Celso Alvizuri et al.

We cross validate a numerical solver for wave propagation in 3-D elastic media written on Graphical Processing Units (GPUs) against the semi-analytical solver for earthquakes in axisymmetric media, using seismic full moment tensors as earthquakes sources, variable earthquake source durations, and comparing observed with synthetic seismograms. The GPU-based solver is based on a numerical formulation of elastodynamic wave equation and can capture isotropic and anisotropic media. The algorithm simulates wave propagation in elastic media in three dimensions and at very high spatial and temporal resolution, and can compute entire wavefields within seconds (Alkhimenkov et al., 2021). For example, the multi-GPU code for elastic wave propagation can compute the entire wavefield of a 1000^3 model (1 billion grid cells) in 40 seconds. We achieve a close-to-ideal parallel efficiency (98% and 96%) on weak scaling tests up to 128 GPUs by overlapping MPI communication and computations. Seismic full moment tensors are routinely used to model a range of seismic processes, natural and anthropogenic, including earthquakes (shear slip), volcanic events, explosions, cavity collapses, landslides, etc. The analytical solver is based on a Thompson-Haskell propagator matrix for layered axisymmetric media (Zhu and Rivera, 2002), with moment tensors as seismic sources, with seismic sources at depths 10s of km below the surface and seismic stations at distances over 2000 km, and has been successfully used in various earthquake source studies (e.g. Alvizuri et al 2018). We validate the GPU-based wave propagation solver through numerical experiments in homogeneous and in layered media, and with observed and synthetic seismograms for an M4.6 earthquake in Linthal, Switzerland on 2017-03-06 with seismic stations at distances up to 30 km. The seismograms from the numerical solver match the analytic and observed seismograms (within frequencies 0.02-0.10 Hz). In future work we will apply the solver to study earthquake source generation, wave propagation in anisotropic media, and seismic source determination.

Alkhimenkov, Y., Räss, L., Khakimova, L., Quintal, B., & Podladchikov, Y., 2021. Resolving wave propagation in anisotropic poroelastic media using graphical processing units (CPUs), J. Geophys. Res., 126, doi:10.1029/2020JB021175.
Alvizuri, C., Silwal, V., Krischer, L., & Tape, C., 2018. Estimation of full moment tensors, including uncertainties, for nuclear explosions, volcanic events, and earthquakes, J. Geophys. Res. Solid Earth, 123, 5099–5119, doi:10.1029/2017JB015325.
Zhu, L. & Rivera, L. A., 2002. A note on the dynamic and static displacements from a point source in multilayered media, Geophys. J. Int., 148, 619–627, doi:10.1046/j.1365-246X.2002.01610.x.

How to cite: Alvizuri, C., Alkhimenkov, Y., and Podladchikov, Y.: Cross-validation of a GPU-based wave propagation solver and application to seismic waveform modeling of an M4.6 earthquake in Linthal, Switzerland, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-11592, https://doi.org/10.5194/egusphere-egu22-11592, 2022.


Taufiq Taufiqurrahman et al.

Broadband earthquake ground motion simulations (>1 Hz) are of great interest to seismologists and the earthquake engineering community. The evolution of the earthquake ruptures related to the 2016 Mw 6.2 Amatrice earthquake and the uniquely dense seismological recordings provide an opportunity to understand better the processes controlling earthquake dynamics, strong ground motion, and the relation between earthquakes. We here propose a novel approach to design data-driven broadband (up to 5 Hz) dynamic rupture scenarios from 0.5-1 Hz Bayesian dynamic finite-fault inversion (Gallovič et al., 2019). We analyze the effects of enhancing the best-fitting smooth dynamic source inversion result by subsequent adding of complexity such as non-planar fault geometry (i.e., fault listricity and surface roughness), topography, inelastic off-fault rheology, and visco-elastic attenuation. We utilize the open-source software package SeisSol (www.seissol.org), suited explicitly for incorporating such geometrical complexities and high-resolution simulations performed on modern supercomputers. The obtained scenarios reproduce synthetics resembling the observations in terms of velocity and accelerations waveforms and Fourier-amplitude-spectra (FAS) up to 5 Hz. The simulated peak ground velocity (PGV) maps show de-amplification of ground motion amplitudes on the foot-wall and amplification on the hanging-wall as a consequence of the wave-focusing effect caused by the listric fault curvature. This effect is seen mainly for distances up to 10 km from the fault. Our study suggests that the complexity of the earthquake source should not be neglected for the seismic hazard assessment for regions adjacent to active faults.

How to cite: Taufiqurrahman, T., Gabriel, A.-A., Ulrich, T., Valentová, L., and Gallovič, F.: Broadband Dynamic Rupture and Ground Motion Simulations (up to 5 Hz) of the 2016 Mw 6.2 Amatrice, Italy Earthquake, EGU General Assembly 2022, Vienna, Austria, 23–27 May 2022, EGU22-11610, https://doi.org/10.5194/egusphere-egu22-11610, 2022.