Broadband Dynamic Rupture Modeling With Fractal Fault Roughness, Frictional Heterogeneity, Viscoelasticity and Topography: The 2016 Mw 6.2 Amatrice, Italy Earthquake

Advances in physics‐based earthquake simulations, utilizing high‐performance computing, have been exploited to better understand the generation and characteristics of the high‐frequency seismic wavefield. However, direct comparison to ground motion observations of a specific earthquake is challenging. We here propose a new approach to simulate data‐fused broadband ground motion synthetics using 3D dynamic rupture modeling of the 2016 Mw 6.2 Amatrice, Italy earthquake. We augment a smooth, best‐fitting model from Bayesian dynamic rupture source inversion of strong‐motion data (<1 Hz) with fractal fault roughness, frictional heterogeneities, viscoelastic attenuation, and topography. The required consistency to match long periods allows us to quantify the role of small‐scale dynamic source heterogeneities, such as the 3D roughness drag, from observational broadband seismic waveforms. We demonstrate that 3D data‐constrained fully dynamic rupture synthetics show good agreement with various observed ground‐motion metrics up to ∼5 Hz and are an important avenue toward non‐ergodic, physics‐based seismic hazard assessment.

The 24 August 2016, M w 6.2 Amatrice earthquake (Chiaraluce et al., 2017;Michele et al., 2020) is the first in the Amatrice-Visso-Norcia earthquake sequence in the Central-Northern Apennine system of NW-SE aligned normal faults. It was the sequence's most destructive event, causing extensive damage to surrounding buildings and infrastructure (Michele et al., 2016). The earthquake was recorded by a remarkably dense network of strong ground motion instruments, including 20 near-source stations within a radius of 50 km from the earthquake epicenter ( Figure 1, Table S1 in Supporting Information S1). The two closest stations, in Amatrice (AMT) and Norcia (NRC), are located only a few kilometers away from the fault.
The source process of the Amatrice event has been imaged using seismic data (Pizzi et al., 2017;, geodetic data (e.g., Cheloni et al., 2017;Walters et al., 2018), or both (Cirella et al., 2018;Kheirdast et al., 2021), suggesting pronounced source heterogeneities. However, kinematic finite-fault inversions are challenged by inherent non-uniqueness (Gallovič & Ampuero, 2015;Mai et al., 2016;Ragon et al., 2018;Shimizu et al., 2020;Tinti et al., 2021). Dynamic source inversions recovering friction parameters and the initial state of fault stress offer a data-driven source description compatible with earthquake physics (e.g., Peyrat & Olsen, 2004) but require a sufficiently simple dynamic rupture model to reduce the computational cost of the forward problem. A Bayesian dynamic inversion using the Parallel Tempering Monte Carlo algorithm (Gallovič et al., 2019a;Sambridge, 2013) Table S1 in Supporting Information S1). Mesh elements are colored by shear wave velocity (V s ). (b, c) Fault slip for the smooth Bayesian dynamic source inversion reference model (b) and the broadband dynamic rupture model (c). Black curves represent rupture front contours every 1 s.

10.1029/2022GL098872
3 of 13 was applied to the Amatrice earthquake, utilizing band-pass filtered (between 0.05 and 0.5-1 Hz) strong ground motion data by Gallovič et al. (2019b). Assuming a 1D medium with planar topography, the best-fitting model was used to predict ground motions up to higher frequencies than considered in the inversion (up to ∼5 Hz). Yet this approach poorly matched the high-frequency content, presumably most sensitive to unresolvable small-scale features of the rupture process.
We here propose a new approach to simulate data-fused broadband ground motion synthetics using 3D dynamic rupture modeling. Our starting point is the best-fitting model from the Bayesian dynamic source inversion of the Amatrice earthquake (Figure 1b, hereafter named "reference model"). We self-consistently augment this smooth reference model by adding fault roughness, small-scale frictional heterogeneities, viscoelastic attenuation, and topography yielding realistic high-frequency radiation without disrupting the large-scale characteristics of the reference model. The synthetic near-field ground motions show good agreement with various observed ground-motion metrics up to frequencies of ∼5 Hz.

Ingredients for Broadband Dynamic Rupture Modeling
While a wide range of mechanisms may enhance high-frequency radiation, the here selected processes have been proposed to be first-order relevant (Bruhat et al., 2020;Fang & Dunham, 2013;Pitarka et al., 2021;Shi & Day, 2013;Takemura et al., 2015;Withers et al., 2018) and are reasonably well constrained beyond the case of the Amatrice earthquake. For example, Bydlon and Dunham (2015) find that fault roughness, and not material heterogeneity, dominates the dynamic rupture process. Ma et al. (2008) discuss that large-scale topography can have locally stronger effects on modeled ground motions than 3D subsurface structure. Contrary to the latter, the topography is typically well-constrained by observations and readily available in high resolution. Viscoelastic attenuation is important to capture the decay of seismic energy with increasing propagation distance (e.g., Hu et al., 2022;Wollherr et al., 2019). Computational advances (e.g., Heinecke et al., 2014;Premus et al., 2020) now allow us to show in fully dynamic rupture models of a real earthquake that fault roughness, frictional heterogeneity, topography, and attenuation have complementary effects. Dynamic rupture source complexity enhances high-frequency generation, while topography effects elongate the synthetic coda signals, together yielding more realistic ground motions. We build our model upon Bayesian dynamic rupture inversion of the 2016 Amatrice earthquake following the approach of Gallovič et al. (2019b) with the improved forward solver FD3D_TSN (Premus et al., 2020), which was verified in a suite of dynamic rupture benchmarks (Harris et al., 2018). The inversion is performed for a 30 km long and 14 km wide planar fault governed by a slip-weakening friction law (Ida, 1972;Palmer et al., 1973). The dynamic rupture slip rate functions along the fault are convolved with pre-calculated Green's functions representing impulse responses of the medium. In this step, the fault is dipping at 45°, embedded in the 1D velocity structure of Ameri et al. (2012) with a flat free surface. The dynamic models are characterized by three spatially heterogeneous parameters: (a) the initial shear stress along dip τ i , (b) the friction drop, μ s -μ d , with μ s and μ d the static and dynamic friction coefficient respectively, and (c) the slip-weakening distance D c . Yielding occurs when the shear stress τ reaches the fault strength τ s = μ s σ n , where σ n is assumed as linearly depth-dependent normal stress with a gradient of 8.52 MPa/km and a minimum value of 0.1 MPa. The initial along-strike shear stress τ strike0 is assumed to be zero. The dynamic friction coefficient μ d is fixed to 0.4, and frictional cohesion of 0.5 MPa is assumed everywhere on the fault. The best-fitting model from the Bayesian inversion represents the reference model of this study. We then perform high-resolution enhanced 3D dynamic rupture simulations using the open-source software package SeisSol (https://github.com/SeisSol/SeisSol), resolving seismic wavefield up to 5 Hz (locally up to 10 Hz) within 50 km distance of the fault using an unstructured, statically adaptive mesh consisting of 80 million tetrahedral elements (Figure 1a, Text S1 in Supporting Information S1).

Fault Roughness
The reference model's dynamic parameters are first bilinearly interpolated from their 1.75 km along-dip and 2.3 km along-strike reference resolution into a denser 25 m sampled grid (see Figure 2, column a). Next, we adapt the fault morphology to adhere to a band-limited self-similar (Hurst exponent H = 1) fractal surface. The amplitude-to-wavelength ratio α of natural faults ranges between 10 −4 and 10 −2 (Power & Tullis, 1991), and we here use α = 10 −2 allowing direct comparability to earlier studies (Bruhat et  Shi & Day, 2013;Withers et al., 2018). This nearly upper limit of roughness may be related to the presumably immature fault system hosting the Amatrice earthquake (Pizzi et al., 2017), as suggested from regional slow long-term slip rates (Galadini and Galli, 2003), the young age of post-orogenic extension in the Apennines and the decoupling effect of multiple décollement levels and strong rheological contrasts (e.g., Ascione et al., 2013;Tondi et al., 2020).
The fractal surface wavelengths are band-limited between λ min and λ max . Choosing λ min = 200 m balances resolution requirements and computational cost for our setup and aligns with previous 3D fault roughness studies (Fang & Dunham, 2013;Shi & Day, 2013). Our choice of λ max = 2 km is motivated by the ∼2 km spatial resolution of the dynamic parameters in the reference model.

3D Roughness Drag and Heterogeneous Initial Stresses
Shear and normal stresses are dynamically perturbed by fault roughness during rupture propagation. The general scaling of the "roughness drag" (Dunham et al., 2011), an additional shear resistance to slip, was derived for a 1D rough fault in a 2D quasi-static boundary perturbation analysis by Fang and Dunham (2013) with Δ being the fault slip, λ min the minimum roughness wavelength, and G* = G/(1 − ν), where G and ν are shear modulus and Poisson's ratio, respectively, and λ max >> λ min .
To preserve the overall characteristics of the reference scenario while incorporating fault roughness, we must compensate the roughness drag in the initial loading by increasing the reference initial shear tractions as dip = dip0 + drag 3D and strike = strike0 + drag 3D . We numerically approximate the roughness drag (following Dunham et al., 2011, but for the first time based on 3D dynamic rupture models) as where C is a dimensionless coefficient. In Text S2 in Supporting Information S1, we demonstrate in numerical experiments that C can be approximated from τ drag 2D using characteristics of the reference model slip distribution. For our choice of n = 4 elements to resolve λ min = 200 m, we obtain C of ∼0.44. The average value of τ drag 3D across the rupture area is ∼1.4 MPa. Our lower τ drag 3D is intuitive since the material off a 3D roughness feature can also deform perpendicular to slip in distinction to the 2D case. We caution that while τ drag 2D is defined for 2D in-plane quasi-static rupture, proper analytical treatment of the 3D (dynamic) roughness drag is very complex. The formulation of the roughness drag is based on stationary statistics and thus applies in the slipping (or slipped) region well behind the rupture front, challenging analytical extensions to account for rupture front curvature, inertia, or slip gradients in the vicinity of the rupture front. This motivates our empirical approach (see also Text S2 in Supporting Information S1). We account for the 3D roughness drag while preserving the smooth reference initial stress distribution by loading the rough fault with a heterogeneous regional stress tensor: we first adapt the smooth reference initial fault loading to balance roughness drag, then expose the now rough fault to the adapted loading (Text S3 in Supporting Information S1). As a result, the broadband model features roughness-induced small-scale fluctuations of the initial shear and normal tractions (Figure 2b), consisting of both releasing and restraining slopes that bring the fault closer and farther from failure, respectively ( Figure S2 in Supporting Information S1).

Frictional Heterogeneity
We perturb the smooth variation of the reference characteristic frictional slip weakening distance D c 0 , the spatially most variable dynamic parameter in the reference Bayesian dynamic inversion. The relative standard deviation of D c is on the order of 50% (Gallovič et al., 2019b), highlighting its importance as a proxy for unaccounted geometrical and geological features. We use a band-limited fractal distribution. We prescribe D c = max (0.14 m, D c 0 (1 + ε)), where 0.14 m is the minimum value of D c 0 , and ε follows a fractal distribution of amplitude-to-wavelength ratio α = 10 −2 generated from a different random seed than the one used for the fault roughness. Including small-scale heterogeneous D c in our broadband model is a proxy for frictional or stress asperities that have been shown to be important for high-frequency radiation in previous work (e.g., Galvez et al., 2020;Ripperger et al., 2008). Heterogeneous D c contributes to radiating additional energy due to fault-local acceleration and deceleration.

Topography and Viscoelasticity
In our broadband dynamic rupture, the flat free surface used in the inversion is superseded by high-resolution topography data sampled to 150 m resolution (Farr et al., 2007). The modeled 3D domain spans 300 × 300 km horizontally and extends to a depth of 150 km to avoid any undesired reflections from the (imperfectly) absorbing boundaries. We incorporate the 1D velocity model, with V p = 1.86V s , and viscoelastic attenuation, with Q p = 2Q s , inferred by Ameri et al. (2012), see Table S2 in Supporting Information S1. Accounting for topography and viscoelasticity is complementary to including dynamic source heterogeneity. Realistic 3D topography scattering redistributes seismic energy to later arrival times enhancing synthetic seismogram coda signals, an effect that cannot be obtained when considering source complexity only.

Rupture Dynamics
We compare the broadband dynamic rupture model, incorporating fault roughness, small-scale D c variation, and topography to the reference model with a planar fault, a flat free surface, and smoothly varying initial conditions, in terms of fault slip (Figures 1b and 1c), slip rate space-time evolution (Figures 2c and 2d), and moment rate release ( Figure S3 in Supporting Information S1). Figures 2c and 2d demonstrate similar large-scale slip evolution. The seismic moment of the broadband model is 2.8 × 10 18 Nm, corresponding to M w = 6.24, which is comparable to the reference model with 2.6 × 10 18 Nm seismic moment (M w = 6.20). We highlight that both models recover the remarkably weak and slow nucleation phase (Gallovič et al., 2019b;, as was also inferred for the Norcia earthquake . The nucleation is followed by bilateral rupture, which is slower toward the NW than toward the SE in both models. At smaller scales, the broadband model features decoherence of rupture fronts (Shi & Day, 2013). Locally fluctuating rupture speeds are due to acceleration and deceleration at releasing and restraining slopes, heterogeneous initial shear and normal traction, and D c heterogeneity. Peak slip rates are increased by ∼15% in the broadband model, while both models feature pulse-like ruptures, and rise time remains largely unaffected.

Figures 1b and 1c and
Comparisons of moment rate releases ( Figure S3a in Supporting Information S1), moment rate spectra ( Figure  S3b in Supporting Information S1), and the second time-derivative of moment rate releases ( Figure S3c in Supporting Information S1) illustrate the effects of the fault roughness, heterogeneous loading, and D c on the high-frequency rupture radiation. While the two distinct episodes of moment rate release are well recovered, its first peak is about 20% higher in the broadband model than the reference model ( Figure S3a in Supporting Information S1), reflecting the increase in negative strength excess in the nucleation region required for broadband self-sustained spontaneous dynamic rupture across the rough, frictionally heterogeneous fault. Figure 3 compares the observed three-component velocity and acceleration waveforms recorded at the 20 strong-motion stations (Figure 1, Luzi et al., 2016) with synthetics from the broadband dynamic rupture model. The overall agreement in terms of waveform shape and duration is good for both velocity and acceleration waveforms. The synthetic amplitudes fit velocity recordings very well at most stations ( Figure 3a). Nevertheless, the modeled acceleration amplitudes are significantly underestimating some station components (Figure 3b, e.g., NRC, MNF, FOS, ASP). Analogous plots for the reference model and the rough, heterogeneous fault model without topography are shown in Figures S4 and S5 in Supporting Information S1, respectively. We isolate the effect of using small-scale heterogeneous D c in comparison to the smoothly varying D c from Gallovič et al. (2019b) in Figures S3 (moment rates) and Figure S6 (synthetics) in Supporting Information S1.

Ground Motions
To highlight the role of the fault roughness, frictional heterogeneity, and topography, Figure 4 compares observed EW velocity and acceleration waveforms and Fourier Amplitude Spectra (FAS) with synthetics of three models: the reference, the broadband rough fault model with topography, and the broadband rough fault model without topography. All components and stations are shown in Figures S7-S13 in Supporting Information S1. The synthetic waveforms of the broadband models match long-period data (0.05-0.5 Hz) equally well as the reference model. Nevertheless, the reference model provides waveforms clearly depleted at high frequencies. A general trend is that fault roughness and topography enhance and elongate waveforms at high frequencies, respectively, although not to the same extent at all stations. The increase in high-frequency content in the broadband waveforms  (Figure 1), ordered by epicentral distance. Synthetics are from the broadband dynamic rupture scenario incorporating fault roughness, D c heterogeneity, and topography. Both observed and synthetic waveforms are scaled by their maximum value, which is indicated on the right-hand side of each plot. Velocity waveforms are scaled by the maximum value of the observed records at each station, while acceleration waveforms are scaled component-wise.
without topography (green) is clearly limited in duration compared to the same model incorporating topography (red). High-frequency ground motions are amplified early-on by fault roughness, while topography-induced scattered waves prolong their duration. The combination of both effects is most pronounced in the central and SE part of the hanging wall region (see, e.g., stations PZI1, LSS, and SPD in Figure 4 and Figure S8 in Supporting Information S1, or MSC, AMT, and ANT in Figures S7, S10, S11, and S13 in Supporting Information S1). At some stations (e.g., stations FOS and ASP), our broadband synthetic spectra are improved yet underestimate the observed spectra at frequencies higher than 1 Hz. We may speculate that larger station-source distances, as for FOS and ASP, and local subsurface complexity such as basin edge effects render more realistic 3D velocity models than used here, important for physics-based broadband modeling. Similarly, fully considering site effects may amplify high-frequency waveform spectra.
Animations of the three components of the velocity wavefield for the reference and broadband models are shown in Movies S2-S7 in Supporting Information S1. They illustrate how seismic waves are both reflected and scattered upon propagating across sharp topographical features like mountains and hills, which explains the prolonged duration of the seismic signal for several receivers (e.g., stations LSS and SPD, Figure 4). Although viscoelastic attenuation is generally important to capture the decay of seismic energy during topography scattering, a comparison between our models with and without attenuation in Figure S14 in Supporting Information S1 shows relatively minor effects (e.g., for stations SPD, PZI1, TERO).
Figures S15 and S16 in Supporting Information S1 quantify the fit of the synthetic ground motions of the broadband and reference models with observations using Goodness-of-fit (GOF) metrics (Olsen & Mayhew, 2010), including peak ground velocity and displacement, spectral acceleration, FAS, energy duration, and cumulative energy (Text S4 in Supporting Information S1). The broadband model with topography fits the observations  Konno and Ohmachi (1998). The observed data are tapered with a 35 s cosine window. 9 of 13 better (GOF 45-65) than the reference model  and the broadband model without topography (GOF 40-60). Figure S17 in Supporting Information S1 details the model bias and standard deviation over the 0.5-10 s period range, averaged over 20 stations used in this study. In general, a near-zero model bias over a specific period suggests that simulated ground motions match observations reasonably well at given frequencies. The reference model fits the observations only at periods longer than 2 s. Compared to the reference model, the fit of the broadband model without topography ( Figure S17b in Supporting Information S1) is improved (30%-40% lower bias) at periods shorter than 2 s. The broadband model with topography shows an even better fit (40%-50% lower bias than the reference model, Figure S17c in Supporting Information S1). However, while both broadband models preserve a perfect fit at periods longer than 2 s, some synthetics still underestimate the observations at periods shorter than 2 s, as discussed above."

Discussion
Recorded broadband ground motions are widely used in earthquake engineering to inform the performance-based design of structures. Typically, generic strong-motion waveforms that fit specific ground motion metrics are selected from a strong-motion database for that purpose (Iervolino et al., 2010). Also, probabilistic seismic hazard analyses often rely on such so-called ergodic ground-motion models (GMMs, e.g., Petersen et al., 2019). Yet, these may not reflect the conditions of a specific region of interest. Regional synthetic broadband ground motions from 3D dynamic rupture inversions, which offer a physically consistent representation of earthquakes, can sample conditions that are not sufficiently constrained in empirical models toward the development of non-ergodic, physics-based GMMs Graves et al., 2011;Moschetti et al., 2017;Wirth et al., 2018;Withers et al., 2020). Our proposed broadband dynamic rupture models can be extended to account for other distinctive regional characteristics, such as a listric or segmented fault geometry, 3D velocity models including low velocity layers and basins, and fault zone plasticity (Roten et al., 2014). They may also inform PSHA-targeted kinematic rupture generators while inherently ensuring realistic scaling of earthquake characteristics (e.g., Savran & Olsen, 2020). Our models emphasize the need to include (a) small-scale source characteristics to enhance the high-frequency source radiation during the rupture propagation, and (b) topography to increase the duration due to scattering. The duration of the latter effect is controlled by viscoelastic attenuation.
We carefully analyze the effects of adding roughness to a flat fault model. We counterbalance the consequent 3D roughness drag by increasing initial shear traction by τ drag 3D (Equation 2), calculated using the spatially variable slip amplitude Δ of the reference model. We explored an alternative model (not presented here), with constant Δ equal to the peak slip of the reference model (1.14 m). It generates a higher average τ drag 3D of about 3.3 MPa (cf. 1.4 MPa, Section 2.2). It may be possible to identify alternative satisfying models based on constant Δ. Nevertheless, the here presented approach of constraining Δ by the spatially variable reference fault slip appears superior due to its simpler and better constrained parametrization.
Our proposed approach to enhance smooth slip models for broadband dynamic-rupture simulations is independent of the data and generation process of the starting smooth source model. Besides the dynamic source inversion used here, kinematic slip models can serve as a reference model as well (e.g., Ma et al., 2008). However, the choices required for conversion in the latter approach may lead to several families of plausible dynamic rupture scenarios that need additional (e.g., geological) constraints to be distinguished .
Although our rough fault model with topography improves the waveform fit at high frequencies, some synthetics still underestimate the observations. Additionally, a good fit of shape, duration and amplitude characteristics of broadband velocity waveforms up to 5 Hz, does not necessarily translate in capturing broadband acceleration amplitudes equally well, a challenge identified in kinematic source modeling as well (e.g., Irikura & Miyake, 2021). More complete matching of the observed records at periods <2 s may in future be enabled by: (a) considering smaller length scales (λ min ) of fault surface roughness, potentially further increasing high-frequency radiation at the cost of increased computational demands; (b) incorporating larger-scale non-planar fault geometry such as a listric fault geometry which has been, for instance, suggested from satellite data (Tung & Masterlark, 2018) and which may modulate peak ground velocities as a consequence of curvature focusing effect (Passone & Mai, 2017); (c) probing and quantifying the variability of the predicted shaking using alternative models from the Bayesian ensemble of the dynamic rupture inversion, (d) incorporating a more realistic Earth model to capture path and local site conditions, that is, 3D velocity models, small-scale scattering media (Bydlon & Dunham, 2015;Imperatori & Mai, 2013), site corrections (Rodgers et al., 2020) or non-linear soil effects (Roten et al., 2012) and (e) off-fault damage (e.g., Okubo et al., 2019;Yamashita, 2000). In particular, low-velocity sedimentary basins (Lee et al., 2009;Pischiutta et al., 2021) can significantly amplify the amplitude and duration of ground motions, which may lead to improved synthetics for stations with strong site-effects, for example, CLF with site-class D (Table S1 in Supporting Information S1). However, accounting for these mechanisms in sufficient detail (at scales down to ∼100 m or less) in observationally constrained 3D broadband dynamic rupture models pose additional computational and observational challenges and their relevance for matching real earthquake recordings remains open to debate.

Conclusions
We present a novel approach for broadband dynamic rupture modeling constrained from low-frequency data toward generating physics-based, non-ergodic ground motion synthetics validated by observations. We generate broadband dynamic rupture models of the 2016 M w 6.2 Amatrice earthquake by combining large-scale heterogeneous stress and frictional parameters, inferred from the best-fitting model of a Bayesian dynamic rupture inversion, with small-scale self-similar fault roughness and frictional (slip weakening distance) heterogeneity, topography, and viscoelastic seismic attenuation. We empirically quantify the 3D roughness drag governing rupture dynamics on small scales by counterbalancing its effective dynamic stress perturbations. We obtain dynamic rupture scenarios that successfully reproduce the low-frequency (<1 Hz) source characteristics of the inverted dynamic model. The combined small-scale heterogeneities of the fault geometry, frictional strength and loading enhance high-frequency source radiation. Topography elongates the synthetic waveforms by enhancing coda effects. In combination, we obtain more realistic high-frequency (up to ∼5 Hz) synthetics comparable with observed strong motion records. Our work demonstrates 3D physics-based, broadband earthquake ground-motion simulations that are tightly constrained by data-driven dynamic earthquake source inversion and allows us to quantify the first-order role of large-and small-scale dynamic source heterogeneities in the broadband seismic wavefield. Challenging future developments may focus on including 3D and site-specific wave propagation effects toward realistic fully physics-based acceleration synthetics suitable for engineering applications.