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.

Convener: Alice-Agnes GabrielECSECS | Co-conveners: Jean-Paul Ampuero, Hideo Aochi
vPICO presentations
| Mon, 26 Apr, 15:30–17:00 (CEST)

vPICO presentations: Mon, 26 Apr

Chairpersons: Alice-Agnes Gabriel, Jean-Paul Ampuero, Hideo Aochi
František Gallovič, Jiří Zahradník, Vladimír Plicka, Efthimios Sokos, Christos Evangelidis, Ioannis Fountoulakis, and Fatih Turhan

The 2020 Mw 6.8 Elazığ (Sivrice) earthquake occurred on the Pütürge segment of the East Anatolian Fault zone. This strike-slip segment is situated between strong earthquakes that happened 100–150 years ago, and, since that time, the segment remained with eight Mw5-6 events, but with no Mw 6+. We relocate the mainshock and aftershock sequence and infer basic characteristics of the event using the ISOLA multiple point source approach and backpropagation of S-waveforms from local strong-motion recordings. Together with clear secondary P wave onsets identified in the recordings, the results suggest complex rupture propagation with reversal of the rupture propagation. We apply a recently developed Bayesian dynamic source inversion with slip-weakening friction and spatially inhomogeneous stress and friction parameters to gain better insight into the rupture process. Using high-quality near field recordings in the low-frequency range (<0.3Hz), we obtain a complex dynamic rupture model explaining the weak rupture initiation followed by a cascade of at least three rupture episodes, including the rupture reversal. The dynamic model explains significant features of the recordings even in a broader frequency range interesting for seismic engineering applications (<2.5Hz), e.g., a directivity pulse associated with rupturing the event’s strongest asperity with 4 m of slip and local stress drop of 40 MPa. We show that by reducing the initial stress in the top 10 km by 10%, the rupture fails to develop into the larger event, finishing as an Mw 5.8 earthquake. Considering the latter experiment corresponds to an earlier state of the fault in the seismic cycle, we hypothesize that the interseismic M5+ events on the Pütürge segment were undeveloped rudiments of potentially large events. Thus, the fault seems to have been ready for the Mw6.8 earthquake only by the time of the earthquake occurrence in 2020. This suggests that at the time of the Elazığ earthquake initiation, the final Mw6.8 magnitude was not determined, making it a treacherous case for early warning systems.

How to cite: Gallovič, F., Zahradník, J., Plicka, V., Sokos, E., Evangelidis, C., Fountoulakis, I., and Turhan, F.: Complex rupture dynamics of the 2020 Mw 6.8 Elazığ (Sivrice) earthquake, Turkey, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-14875,, 2021.

Hoby Razafindrakoto, Fabrice Cotton, Dino Bindi, Marco Pilz, and Robert Graves

Over the past decade, there is growing consensus that physics-based simulations can be utilized in engineering applications. However, for the simulations to be accepted, they need to be calibrated and validated. This study presents the results of ground motion simulation calibration and validation using earthquakes that occurred in the Upper Rhine Graben with a modified version of the Graves-Pitarka (GP) hybrid ground-motion simulation methodology implemented on the Southern California Earthquake Center Broadband Platform, which uses an improved high-frequency computation. To calibrate the HF simulation, we take advantage of the growth of seismological data (including weak motions) in the region and the ability to evaluate critical seismic parameters such as anelastic attenuation, stress drop, and site effects through spectral decomposition methods (separate site-source-propagation from the datasets).  Hence in the simulation, the adopted anelastic attenuation and stress parameter are defined based on the spectral decomposition results. The additional modification of the standard GP method is the incorporation of compressional wave in the HF motion.

Results are compared with observations and simulations from the unmodified GP approach; we also use a range of ground motion intensity measures as summary statistics. We found that in general, the modification in the HF part (e.g., incorporation of compressional waves) was necessary to improve the fit with observations. Our findings also validate the fact that parameters from the spectral decomposition are giving well-calibrated time-histories (in terms of frequency and amplitude) when used as input parameters of the broadband simulations. The findings in this study support the incorporation of scenario-based ground motion simulations for use in the characterization of seismic hazard and other engineering applications. For simulation of future earthquakes, instead of using event-specific stress-drops, we use the average stress-drops taken from the distribution of the stress drops derived from spectral decomposition.

How to cite: Razafindrakoto, H., Cotton, F., Bindi, D., Pilz, M., and Graves, R.: Embedding Spectral Decomposition Results in Broadband Simulation: Application to Rhine Graben Area, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8882,, 2021.

Jayalakshmi Sivasubramonian and Paul Martin Mai

We analyze the effect of earthquake source parameters on ground-motion variability based on near-field wavefield simulations for large earthquakes. We quantify residuals in simulated ground motion intensities with respect to observed records, the associated variabilities are then quantified with respect to source-to-site distance and azimuth. Additionally, we compute the variabilities due to complexities in rupture models by considering variations in hypocenter location and slip distribution that are implemented a new Pseudo-Dynamic (PD) source parameterization.

In this study, we consider two past events – the Mw 6.9 Iwate Miyagi Earthquake (2008), Japan, and the Mw 6.5 Imperial Valley Earthquake, California (1979). Assuming for each case a 1D velocity structure, we first generate ensembles of rupture models using the pseudo-dynamic approach of Guatteri (2004), by assuming different hypocenter and asperities locations (Mai and Beroza, 2002, Mai et al., 2005; Thingbaijam and Mai, 2016). In order to efficiently include variations in high-frequency radiation, we adopt a PD parameterization for rupture velocity and rise time distribution in our rupture model generator. Overall, we generate a database of rupture models with 50 scenarios for each source parameterization. Synthetic near-field waveforms (0.1-2.5Hz) are computed out to Joyner-Boore distances Rjb ~ 150km using a discrete-wavenumber finite-element method (Olson et al., 1984). Our results show that ground-motion variability is most sensitive to hypocenter locations on the fault plane. We also find that locations of asperities do not alter waveforms significantly for a given hypocenter, rupture velocity and rise time distribution. We compare the scenario-event simulated ground motions with simulations that use the rupture models from the SRCMOD database (Mai and Thingbaijam, 2014), and find that the PD method is capable of reducing the ground motion variability at high frequencies. The PD models are calibrated by comparing the mean residuals with the residuals from SRCMOD models. We present the variability due to each source parameterization as a function of Joyner-Boore distance and azimuth at different natural period.

How to cite: Sivasubramonian, J. and Mai, P. M.: Spatial Variability of Near-field Ground Motions from Pseudo-Dynamic Rupture Simulations, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9463,, 2021.

Jagdish Chandra Vyas, Martin Galis, and Paul Martin Mai

Geological observations show variations in fault-surface topography not only at large scale (segmentation) but also at small scale (roughness). These geometrical complexities strongly affect the stress distribution and frictional strength of the fault, and therefore control the earthquake rupture process and resulting ground-shaking. Previous studies examined fault-segmentation effects on ground-shaking, but our understanding of fault-roughness effects on seismic wavefield radiation and earthquake ground-motion is still limited.  

In this study we examine the effects of fault roughness on ground-shaking variability as a function of distance based on 3D dynamic rupture simulations. We consider linear slip-weakening friction, variations of fault-roughness parametrizations, and alternative nucleation positions (unilateral and bilateral ruptures). We use generalized finite difference method to compute synthetic waveforms (max. resolved frequency 5.75 Hz) at numerous surface sites  to carry out statistical analysis.  

Our simulations reveal that ground-motion variability from unilateral ruptures is almost independent of  distance from the fault, with comparable or higher values than estimates from ground-motion prediction equations (e.g., Boore and Atkinson, 2008; Campbell and Bozornia, 2008). However, ground-motion variability from bilateral ruptures decreases with increasing distance, in contrast to previous studies (e.g., Imtiaz et. al., 2015) who observe an increasing trend with distance. Ground-shaking variability from unilateral ruptures is higher than for bilateral ruptures, a feature due to intricate seismic radiation patterns related to fault roughness and hypocenter location. Moreover, ground-shaking variability for rougher faults is lower than for smoother faults. As fault roughness increases the difference in ground-shaking variabilities between unilateral and bilateral ruptures increases. In summary, our simulations help develop a fundamental understanding of ground-motion variability at high frequencies (~ 6 Hz) due small-scale geometrical fault-surface variations.

How to cite: Vyas, J. C., Galis, M., and Mai, P. M.: High-frequency ground-motion variability for rough-fault ruptures  , EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-11998,, 2021.

Nesrin Yenihayat, Eser Çaktı, and Karin Şeşetyan

One of the major earthquakes that resulted in intense damages in Istanbul and its neighborhoods took place on 10 July 1894. The 1894 earthquake resulted in 474 losses of life and 482 injuries. Around 21,000 dwellings were damaged, which is a number that corresponds to 1/7 of the total dwellings of the city at that time. Without any doubt, the exact loss of life was higher. Because of the censorship, the exact loss numbers remained unknown. There is still no consensus about its magnitude, epicentral location, and rupture of length. Even though the hardness of studying with historical records due to their uncertainties and discrepancies, researchers should enlighten the source parameters of the historical earthquakes to minimize the effect of future disasters especially for the cities located close to the most active fault lines as Istanbul. The main target of this study is to enlighten possible source properties of the 1894 earthquake with the help of observed damage distribution and stochastic ground motion simulations. In this paper, stochastic based ground motion scenarios will be performed for the 10 July 1894 Istanbul earthquake, using a finite fault simulation approach with a dynamic corner frequency and the results will be compared with our intensity map obtained from observed damage distributions. To do this, in the first step, obtained damage information from various sources has been presented, evaluated, and interpreted. Secondly, we prepared an intensity map associated with the 1894 earthquake based on macro-seismic information, and damage analysis and classification. For generating ground motions with a stochastic finite fault simulation approach, the EXSIM 2012 software has been used. Using EXSIM, several scenarios are modeled with different source, path, and site parameters. Initial source properties have been obtained from findings of our previous study on the simulation of the 26 September 2019 Silivri (Istanbul) earthquake with Mw 5.8. With the comparison of spatial distributions of the ground motion intensity parameters to the obtained damage and intensity maps, we estimate the optimum location and source parameters of the 1894 Earthquake.

How to cite: Yenihayat, N., Çaktı, E., and Şeşetyan, K.: 10 July 1894 Istanbul Earthquake: Comparing Damages and Ground Motion Simulations, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-16183,, 2021.

María del Puy Papí Isaba, Christa Hammerl, Maurizio Mattesini, and Vicenta María Elisa Buforn Peiró

Tyrol is one of the provinces with the highest seismicity in Austria. Most of the stronger historical earthquakes occurred around Innsbruck and Hall in Tirol (1572, 1670, 1689).

Within the framework of the project[1] “Historical and recent earthquake activity in Tyrol - sources, data, seismological analysis”, a study was carried out from 2014-2020, which mainly deals with historical earthquakes in Tyrol up to 1900 but also in detail with damaging earthquakes in Tyrol in the 20th century. The project’s purpose was to create a new earthquake catalog for Tyrol, which for the first time also includes Macroseismic/Intensity Data Points (M/IDPs).

An essential aspect of this study is that the sources and literature references used for all Tyrolean earthquakes up to 1900 are largely documented. Furthermore, selected damaging earthquakes of the 20th century are reported in detail. Numerous Tyrolean archives, such as the Tyrolean Provincial Archives, and the City Archives of Hall in Tyrol, were searched for contemporary earthquake sources. Likewise, the seismic archive of the Austrian Seismological Service at ZAMG (Zentralanstalt für Meteorologie und Geodynamik) contains a wealth of valuable recent information, such as the questionnaires on earthquakes of the entire 20th century.

The very time-consuming research and documentation are followed by the conversion of the written information into earthquake parameters. Briefly outlined, this comprises the following working steps: Interpretation of the sources, assignment of geographical coordinates to the pieces of evidence, evaluation of the intensity according to the European Macroseismic Scale (EMS-98), (re)calculation of the focal parameters of all damaging earthquakes and numerous newly found earthquakes.

The latter is the content of this presentation, namely to (re-)evaluate the focal parameters for historical and recent earthquakes in Tyrol for the first time using the intensity prediction equations (IPE) with the Grid Search (GS) technique. GS has been widely used in many Machine Learning types of research when it comes to hyperparameter optimization, which in this study corresponds to the earthquake focal parameters.

We used IDPs whose intensities were mostly assessed from contemporary historical sources, such as annals, chronicles, questionnaires, newspapers, etc.

A total of 1750 M/IDPs for 35 damaging earthquakes from the Austrian Earthquake Catalogue (AEC2020) could be determined based on the historical sources. The focal parameters for these earthquakes were reevaluated by means of the IPE and GS. 

Likewise, 726 new M/IDPs from a total of 154 non-damaging earthquakes not yet included in the AEC2020 were determined. For 38 of them, it was possible to calculate new sets of focal parameters.

Problems encountered, accuracy, and error of the results will be introduced in the presentation. 


[1]The project was funded by TIWAG-Tiroler Wasserkraft AG, ASFINAG Alpenstraße GmbH, Fachgruppe der Seilbahnen Tirol, Verbund Hydro Power GmbH, Amt der Tiroler Landesregierung - Abteilung Allgemeine Bauangelegenheiten, Landesgeologi, ÖBB Infrastruktur AG and ZAMG.

How to cite: Papí Isaba, M. P., Hammerl, C., Mattesini, M., and Buforn Peiró, V. M. E.: Can Machine Learning help us in [re]assessing historical earthquakes?, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-597,, 2021.

Viktoriya Yarushina and Alexander Minakov

The microseismic events can often be characterized by a complex non-double couple source mechanism. Recent laboratory studies recording the acoustic emission during rock deformation help connecting the components of the seismic moment tensor with the failure process. In this complementary contribution, we offer a mathematical model which can clarify these connections. We derive the seismic moment tensor based on classical continuum mechanics and plasticity theory. The moment tensor density can be represented by the product of elastic stiffness tensor and the plastic strain tensor. This representation of seismic sources has several useful properties: i) it accounts for incipient faulting as a microseismicity source mechanism, ii) it does not require a pre-defined fracture geometry, iii) it accounts for both shear and volumetric source mechanisms, iv) it is valid for general heterogeneous and anisotropic rocks, and v) it is consistent with elasto-plastic geomechanical simulators. We illustrate the new approach using 2D numerical examples of seismicity associated with cylindrical openings, analogous to wellbore, tunnel or fluid-rich conduit, and provide a simple analytic expression of the moment density tensor. We compare our simulation results with previously published data from laboratory and field experiments.  We consider three special cases corresponding to "dry" isotropic rocks, "dry" transversely isotropic rocks and "wet" isotropic rocks. The model highlights theoretical links between stress state, geomechanical parameters and conventional representations of the moment tensor such as Hudson source type parameters. 

How to cite: Yarushina, V. and Minakov, A.: Can plasticity explain microseismic source mechanisms?, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-2573,, 2021.

Huihui Weng, Jean-Paul Ampuero, and Loes Buijze

The induced seismicity in the Groningen gas field, The Netherlands, has led to intense public concerns and comprehensive investigations. One of the main challenges for assessing future seismic hazard in the Groningen gas field is to estimate the maximum possible earthquake magnitude (Mmax) that could be induced by gas extraction. Previous methods are strongly rooted in empirical and statistical approaches that are inherently limited by the scarcity of data. Here, we combine a physics-based dynamic rupture model based on the 3D theory of fracture mechanics with field-based and lab-based constraints to estimate Mmax in the Groningen gas field. If earthquakes in the reservoir have a rupture depth extension constrained by the reservoir thickness, the largest earthquakes should develop a large aspect ratio (longer horizontally than vertically). The model is thus an extension of the 3D theoretical rupture model on long faults with uniform stress and strength developed by Weng & Ampuero (2019), in which we have incorporated spatial heterogeneities, such as along-strike variable fault width, depth-dependent initial stresses and friction properties. The essential parameters that control rupture propagation and earthquake magnitude are the stored elastic energy and the fracture energy. Our method requires estimates of the stored elastic energy on reservoir faults as a result of the stresses induced by differential reservoir compaction during depletion. The fracture energy is constrained by laboratory experiments and theoretical frictional models. Coupling physics-based rupture models with field and lab observations provides an estimate of Mmax in the Groningen gas field and serves as a practical step toward physics-based seismic hazard assessment for other gas fields in the world.



Weng, H. and J. P. Ampuero (2019). "The Dynamics of Elongated Earthquake Ruptures." Journal of Geophysical Research: Solid Earth.

How to cite: Weng, H., Ampuero, J.-P., and Buijze, L.: Physics-Based Estimates of the Maximum Magnitude of Induced Earthquakes in the Groningen Gas Field, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-6144,, 2021.

Waheed Gbenga Akande, Quan Gan, David G. Cornwell, and Luca De Siena

Modelling volcanic processes at active volcanoes often requires a multidisciplinary approach, which adequately describes the complex and ever-dynamic nature of volcanic unrests. Campi Flegrei caldera (southern Italy) is an ideal laboratory where numerical modelling of injection-induced seismicity could be tested to match the observed seismicity. In the current study, thermal-hydraulic-mechanical (THM) effects of hot-water (fluid) injections were investigated to ascertain whether the observed seismicity (past and ongoing seismic swarms) could be quantitatively reproduced and modelled in isothermal or non-isothermal approximations. Fluid-flow modelling was carried out using a coupled TOUGHREACT-FLAC3D approach to simulate THM effects of fluid injections in a capped reservoir, where the sealing formation serves as a geological interface between supercritical reservoir and fractured shallow layers of the caldera. Results from previous seismic, deformation, tomographic and rock physics studies were used to constrain the model for realistic volcano modelling. The results indicated that fluid injections generated overpressure beneath the caprock and subjected it to different stress regimes at its top and bottom, and this prompted deformation. Thus, caprock deformation, triggered by injection-induced basal compressional forces and top extensional fractures, is a critical factor determining the required timing for pressure build-up and fault reactivation, and magnitudes of seismicity. Higher fluid injection rates and temperature contrasts, heterogeneity due to fault and its contrast with the host rock, and caprock hydraulic properties were among the identified secondary factors modulating fault reactivation and seismicity. Simulation results revealed that seismicity can be better modelled in isothermal (HM) approximations. A comparative study of the THM-modelled seismicity and 4-month-long (August 5th to December 5th, 2019) seismic monitoring data recorded at the Osservatorio Vesuviano showed that our model reproduced the magnitudes and depths (~2.5 Ms within 2 km) at the onset of the ongoing unrests on October 5th, 2019. However, the model could not adequately reproduce the highest magnitude (3.3 Ms at 2.57 km) seismicity on April 26th, 2020 observed since 1984 major unrests.


How to cite: Akande, W. G., Gan, Q., Cornwell, D. G., and De Siena, L.: Numerical Modelling of Injection-induced Seismicity at Campi Flegrei Caldera, Southern Italy, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-9064,, 2021.

Lise Alalouf and Yajing Liu

Subduction zones are where the largest earthquakes occur. In the past few decades, scientists have also discovered the presence of episodic aseismic slip, including slow slip events (SSEs), along most of the subduction zones. However, it is still unclear how these SSEs can influence megathrust earthquake ruptures. The Costa Rica subduction zone is a particularly interesting area because a SSE was recorded 6 months before the 2012 Mw7.6 earthquake in the Nicoya Peninsula, suggesting a potential stress transfer from the SSE to the earthquake slip zone. SSEs beneath the Nicoya Peninsula were also recorded both updip and downdip the seismogenic zone, making it a unique area to study the complex interaction between SSEs and earthquakes.

As most of the shallow SSEs were recorded around the Nicoya Peninsula, we chose to start using a 1D planar fault embedded in a homogeneous elastic half-space, with different dipping angles following several geometric profiles of the subduction fault beneath the Nicoya Peninsula section of the Costa Rica margin. This 1D modelling study allows us to better investigate the interaction between shallow and deep SSEs and megathrust earthquakes with high numerical resolution and relatively short computation time. The model provides information on the long-term seismic history by reproducing the different stages of the seismic cycle (interseismic slip, shallow and deep episodic slow slip, and coseismic slip).

We study the influence of the variation of numerical parameters and frictional properties on the recurrence interval, maximum slip velocity and cumulative slip of SSEs (both shallow and deep) and earthquakes and their interaction with each other. We then compare our results with GPS and seismic observations (i.e. cumulative slip, characteristic duration, moment rate, depth and size of the rupture, equivalent magnitude) to identify an optimal set of model parameters to understand the interaction between various modes of subduction fault deformation.

How to cite: Alalouf, L. and Liu, Y.: Modeling the Interaction between Slow Slip Events and Earthquake Ruptures in the Nicoya Peninsula, Costa Rica. , EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-6324,, 2021.

Hideo Aochi and Kenichi Tsuda

Dynamic rupture simulation of an earthquake mostly aims at a characteristic event, which may rupture the entire seismogenic zone of a fault system, perhaps reaching the ground surface. However, hazardous earthquakes sometimes occur along a part of the depths of a fault. Many questions arise why only this particular depth does rupture and whether the surrounding part remains hazardous. Previously, Aochi (GJI, 2018) has considered a depth-dependent stress accumulation for emphasizing the difference of reverse and normal faults under the hypothesis that stress is sufficiently and uniformly charged at all depths. We probably need to revise this hypothesis and the partially charged fault along depth would be more suitable for explaining the given question. By developing the previous simulations by Aochi (GJI, 2018), we carry out numerical simulations for demonstrating the importance of the depth-dependent stress accumulation.  

How to cite: Aochi, H. and Tsuda, K.: On the depth-dependent stress accumulation for earthquake generation process, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8682,, 2021.

Thomas Ulrich, Alice-Agnes Gabriel, Yann Klinger, Jean-Paul Ampuero, Percy Galvez, and Bo Li

The Dead-Sea Transform fault system, a 1200 km-long strike-slip fault forming the tectonic boundary between the African Plate and the Arabian Plate, poses a major seismic hazard to the eastern Mediterranean region. The Gulf of Aqaba, which terminates the Dead Sea fault system to the South, results from a succession of pull-apart basins along the Dead-Sea Transform fault system. The complexity of the fault system in the Gulf has been recently evidenced by Ribot et al. (2020), who compiled a detailed map of its fault traces, based on a new multibeam bathymetric survey of the Gulf. Part of the Gulf of Aqaba was ruptured by an Mw 7.3 earthquake in 1995. Teleseismic data analysis suggests that it may have been a multi-segment rupture (Klinger et al., 1999). This event occurred offshore, in a poorly instrumented region, and therefore the exact sequence of faults that ruptured is not precisely known. The detailed fault mapping of Ribot et al. (2020) offers a fresh view of this earthquake. In particular, it identifies many oblique faults between the major strike-slip faults, which may have linked these segments.

Relying on this new dataset, on a new back-projection study, and on 3D dynamic rupture modeling with SeisSol (, we revisit the 1995 Aqaba earthquake. Using back projection, we identify 2 strong radiators, which we associate with 2 step-overs. Using 3D dynamic rupture modeling, we propose scenarios of the 1995 earthquake, compatible with the various dataset available. Our modeling allows constraining the regional state of stress in the region, acknowledging transtension, offers constraints on the nucleation location and confirms the role of the oblique faults in propagating the rupture to the North. It offers new constraints on the regional seismic hazard, in particular on the expected maximum moment magnitude.

Finally, we explore the dynamics of the Gulf of Aqaba fault system using earthquake cycle modeling. For that purpose, we rely on QDYN (, a boundary element software, which simulates earthquake cycles under the quasi-dynamic approximation on faults governed by rate-and-state friction and embedded in elastic media. We inform our parameterization of the earthquake cycle modeling using the previously described datasets and modeling results. Recently Galvez et al. (2020) demonstrated the capability of the method to model the dynamics of complex fault system in 3D. Here new code developments are required to adapt the method to the Gulf of Aqaba fault system, e.g. to allow accounting for normal stress changes and for variations in the fault rake.

Overall, we aim to better understand how large earthquakes may nucleate, propagate, and interact across a complex transform fault network. Our findings, e.g. on fault segmentation or the conditions that promote larger earthquakes, will have important implications for other large strike-slip fault systems worldwide.

How to cite: Ulrich, T., Gabriel, A.-A., Klinger, Y., Ampuero, J.-P., Galvez, P., and Li, B.: Exploring the dynamics of the Dead-Sea Transform fault using data-integrated numerical models, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-7496,, 2021.

Meng Li, Casper Pranger, and Ylona van Dinther

Numerical models are well-suited to overcome limited spatial-temporal observations to understand earthquake sequences, which is fundamental to ultimately better assess seismic hazard. However, high-resolution numerical models in 3D are computationally time and memory consuming. This is not optimal if the aspects of lateral or depth variations within the results are not needed to answer a particular objective. In this study we quantify and summarize the limitations and advantages for simulating earthquake sequences in all spatial dimensions.


We simulate earthquake sequences on a strike-slip fault with rate-and-state friction from 0D to 3D using both quasi-dynamic and fully dynamic approaches. This cross-dimensional comparison is facilitated by our newly developed, flexible code library Garnet, which adopts a finite difference method with a fully staggered grid. We have validated our models using problems BP1-QD & FD and BP4-QD & FD of the SEAS (Sequences of Earthquakes and Aseismic Slip) benchmarks from the Southern California Earthquake Center.


Our results demonstrate that lower-dimensional/quasi-dynamic models are qualitatively similar in terms of earthquake cycle characteristics to their higher-dimensional/fully-dynamic counterparts, while they could be hundreds to millions times faster at the same time. Quantitatively, we observe that certain earthquake parameters, such as stress drop and fracture energy release, can be accurately reproduced in each of these simpler models as well. However, higher dimensional models generally produce lower maximum slip velocities and hence longer coseismic durations. This is mainly due to lower rupture speeds, which result from increased energy consumption along added rupture front directions. In the long term, higher dimensional models produce shorter recurrence interval and hence smaller total slip, which is mainly caused by the higher interseismic stress loading rate inside the nucleation zone. The same trend is also observed when comparing quasi-dynamic models to fully dynamic ones. We extend a theoretical calculation that to first order approximates the aforementioned physical observables in 3D to all other dimensions. These theoretical considerations confirm the same trend as what is observed for stress drop, recurrence interval and total slip across dimensions. These findings on similarities and differences of different dimensional models and a corresponding quantification of computational efficiency can guide model design and facilitate result interpretation in future studies.

How to cite: Li, M., Pranger, C., and van Dinther, Y.: Characteristics of earthquake sequences: comparison from 0D to 3D, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-8833,, 2021.

Eyup Sopaci and Atilla Arda Özacar

The 30 October 2020 Samos Earthquake (Mw=7.0) ruptured a north-dipping offshore normal fault, north of the Samos Island with an extensional mechanism. Aftershocks mainly occurred at the western and eastern ends of the rupture plane in agreement with the Coulomb static stress changes. Mechanism of aftershocks located west of the rupture area supported activation of the neighboring strike-slip fault almost instantly. In addition, a seismic cluster including events with magnitudes reaching close to 4 has emerged fifty hours later at the SE side of Samos Island. This off-plane cluster displays a clear example of delayed seismic triggering that produced small magnitude earthquakes at nearby active faults. In this study, numerical simulations are conducted using rate-and-state friction dependent quasi-static&full-dynamic spring slider model with shear-normal stress coupling to mimic the instant and delayed seismic triggering observed after this event. Coulomb static stress changes and seismic waveforms recorded at nearby strong-motion stations are used as static and dynamic triggers during simulations. According to our results, earthquakes with Mw<3.5 can be triggered almost instantly at the rupture edge and failure time of earthquakes with Mw>3.5 advances for both strike-slip and normal faults which may explain the delayed triggering observed SE of Samos Island. Moreover, simulations revealed that the shear-normal stress coupling increases the triggering potential.

How to cite: Sopaci, E. and Özacar, A. A.: Simulation of Instant and Delayed Seismic Triggering Observed After the 30 October 2020 Samos Earthquake at Nearby Faults, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-13090,, 2021.

Ruth Harris, Michael Barall, David Ponce, Diane Moore, Russell Graymer, David Lockner, Carolyn Morrow, Gareth Funning, Christos Kyriakopoulos, and Donna Eberhart-Phillips

The Rodgers Creek-Hayward-Calaveras-Northern Calaveras fault system in California dominates the hazard posed by active faults in the San Francisco Bay Area. Given that this fault system runs through a densely populated area, a large earthquake in this region is likely to affect millions of people. This study produced scenarios of large earthquakes in this fault system, using spontaneous (dynamic) rupture simulations. These types of physics-based computational simulations require information about the 3D fault geometry, physical rock properties, fault friction, and initial stress conditions. In terms of fault geometry, the well-connected multi-fault system includes the Hayward fault, at its southern end the Central and Northern Calaveras faults, and at its northern end the Rodgers Creek fault. Geodetic investigations of the fault system’s slip-rate pattern provide images of where the fault surfaces at depth are creeping or locked interseismically, and this helped us choose appropriate initial stress conditions for our simulations. A 3D geologic model of the fault system provides the 3D rock units and fault structure at depth, while field samples from rocks collected at Earth’s surface provide frictional parameters. We used this suite of information to investigate the behavior of large earthquake ruptures nucleating at various positions along this partially creeping fault system. We found that large earthquakes starting on the Hayward fault or on the Rodgers Creek fault may be slowed, stopped, or unaffected in their progress, depending on how much energy is released by the creeping regions of the Hayward and Central Calaveras faults during the time between large earthquakes. Large earthquakes starting on either the Hayward fault or the Rodgers Creek faults will likely not rupture the Northern Calaveras fault, and large earthquakes starting on either the Northern Calaveras fault or the Central Calaveras fault will likely remain confined to those fault segments.

How to cite: Harris, R., Barall, M., Ponce, D., Moore, D., Graymer, R., Lockner, D., Morrow, C., Funning, G., Kyriakopoulos, C., and Eberhart-Phillips, D.: Dynamic Rupture Scenarios of Large Earthquakes on the Rodgers Creek-Hayward-Calaveras-Northern Calaveras Fault System, California, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-13830,, 2021.

P. Martin Mai, Jagdish Vyas, Alice-Agnes Gabriel, and Thomas Ulrich

Frictional heat generated in the fault core during earthquake rupture can raise the fluid pressure in the slip zone. Such increase of fluid pressure decreases the effective normal stress and thereby lowers the frictional strength of the fault. Therefore, thermal pressurization (TP) of pore fluid affects earthquake rupture processes including nucleation, propagation, and arrest. While the effects of pore pressure and fluid flow rate on dynamic weakening of faults are qualitatively understood, a detailed analysis of how TP affects  earthquake rupture parameters is needed to further deepen our understanding. 

In this study, we investigate the role of two key TP parameters -- hydraulic diffusivity and shear-zone half-width -- earthquake dynamics and kinematic source properties (slip, peak slip-rate, rupture speed and rise time). We conduct  a suite of 3D dynamic rupture simulations applying a rate-and-state dependent friction law (with strong velocity weakening) coupled with thermal-pressurization of pore fluids. Simulations are carried out with the open source software SeisSol ( The temporal evolution of rupture parameters over ~1’000 randomly  distributed on-fault receivers is statistically analyzed in terms of  mean variations of rupture parameters and correlations among rupture parameters. 

Our simulations reveal that mean slip decreases with increasing hydraulic diffusivity, whereas mean peak slip-rate and rupture speed remain nearly constant. On the other hand, we observe only a slight decrease of mean slip with increasing shear-zone half-width, whereas mean peak slip-rate and rupture speed show clear decrease. The faster diffusion of pore pressure as hydraulic diffusivity increases promotes faster increase of the effective normal stress (and fault strength) behind the main rupture front, reducing the rise time and, therefore also affecting mean slip. An increase in shear-zone half- width represents a heat source distributed over larger fault normal distance causing a second-order effect on mean slip. Additionally, our simulations reveal correlations among rupture parameters: 1) slip has weak negative correlation with peak slip-rate and negligible correlation with rupture speed, but a positive correlation with rise time, 2) peak slip-rate has a strong positive correlation with rupture speed, but a strong negative correlation with rise time, 3) rupture speed has strong negative correlation with rise time. We observe little or negligible effects of variations of hydraulic diffusivity and shear-zone half- width on the correlations between rupture parameters. Overall, our study builds a fundamental understanding on how thermal pressurization of pore fluids affects dynamic and thereby kinematic earthquake rupture properties. Our findings are thus important for the earthquake source modeling community, and particularly, for assessing seismic hazard due to induced events in geo-reservoirs.

How to cite: Mai, P. M., Vyas, J., Gabriel, A.-A., and Ulrich, T.: Earthquake rupture properties in presence of thermal -pressurization of pore fluids, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15202,, 2021.

Jorge Nicolas Hayek Valencia, Dave A. May, and Alice-Agnes Gabriel

Faults in earthquake rupture dynamic simulations are typically treated as infinitesimally thin planes with distinct on- versus off-fault rheologies. These faults are prescribed and can be explicitly accounted for with hexahedral or unstructured tetrahedral meshing approaches.  
We present a diffuse interface alternative to dynamic rupture modelling on non-mesh aligned faults and, by design, permits modelling of non-planar faults and time-dependent fault geometries. We use se2dr, a spectral finite element (continuous Galerkin) method with a non-mesh aligned embedded diffuse discontinuity for dynamic rupture simulations.

Natural fault systems are characterised by fault zone complexity, e.g. the frictional strength and spatio-temporal slip localisation may change drastically from the outer damage zone to the fault core. Complex volumetric failure patterns are observed in well-recorded large complex earthquakes (e.g., the 2016 Mw7.8 Kaikōura event, Klinger et al. 2018), small events (e.g.,  in the San Jacinto Fault Zone, Cheng et al. 2018), and laboratory-scale experiments (e.g., in high-velocity friction experiments, Passelègue et al., 2016).

We develop a diffuse description of fault slip to better understand complex volumetric failure patterns and the mechanics of slip in diffuse fault zones. The fault is defined via a signed distance function (s(x)), which is in turn used to define a fault indicator function with compact support H. If s(x) > H the material behaves as a pure elastic solid - otherwise the tangential stress is governed by a frictional sliding law.
Our approach is implemented on a structured hexahedral mesh using a spectral finite element (continuous Galerkin) method for wave propagation using PETSc. Our diffuse fault SEM method is inspired by the stress-glut method of Andrews, 1999.  A non-mesh aligned embedded diffusive discontinuity allows for complex dynamic rupture simulations. We present 2D numerical experiments of kinematically driven rupture and spontaneous dynamic rupture on non-planar and non-mesh aligned complex fault geometries. The method can be used to model earthquake rupture dynamics on specifically complex and evolving fault faults such as the San Jacinto, CA, fault, or shallowly dipping megathrusts and splay faulting structures in subduction zones.

How to cite: Hayek Valencia, J. N., May, D. A., and Gabriel, A.-A.: Non-planar dynamic rupture modelling across diffuse, deforming fault zones using a spectral finite element method with a non-mesh aligned embedded diffuse discontinuity, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15686,, 2021.

Abolfazl komeazi, Yi Zhang, Luca de Siena, Georg rümpker, Boris Kaus, Miriam Reiss, César Castro, Arne Spang, Philip Hering, Andreas Junge, and Tobias Baumann

What is the effect of crustal melt accumulation on the seismic wavefield? Can we reproduce the dispersion, scattering and associated stress-anisotropy with modeling tools? By performing numerical experiments of seismic wave propagation in a synthetic and geodynamically-consistent volcanic system we can test our ability to model the seismic wavefield and to reconstruct the target “magma chamber”.

We built a synthetic volcano based on recent seismic observations at the Oldoinyo Lengai volcanic complex. The velocity model  is based on a geodynamic model that provides shear modulus, Poisson's ratio, and density. The isotropic P- and S-wave velocities can be computed directly from these parameters. To test a more realistic depth dependence, we introduced a reference 1D velocity model for Northern Tanzania and expanded this to 3D. Then, we inserted variations in the rock parameters mimicking a magma chamber and resolved it using the Fast Marching Travel Time tomography code.

To further our understanding, we also added  3D anisotropy and random velocity fluctuations to the system, acting both as synthetic input for future applications and testing of seismic techniques (e.g., shear wave splitting analysis) and as noise for the travel time tomography. For the waveform modeling we used the velocity-deviatoric stress-isotropic pressure equations together with perfectly matched layers. Also, we encoded the boundary condition between solid and air in this formulation. The 25 receivers with their real geographic locations were placed for inversion sensitivity analysis. In particular, the ability to reconstruct the magma chamber and the effect of anisotropy and velocity fluctuations at frequencies up to 5 Hz are evaluated. The results are compared with a parallel forward modeling and inversion of synthetic MT data. To confirm our results and as an additional test, we also employ adjoint tomography based on spectral element method to implement a forward waveform modeling and inversion using the tools provided in the SPECFEM3D_Cartesian package.

The results present a better idea of how to construct a realistic synthetic volcano in the future. By combining multiple seismic forward models and inversion approaches, this study yields insights into the sensitivity of the seismic wavefield to geodynamically-consistent volcanic structures.

How to cite: komeazi, A., Zhang, Y., de Siena, L., rümpker, G., Kaus, B., Reiss, M., Castro, C., Spang, A., Hering, P., Junge, A., and Baumann, T.: Synthetic seismic modeling and inversion for the Oldoinyo Lengai volcanic complex., EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-16028,, 2021.

Robert Viesca

In models of faults as elastic continua with a frictional interface, earthquake nucleation is the initiation of a propagating dynamic fault rupture nucleated by a localized slip instability. A mechanism capturing both the weakening process leading to nucleation as well as fault healing between events, is a slip rate- and state-dependent friction, with so-called direct effect and evolution effects [Dieterich, JGR 1979; Ruina, JGR 1983]. While the constitutive representation of the direct effect is theoretically supported [e.g., Nakatani, JGR 2001; Rice et al., JMPS 2001], that of the evolution effect remains empirical and a number of state-evolution laws have been proposed to fit lab rock friction data [Ruina, JGR 1983; Kato and Tullis, GRL 2001; Bar-Sinai et al., GRL 2012; Nagata et al., JGR 2012]. These laws may share a common linearization about steady-state, such that a linear stability analysis of steady, uniform sliding yields a single critical wavelength for unstable growth of perturbations [Rice and Ruina, JAM 1983]. However, the laws’ differences are apparent at later, non-linear stages of instability development.

Previously, we showed that instability development under aging-law state evolution could be understood in terms of dynamical systems [Viesca, PR-E 2016, PRS-A 2016]: the non-linear acceleration of slip occurs as the attraction of a fault’s slip rate to a fixed point, corresponding to slip rate diverging with a fixed spatial distribution and rate of acceleration. Here we show that this framework can also be applied to understand slip instability development under all commonly used evolution laws, including the so-called slip and Nagata laws. To do so, we develop an intermediate state evolution law that transitions between the slip and aging laws with the adjustment of a single parameter. We show that, to within a variable transformation, the intermediate law is equivalent to the Nagata law and that fixed-point blow-up solutions exist for any value of the transition parameter. We assess these fixed-points’ stability via a linear stability analysis and provide an explanation for previously observed behavior in numerical solutions for slip rate and state evolution under various evolution laws [Ampuero and Rubin, JGR 2008; Kame et al., 2013; Bar-Sinai et al., PR-E 2013; Bhattacharya and Rubin, JGR 2014].

How to cite: Viesca, R.: State evolution laws and earthquake nucleation, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-16502,, 2021.

Breakout live text chats for every vPICO