Physics-based broadband characterization of weak earthquakes

Domain: earth_science
Published: 2025-12-08
Editorial outcome: major_revision

Abstract
Analyses of abundant small earthquakes have the potential to map the Earth’s small-scale stress and frictional properties. However, standard seismological characterizations lack the resolution to capture the physical complexity of earthquake rupture, particularly its spatial and temporal heterogeneity. Advanced dynamic rupture models, which integrate elastodynamic simulations with frictional laws, have so far been applied to large earthquakes. Here, we develop a Bayesian dynamic inversion from station-specific (apparent) source spectra, extending up to 25 Hz, for slip-weakening friction, assuming heterogeneous parameterization along a finite-extent planar fault with resolution down to ~100-m scale. The approach is demonstrated on two Mw ~ 4 earthquakes in Central Italy with distinct spectral behavior, one directive and one nondirective. Results show that the inversion resolves key mean source parameters, pinpoints possible pitfalls in standard estimates, and infers a power-law decay of stress and friction heterogeneity spectra down to the smallest scales. Such advanced studies promise to unravel so-far elusive small-scale characteristics of earthquake ruptures.

Introduction
Earthquakes are presumably self-similar phenomena, as suggested by observations indicating that stress drop remains independent of earthquake magnitude across a broad range of event sizes1,2. While large earthquakes are routinely studied by detailed kinematic and dynamic slip inversions thanks to seismological and geodetic observations at short distances, analyses of smaller events typically rely on fitting spectral source models to processed P- or S-wave traces. Rare kinematic slip inversions of small earthquakes generally demonstrate spatial variability in slip distribution like that well known for large earthquakes3–5. A similar conclusion was reached by Yoshida6 from the observed complexity of the shapes of the direct P waves of small events observed in near-source borehole data in Japan.

The temporal evolution of earthquake source radiation as observed by various receivers is described by an apparent (station-dependent) source time function (ASTF)1. ASTF or their amplitude spectra (so-called apparent spectra, AS) are commonly used to estimate earthquake source characteristics (rupture size, stress drop, etc.) under the assumption of simplified7,8 source models. These models generally assume an omega-squared source spectrum with a corner frequency depending on the seismic moment and stress drop. However, this is an oversimplification because the corner frequency (as well as the whole spectrum) generally depends on the shape of the ASTF9–12. Empirical source spectra often exhibit high-frequency decay that deviates from the classical omega-squared model. This behavior is sometimes characterized using the so-called κ-source parameter, which quantifies attenuation at high frequencies13.

Due to the rupture directivity effect, the ASTF and AS can exhibit spatially strongly varying apparent duration and corner frequencies, respectively14–18. Azimuthal variations in corner frequency and spectral decay from smooth dynamic ruptures have been extensively explored by Kaneko and Shearer19, who suggested that properly accounting for the spatial dependence of source spectral properties could provide valuable constraints on source parameters. Yoshida20 demonstrated that asymmetrical rupture propagation may result in biased stress drop estimates if the directivity effect on the source spectra is neglected. Pacor et al.16 and Colavitti et al.14 demonstrated that the empirical AS as a whole (not only the corner frequency) can exhibit complex azimuthal dependence even for weak directive events in Central Italy, assigning it tentatively to rupture process heterogeneity. Schliwa and Gabriel21 demonstrated by numerical simulations that the local rupture heterogeneities are imprinted in the spatial distributions of apparent corner frequencies. Calderoni and Abercrombie22 point out the trade-off between high stress drop and increased rupture complexity, as both factors contribute to stronger high-frequency ground motions. Nevertheless, the stress drop and rupture complexity likely affect the radiated wavefield with distinct azimuthal and frequency dependence. Therefore, an inversion utilizing a suitable finite-extent physics-based rupture model might resolve and quantify such competing source effects.

In dynamic source inversions, frictional parameters on a fault are inferred from observed waveforms, obtaining a rupture model constrained by both the data and the governing physics (elastodynamic equation and friction law). Bayesian dynamic rupture inversions of regional waveforms have already been performed for several large events23–26. The waveform inversions are generally limited by the accuracy of Green’s functions describing the path and site effects, restricting their application to low-frequency waveforms. As a result, they are unsuitable for analyzing small earthquakes and unraveling small-scale source properties. Therefore, an inversion of AS corrected for path and site effects in a broad frequency range represents a tempting alternative to constrain a dynamic model over a wide range of scales. Yet, up to now, there has been no attempt to invert the AS for dynamic rupture parameters directly.

Although earthquake ruptures are generally affected by many physical source phenomena (heat production, off-fault fracturing and plastic deformation, fault non-planarity, inhomogeneous fault rheology, etc.), here we attempt to use a relatively simple dynamic model yet with generally heterogeneous prestress and slip-weakening parameters along a planar fault embedded in a homogeneous space. Whereas theoretical studies often employ smooth models with strong barriers to reproduce high-frequency spectra, we instead adopt a heterogeneous distribution of dynamic parameters to capture rupture complexity. The source heterogeneity is suggested across multiple scales in, e.g., laboratory experiments27–29, field observations of natural faults30–32 or strong motion simulations33–38. Note that, unlike kinematic source models, the rupture dynamics guarantee that the inferred models are physically consistent, thus serving as a physical constraint in the inverse modeling.

As input data for the inversion, we use apparent (station-dependent) source spectra (AS) up to 25 Hz obtained by the Generalized Inversion Technique (GIT39–42). GIT is a non-parametric approach that decomposes recorded spectra into source, path, and site contributions, thereby providing robust corrections for path and site effects and yielding reliable source spectra for rupture characterization. Since the spectra do not include the phase information, we account for the associated uncertainties within the Bayesian framework. To our knowledge, a dynamic inversion of apparent amplitude spectra has not been attempted previously; therefore, we perform a synthetic test to guide the interpretation of the real-data inversion.

In this work, we focus on small M4 earthquakes in Central Italy to exploit the full frequency range while making the dynamic simulation computationally tractable. The Central Italy region, one of the most seismically active in Italy, has produced tens of thousands of well-recorded events (by up to hundreds of stations) and thus has been subject to numerous source studies addressing stress drop scaling, though all assume simplified source spectral model22,43–45. The present study develops a methodology that relaxes the assumption of simplified source radiation while employing similar input data. We select two events yet with distinct AS behavior to be compared in terms of their source characteristics and variability within the Bayesian framework. Assuming a heterogeneous distribution of model parameters, we recover complex rupture propagation down to the smallest scales ~100 m that is both physics and data constrained. In addition, we compare our results with those determined by common seismological practice based on standard source spectral models and finally, investigate the power spectral densities of the heterogeneities in the inferred dynamic rupture parameters.

Results
The two example events occurred during 2016-2017 Central Italy sequence: an Mw 4.2 aftershock of the 2016-10-30 Norcia mainshock and an Mw 4.5 aftershock of the 2016-08-24 Amatrice earthquake. Figure 1b shows the AS for the two events as a function of frequency with color-coded station azimuth to emphasize the potential systematic dependence and spatial distribution of the AS amplitudes. It also highlights different AS variability between the two events over the whole frequency range. Figure 1a then shows the deviation of the AS from the average at an example frequency of 10 Hz in a map view, confirming no systematic dependence on the distance from the source. The azimuthal coverage is good with 16° and 35° gaps for the directive and nondirective events, respectively, yet with a prevalence of stations in the NW-SE direction (see the rose diagram in Fig. 1a), following the orientation of the Apennines. Visual inspection of Fig. 1 suggests that the Mw 4.2 earthquake exhibits a strong directivity effect towards the south, while no spatial trend is visible for the Mw 4.5 event.

Best rupture models
We examine the maximum a-posteriori (MAP) models of the two events (i.e., having the lowest misfit). Figure 2a displays their kinematic characteristics (slip, stress drop, rupture time, and rupture speed), and Supplementary Movies S1, S2 show their rupture propagation. Models of both earthquakes are complex down to the smallest scales, featuring locally high stress drops and numerous local decelerations and accelerations of the rupture up to super-shear speeds. As expected, the directive Mw 4.2 event has a more asymmetric rupture due to the nucleation point being closer to one side of the ruptured area, unlike the nondirective Mw 4.5 event. The Mw 4.5 earthquake ruptured a larger portion of the fault and had larger stress-drop peaks than the Mw 4.2 event. Figure 2b displays dynamic parameters confirming the significant degree of spatial heterogeneity of the dynamic parameters, which translates to the heterogeneity of the slip distributions and rupture speeds described above.

Figure 3a displays the bias between the observed and synthetic AS, evaluated as the natural logarithm of their ratio. The average bias being close to zero in the entire frequency range demonstrates that the AS’s overall level and azimuthal dependence are fitted well for both the directive and nondirective events. This is particularly interesting considering the distinct azimuthal and frequency dependence of the GIT AS of the two events (Fig. 1), see also the synthetic spectra in Supplementary Fig. S1. The standard deviation of the bias is approximately frequency-independent and similar for both events (0.45-0.46). The bias can be further explored in a map view plotted for selected frequencies in Fig. 3b, demonstrating very weak azimuthal and distance dependence. Note that the ~30% underestimation at the lowest frequencies is related to the seismic moment M0; models with higher M0 overestimate the data by a similar amount in that frequency range.

Ensemble rupture properties
Since the dynamic inversion of AS amplitude spectra has not been performed previously, we must properly inspect the uncertainties revealed by the Bayesian framework, here represented by the rupture model ensembles. In particular, we obtained 4745 and 3667 dynamic rupture models for the Mw 4.2 and Mw 4.5 events, respectively.

Figure 4a shows moment rate (MR) functions calculated as integrals of the slip rates along the fault multiplied by the shear modulus calculated from the medium properties (Table 1). The MRs are color-coded by the posterior PDF values. The MAP model is shown in black and has ~1 s duration for both events with a weak onset due to the prescribed weak nucleation (see “Methods”). The temporal spread of the MRs is relatively wide, including cases with delayed nucleation, mainly because the amplitude spectra lack explicit temporal information.

In addition, we compare the ensemble amplitude spectra of the MR second derivative (so-called moment jerk) with the spectrum averaged over empirical AS and its uncertainty (cf. green curve in Fig. 1b). In empirical studies, the station-averaged AS is considered to represent the MR spectrum to derive source characteristics. Figure 4 shows that while this assumption is reasonable for the nondirective circular Mw 4.5 rupture, where the observed average spectrum matches the ensemble MR spectra within uncertainty, it may be disputable for the directive Mw 4.2 event, where the average AS overestimates the MR spectrum. This affects the derived source characteristics of the earthquakes, as discussed later.

To reveal stable slip features, Fig. 4b displays the average slip distribution (in color scale) for the two events. Both slip models attain relatively simple, smooth, quasi-circular shapes. This suggests that the localized peaks in slip seen in the MAP models in Fig. 2a vary in position over the ensemble models, and the ensemble averaging smooths them out. Figure 4b also shows the variability of the rupture extent over the ensemble in terms of mean and one-sigma standard deviation contours of slip. We note that the individual models exhibit diverse shapes, which produce circular distribution in the ensemble averaging.

Figure 5 shows the ensemble averages of the dynamic rupture properties along the fault and their standard deviations. Both the ensemble-averaged prestress and friction drop are smoothed out, being approximately constant within the average rupture area, with increased variability towards its border due to the variability of the individual rupture extents. The constant values are similar for both events, the prestress and friction drop being in the range of 3-4 MPa and 0.08-0.09 (i.e., the strength of 8-9 MPa), respectively. The slip-weakening distance of both events increases with distance from the nucleation, ρ, approximately as Dc(ρ) = 0.35 [cm] + 3.8ρ [km] and Dc(ρ) = 0.5 [cm] + 4.1ρ [km] for the Mw 4.2 and Mw 4.5 event, respectively. We would like to emphasize that the smooth averaged dynamic model (if able to rupture) would be unable to fit the observed AS data (see Materials and Methods and cf. Supplementary Fig. S2).

Figure 6 shows histograms of the misfit in the ensembles for the two events (Mw 4.5 in red and Mw 4.2 in blue). The misfit attains slightly lower values for the Mw 4.5 rupture, which can be attributed to the lower number of stations (93) than for the Mw 4.2 event (117). The similar misfit values confirm a similar fit of the data despite their very different azimuthal dependence and variability (Fig. 1b). Figure 6 also displays the distribution of seismic moments M0 and the assumed prior distributions. The posterior distributions are narrower than the prior, with slightly larger mean values.

For each dynamic model in the ensembles, we evaluate gross physical quantities describing the rupture process, such as the mean slip, equivalent rupture radius R = sqrt(S/π) (S being the actual ruptured area), duration, average rupture speed, mean Dc, and average stress drop, radiated energy Er, average fracture energy density (or fracture energy G over unit area), radiation efficiency defined as η = Er / (Er + G), and radiated-energy-to-moment ratio Er/M0. Figure 6 displays their marginal posterior distributions as histograms. All displayed parameters (perhaps except for mean Dc) are characterized by unimodal Gaussian-like probability density distributions. Average rupture duration, velocity, and stress drop are similar between the events despite the significant difference in their AS shapes (Fig. 1b). While the similarity of the rupture duration can be easily explained by the rupture geometry (i.e., unilateral for Mw 4.2 and circular for Mw 4.5 event), the similarity in the average velocity and stress drop is clarified by the synthetic test described in Supplementary Text S1. It shows that the average rupture velocity, including its variability, is primarily governed by the dynamic rupture with heterogeneous parameters, and therefore remains unconstrained by the AS inversion. However, the average stress drop is well resolved according to our synthetic test, and thus, the agreement of the stress drop marginal PDFs between the two events can be considered an important result of the AS inversion, suggesting a constant stress drop scaling of these two events in the same tectonic setting of Central Italy.

Figure 6 further shows how the events clearly differ in rupture radius, mean slip, Er, average fracture energy density, and mean Dc. This agrees with the constant stress-drop source scaling laws: the seismic moment differs approximately 4 times, and thus, the mean slip and rupture radius of the smaller event is consistently sqrt(1/4) = 1/2 times smaller. Similarly, the radiated energy is roughly 4 times larger for the Mw 4.5 event. Our results also agree with the empirical and theoretical findings that Dc and fracture energy rate scale linearly with slip46,47, i.e., 1.6 times in our case. Although these parameters follow the scaling laws, the synthetic test (Supplementary Text S1) suggests that the AS constrain their values significantly compared to the assumed prior.

From the source scaling, it also follows that the scaled (non-dimensional) parameters, namely radiation efficiency and radiated-energy-to-moment ratio, seem magnitude-independent (Fig. 6). The inferred radiation efficiency of 0.15-0.25 suggests that approximately 75–85% of the available potential energy was dissipated in the rupture process. The radiation efficiency is relatively small compared to estimates for other earthquakes using various methods2. Nevertheless, similarly low values were also found for large crustal events using the same inversion technique applied to the waveforms and for various friction laws23,24,26. The inferred radiated-energy-to-moment ratio of 0.7-2×10−5 is consistent with the range of estimates for other earthquakes2.

Beyond standard estimates based on corner frequency
In standard seismological practice, the source spectra are used to infer the corner frequency and estimate the stress drop from it. The standard omega-squared model of the source spectrum Ω(f) with corner frequency fc reads
Ω(f) = M0 / (1 + (f/fc)^2)

We fit Eq. (1) to the synthetic MR spectra from our ensemble using L2 minimization in logarithms, treating scalar moment M0 as known. Figure 7 shows the distributions of the estimated corner frequencies as histograms. Interestingly, the directive smaller event has generally larger corner frequencies (0.5–1.5 Hz) than the nondirective one (0.4–1.1 Hz) despite having similar durations (Fig. 6). As demonstrated by Liu et al.10 and Gallovič and Valentová9, the corner frequency is affected by the shape of the MR function, which is more complex in case of the directive event (Fig. 4). As shown in Supplementary Fig. S3, the fit is slightly worse for the nondirective Mw 4.5 event but only above ~3 Hz (i.e., well above the corner frequency) because its spectral decay is faster than the omega-squared. Note that this corner frequency estimate represents an indirect approach because it relies on fitting the MR spectra that can be obtained only by a finite-fault inversion.

We also fit the empirical spectrum of Eq. (1) directly to the apparent source spectra from GIT shown in Fig. 1b to obtain a standard seismological estimate. We calculate the L2 norm over all stations, effectively averaging over all azimuths and distances, and depict them by arrows in Fig. 7. While the resulting corner frequency of the nondirective Mw 4.5 event falls within the values inferred by the dynamic inversion, the direct estimate from the GIT spectra is significantly overestimated for the directive Mw 4.2 event. This is understandable when looking at the great (mainly azimuthal) variability of the GIT spectra shapes, here interpreted as due to the directivity effect.

The corner frequency is generally related to spectral stress drop Δσ as
Δσ = (7/16) * (M0 * fc^3) / (k * β^3)

where k is a model-dependent parameter. Typically used are the Brune7 model with k = 0.37 and the Madariaga8 model with k = 0.21. Histograms in Fig. 7 show the Brune and Madariaga stress drops calculated using Eq. (2) from the indirectly estimated corner frequencies (i.e., from the ensemble MR spectra). Note that both events are characterized by similar stress drops, as was also observed for the stress drop calculated from the inferred dynamic models (Fig. 6). Nevertheless, the Madariaga stress drops are closer to the true stress drops derived from the dynamic models (Fig. 6) compared to those estimated by the Brune model. This supports the suitability of the Madariaga model, as it is based on dynamic simulations and aligns well with the results of Kaneko and Shearer19 for a symmetric circular dynamic rupture using stacked apparent spectra (over station azimuths) and fitted by the omega-squared function.

We also calculate the Madariaga stress drop using the corner frequencies inferred directly by fitting the GIT AS to the omega-square model (shown as arrows in Fig. 7). The stress drop is overestimated significantly for the directive Mw 4.2 event, resulting from the overestimated corner frequency. We attribute this to the directivity of the rupture model and the insufficient averaging over the focal sphere19. This suggests that the simple spectral stress drop estimates from averaged AS (provided physics-based models, such as Madariaga, are utilized) may be biased in practice, especially for events with strong azimuthal variability of the AS. Therefore, the practitioners should first examine the variability of the AS shapes as their large spread may indicate a directive event for which the standardly inferred source properties might be significantly biased (Fig. 4a). A potential improvement may be in such cases represented by fitting apparent (station-dependent) corner frequencies estimated from the AS using asymmetrical rupture models20. However, the most reliable approach, regardless of the event’s directivity, is to derive directly the moment rate functions by, e.g., dynamic or kinematic finite-fault inversion48. Improved stress drops can also be obtained from rupture area estimation by the second moment approach49–51 to avoid using corner frequency and the underlying assumptions on the temporal rupture characteristics9.

Spectral properties of the dynamic rupture spatial parameters
To characterize the spatial properties of the dynamic rupture parameters, we explore the power spectral density (PSD) of the perturbations with respect to the average distribution. We recall that the average distributions are relatively smooth (Fig. 5), in contrast to the heterogeneity of the individual models (e.g., Fig. 2).

Let us select one of the dynamic parameters, g(x,y), where x and y are the fault coordinates. We define the perturbations as d(x,y) = g(x,y) - E(g(x,y)), where E represents the expectation evaluated by averaging the dynamic parameter over the model ensemble. The PSD is defined as
P(kx, ky) = (1/S) * |E(F(d(x,y)s(x,y)))|^2

where F represents the Fourier transform to wavenumbers kx and ky, and S is the fault area. Spatial taper s(x,y) is the slip distribution Δu(x,y) normalized by the square root of its average power, s(x,y) = Δu(x,y) / sqrt(integral(Δu(x,y)^2 dS / S)), to analyze only the ruptured part of the fault. We perform the analysis for g being the initial stress and strength excess (strength minus the initial stress). Since the slip-weakening distance is presumably a scale-dependent parameter that varies over an order even in our models, we consider g = log Dc. This corresponds to the analysis of the relative perturbations with respect to the median of Dc.

Figure 8 shows the resulting PSD estimates as a function of radial wavelength (λ = k^-1 = (kx^2 + ky^2)^-1/2). The PSD exhibits k^-2 spectral decay, corresponding to Hurst exponent H=0, between approximately 200-600 m long wavelengths. The maximum 600 m long wavelength is close to the estimated radius of the ruptures. Thus, the correlation length of the heterogeneities is longer than or similar to the rupture size. The minimum 200 m long wavelength can be associated with the bilinear interpolation of the parameters used in the inversion. Indeed, the spatial discretization is 125 m, and thus, the faster power-law decay at shorter wavelengths can be attributed to the smoothing effect of the spatial interpolation.

The inferred k^-2 power-law decay of PSD corresponds to k^-1 decay of the stress field, which is consistent with the results of dynamic rupture simulations. They followed the theoretical derivation by Andrews52 based on the assumption of geometrical self-similarity of the ruptures. Smith and Heaton53 estimated the falloff exponent of 1.3 for stress from the focal mechanisms of earthquakes in southern California.

The standard deviation of the perturbations (which, by definition, corresponds to the square root of the integral of PSD over wavenumbers) is specified in the respective panels of Fig. 8. The perturbations of the initial stress and strength excess are similar, 3-4 MPa. The inferred variability of log Dc (0.46–0.53) corresponds to the relative perturbation by 50–70% around the median value.

We further evaluate the coherence between the parameter perturbations d1 and d2 defined as
γ = E(integral(d1(x,y)d2(x,y)s(x,y)dS)) / sqrt(E(integral(d1(x,y)^2 s(x,y)dS)) * E(integral(d2(x,y)^2 s(x,y)dS)))

The coherence between each two parameters from the analyzed group (prestress, strength excess, and log Dc) is <0.1. This suggests that the perturbations of the three dynamic rupture parameters are generally statistically independent. We also point out that the overall incoherence in the variations of the dynamic parameters does not imply complete independence because not all random combinations of the parameter perturbations would yield a self-sustaining rupture. Indeed, the coherence analysis is limited to identifying simple, approximately linear-like relationships and may fail to capture more complex, nonlinear dependencies governed by the dynamic model. Moreover, we cannot rule out possible coherence of source features at large scales, such as the correlation of the high stress drop and high fracture energy inferred from detailed dynamic rupture inversions of, e.g., the 2020 Mw 6.8 Elazığ earthquake by Gallovič et al.24.

Let us point out that we use an idealized planar fault in the present application. In reality, rupture processes are affected by a range of physical phenomena not captured in this model, such as heat production, off-fault fracturing and plastic deformation, fault non-planarity, heterogeneous fault rheology, etc. As a result, the dynamic parameters inferred here should be regarded as effective representations of more complex underlying processes. As such, the inferred heterogeneities introduce sufficient complexity in rupture propagation to reproduce the observed AS. Direct association of our results with real physical fault zone complexities is beyond the scope of the present paper, but it can be the aim of further research.

Methods
GIT derived apparent spectra
We utilize AS obtained by the so-called non-parametric Generalized Inversion Technique, GIT39–41,54, a method used to isolate frequency-dependent source, path, and site functions from seismic recordings. The geometric mean of the horizontal components of the Fourier Amplitude Spectra (FAS) of the dominant S-wave in the acceleration time window are smoothed by the Konno and Ohmachi55 algorithm. Then the GIT inversion decomposes the smoothed FAS Aij at station j for event i at source-station distance rij into source, path, and site components as
log Aij(f, rij) = log e S(f) + log P(f, rij) + log G(f) + δWij(f)
where e S correspond to the source term, P is the attenuation function comprising the geometrical spreading and frequency-dependent attenuation, G is the site response, and δWij is the remaining residual term from the GIT spectral fitting. P is considered 1 at distance r = rref = 10 km, which is also the minimum source-station distance permitted in the dataset. One or more stations are carefully selected as reference rock stations56, for which a flat and non-amplified site response is assumed (i.e., Gj = 1). The site response functions at other stations are thus derived relatively to this reference set.

The GIT decomposition was performed over many crustal events and stations in Central Italy42 (see also Data availability), considering amplitude spectra at 56 frequencies equally spaced on a logarithmic scale in the 0.5-25 Hz range. The solution of the system described by Eq. (5) was calculated following a non-parametric inversion scheme, i.e., no specific functional forms related to the source and propagation terms were prescribed. Therefore, the non-parametric inversion scheme permits capturing possible complexities not only in the path and site effect terms but also in the earthquake source term.

Here, we utilize the GIT apparent source term defined as log Sij = log e S(f) + δWij. Thus, Sij corresponds to a station-specific S-wave spectrum recorded at a reference rock station located at the reference distance rref. To obtain the apparent source spectrum ΩGIT(f), one must correct for the residual radiation, wave propagation effects in a half-space from the source to the station at reference distance, and for the attenuation of the shallow layers above the half-space. We assume
Sij(f) = ΩGIT(f) * FSij / (4πρβ^3 * rref) * Ge(f)
where ρ = 3.1 g cm-3, β = 3.5 km s-1 are the density and shear wave speed at the source, respectively. FSij is the radiation pattern considering the same combination of horizontal components as in GIT (i.e., geometric mean). Regarding the S-wave radiation pattern, we assume that wave scattering due to the random crustal heterogeneities disturbs the radiation pattern. We substitute this effect by calculating an average radiation pattern FSij over random perturbations of the source mechanisms by ±30° in strike, dip, and rake. Term Ge(f) accounts for the source-to-surface amplification of the reference rock site. In general, the rock sites experience site amplification with respect to the homogeneous full space due to the free surface effect and the impedance contrast between the shallow S-wave velocity-density profile and the velocity-density at the source region. We assume Ge(f) comprises the free surface effect (factor of 2), Boore and Joyner57’s generic rock site amplification, and shallow attenuation in the form of the high-frequency kappa effect with κ = 0.015 s corresponding to rock site conditions in Central Italy42,56. We note that overall, Ge(f) varies only slightly between 3.4 at 0.5 Hz and 2.6 at 25 Hz, with a maximum of 4.7 at 4 Hz (see Supplementary Fig. S4). After correcting Sij for the terms in Eq. (6), the resulting ΩGIT(f) represent the apparent spectra (AS) that we use as data in the dynamic inversion.

We point out that the AS includes the term δWij from GIT, which encompasses residual contributions from source, path and site-effect components. Nevertheless, the source radiation exhibits consistent, systematic spatial patterns across stations, producing a coherent signal within the AS. Therefore, when using a large number (~100 in our study) of stations in the dynamic rupture inversion, site and path residual effects primarily appear as unexplained posterior data errors and have only a minor impact on the inferred source models (see Supplementary Text S2).

We utilize the GIT-derived AS of two example events that are part of a high-quality catalog related to the latest major seismic sequence in the Central Italy area, which was densely recorded within 10–120 km and with a magnitude range of 3.2–4.5. The extensive collection was made possible thanks to the contributions from the stations of the temporary networks installed by the Italian Department of Civil Protection (DPC), as well as by several academic and research institutions, including the Istituto Nazionale di Geofisica e Vulcanologia (INGV); for more details, see Morasca et al.42, Colavitti et al.14 and Spallarossa et al.58.

The non-parametric GIT decomposition is implemented using the open-source object-oriented Python package GITpy59. The AS are retrieved at ~100 stations (Fig. 1a). Hypocentral parameters are adopted from the Engineering Strong Motion database (see Data Availability) and listed in Table 1. The centroid moment tensors were calculated using BayesISOLA60, and the fault planes were adopted from the resulting mechanisms. Of the two focal planes, we select the steeper ones that correspond to the orientation of the main tectonic structures in Central Italy. Nevertheless, we tried to calculate the AS for both planes, and the corresponding AS were almost the same, with a variance reduction of 0.999. This suggests that our method cannot reliably discriminate the fault plane for these two events. For more discussion, see also Supplementary Text S3.

Dynamic rupture simulation and synthetic AS calculation
As a numerical solver, we utilize open-source code FD3D_TSN designed to simulate 3D dynamic earthquake ruptures61. It employs a high-accuracy, 4th-order finite-difference (FD) method on a Cartesian grid to solve the wave equation in a layered Earth model containing a vertical fault located on one face of the computational domain. The code implements a traction-at-split-node method to model fault friction and utilizes perfectly-matched-layers absorbing boundaries to minimize reflections from the artificial domain margins. FD3D_TSN employs OpenACC directives of the Nvidia Fortran compiler for GPU acceleration for improved computational efficiency, achieving significant (~10x) speedups compared to CPU execution. The code has been validated against established SCEC benchmarks62, ensuring the code’s accuracy in simulating ruptures with heterogeneous material properties. See more details in Premus et al.61.

We assume slip-weakening friction law with spatially variable dynamic parameters. Namely, we invert for prestress above the dynamic friction, characteristic slip-weakening distance Dc, and so-called friction drop (difference between the static and dynamic friction), assuming a constant dynamic friction coefficient. The parameters are defined on a grid of control points 125 m apart, from which they are bilinearly interpolated onto the finer FD grid. The grid size was chosen so that the reliably resolved maximal simulated frequency correspond to the maximum frequency of the employed data, 25 Hz (compare with, e.g., Valentová et al.63). For more details, see e.g., Gallovič et al.23,64. Normal stress is assumed to be 100 MPa and constant over the fault plane. We prescribe a weak nucleation with a radius of 0.2 km and a maximum mean overstress of 1 MPa (Table 1) to avoid a strong effect on the synthetic AS at high frequencies.

Since the inverted events are small and relatively deep, we consider a dynamic rupture on a planar fault 2.5 × 2.5 km2 embedded in a homogeneous space (see Table 1). The FD simulations are challenged by the maximum considered frequency of 25 Hz. As listed in Table 1, the discretization of the simulation domain for the FD simulations is very dense (10 m) with dense temporal sampling (0.4 ms). Still, the dynamic rupture simulation takes less than 10 seconds at Nvidia GTX 4070 graphics card. We tested the numerical convergence of the dynamic rupture simulation by decreasing the spatial and temporal grid size to half (5 m, 0.2 ms) for the MAP models. The models differ by 2–4% in their AS misfit, which is negligible considering the frequency range up to 25 Hz.

After simulating the rupture for a given distribution of dynamic parameters (event i), slip rates Δu_i(ξ, t) are shifted in time following the S-wave propagation in a homogeneous space and integrated over the fault plane S to obtain apparent source time functions (ASTF) for station j,
Ψsynt(t, xj) = μ * integral(Δu_i(ξ, t - |xj - ξ| / β) dξ)
where μ is the shear modulus, xj is the j-th receiver location and ξ denotes the position of the slip rate on the fault. The synthetic AS is then obtained by calculating the Fourier amplitude spectrum of the second-time derivative of Eq. (7), Ωsynt(f) = |F(Ψ€synt)|, where F is the Fourier transform and |:| is the absolute value. As a final step, we apply the same smoothing to the synthetic AS as used in the empirical AS and evaluate it at the identical set of 56 frequencies in the same range of 0.5–25 Hz.

Model parameterization and Bayesian inversion
We use Bayesian formalism for the inverse problem, where the solution is defined as a posterior probability density function (PDF) combining likelihood and prior PDFs64. The misfit defining the likelihood function is assumed to be the L2 norm between the natural logarithms of empirical and synthetic AS values, normalized by the square of data error equal to 4; for more information on the error estimate, see Supplementary Text S2. The prior is assumed to comprise the following constraints: (i) Seismic scalar moment M0 being equal to the value from centroid moment tensor inversion (see Table 1) with Mw uncertainty ±0.05. This can also be interpreted as adding a zero-frequency amplitude to the data. (ii) To ensure a weak nucleation, we assume that the average overstress cannot exceed 1 MPa within a 200 m radius around the hypocenter (Table 1). (iii) The friction drop and prestress are non-negative and limited, and iv) Dc must be larger than 1 mm for the sake of the numerical precision.

For the posterior sampling, we utilize the parallel tempering method65, based on the Markov chain Monte Carlo (MCMC) approach. In this method, the MCMC exploration is performed in parallel at various temperatures, with the possibility of swapping temperatures based on an appropriately adapted Metropolis-Hastings criterion, allowing the exploration of larger areas of parameter space.

The parallel inversion is started from a suite of simple/smooth starting models with elliptical or circular rupture with various rupture speeds and stress drops. First, we created a dynamic rupture model with constant prestress and strength within an ellipse. Outside this region, the prestress is assumed to be half and the strength twice as high than on the inside, ensuring a gradual arrest of the rupture. Slip-weakening distance Dc is considered to increase linearly from the hypocenter at a prescribed rate. This ensures that the nucleation is weak, almost negligible in the ASTF and AS, and the rupture propagates at a constant speed (see below).

We introduce variability into the initial smooth models by multiplying the prestress by a random factor drawn from a Gaussian distribution with a standard deviation of 2. To keep the rupture of these perturbed models propagating and the seismic moment approximately unchanged, we adjust the dimensions of the source model and the Dc rate from nucleation, following scaling relations derived below. We additionally change the Dc rate randomly within ±10% to introduce variability in the rupture speed.

Supplementary Fig. S2 shows two representatives of smooth models of the directive event, one with symmetric and one with asymmetric rupture propagation. The rupture gradually arrests at the edges of the elliptical source region due to the prescribed weak barrier represented by decreased initial stress and increased strength as described above. Both models underestimate the observed AS at frequencies >~5 Hz. Our symmetrical model does not capture the azimuthal dependence of the observed AS. Interestingly, this holds even for the asymmetrical model, confirming the sensitivity of the synthetic AS to the rupture evolution. Note that prescribing a stronger barrier with a fast rupture propagation would result in omega-squared spectra due to the radiation from the sudden stop of the rupture in some directions, as shown in the classical models of Sato and Hirasawa66 and Madariaga8 and extensively studied by Kaneko and Shearer19. Supplementary Text S4 describes the performance of several simulated smooth models surrounded by the strong barrier, suggesting that the rupture must be relatively fast to generate enough high-frequency radiation at its arrest, unlike the heterogeneous models with lower rupture speeds. Although such parameterization warrants further testing against real observations, it would require additional data, e.g., data sensitive to rupture velocity, which is beyond the scope of this work. Here, we employ parameterization with generally heterogeneous dynamic parameters, which naturally incorporates multiple deceleration and acceleration episodes at various scales and is thus considered more realistic when modeling real data (see Results).

To test the resolving power of the Bayesian dynamic source inversion of the amplitude AS, we designed a synthetic test based on the directive Mw 4.2 event (Supplementary Text S1). The target model used to generate AS data for the same set of stations is a directive model with additional random heterogeneities in dynamic parameters for increased complexity. The MCMC inversion is run with the same settings as for the real event. For a detailed description, see Supplementary Text S1.

For each of the real Mw ~ 4 earthquakes, we performed the parallel tempering sampling of the posterior PDF in 6 parallel jobs on 3 GPU Nvidia RTX 4070 cards, each job with 8 chains with 6 random temperatures T between 1 and 100 and 2 temperatures T = 1. Each chain started from a randomly altered initial model. The entire inversion ran for >30 days, during which each chain performed more than 5000 steps, i.e., about 250k models were visited. The posterior analysis is conducted from samples at T = 1 after discarding 20% of the initial steps (burn-in phase) and keeping only every 10th accepted model to manage the storage size of the output models.

Generation of initial elliptical models
To start the MCMC random walk, we use a set of various randomized smooth initial models. We note that most of the model space is filled with models that fail to even nucleate rupture. Our strategy is to generate an ensemble of reasonable models that are characterized by a priori known (approximately) seismic moment M0.

We first create a dynamic rupture model that fits the assumed seismic moment and the low-frequency range of the observed apparent source spectra. Its rupture characteristics stem from the following equations. A simple circular crack model predicts the dependence between the seismic moment, constant stress drop Δσ, and radius a,
M0 = (16/7) * Δσ * a^3

We assume standard values of stress drop (~3 MPa for Central Italy, e.g., Pacor et al.16), and use Eq. (8) to infer the corresponding radius of the source. The nucleation point is set either in the center if no directivity is observed in the empirical spectra or shifted by trial-and-error to either side of the rupture. For characteristic slip-weakening distance Dc, we assume it increases with hypocentral distance ρ (in km) from minimum value Dc0 at rate Dcrate,
Dc = Dc0 + Dcrate * ρ

This way, the nucleation can be kept small, and the rupture propagates at approximately constant speed. In the next section, we derive the scaling relation of Dcrate with Δσ and strength drop Δτs. It reads
Dcrate = (π * Δσ^2 / (μ * Δτs)) * ec(v)

where ec(v) is a non-dimensional parameter between 0.1-1, depending on rupture velocity v, rupture mode and shape, etc. (see below). We found that in our case of having nucleation sizes of diameter ~0.1 km, a reasonable choice of Dc0 is its value predicted by the Dcrate value at the edge of the nucleation, i.e., Dc0 = Dcrate * 0.05 km. This defines an initial rupture model. We then manually adjust and test these parameters to obtain a reasonable fit up to about twice the corner frequency. Note that resulting rupture models usually exhibit strong decay above their corner frequencies due to the smoothness of the dynamic parameters and gradual rupture arrest33, see also Supplementary Fig. S2. This is later remedied by the inversion, permitting small-scale variations in the dynamic parameters.

Scaling of the distance-dependence of Dc
We assume an expanding self-similar anti-plane (mode III) 1D crack with constant speed v, for which Broberg67 derived the stress intensity factor at the crack tip,
KIII(v) = Δσ * sqrt(π * v * t) * E(γ(v))

where Δσ is the stress drop, γ(v) = sqrt(1 - (v/β)^2) and E is the complete elliptical integral of the second kind. For such a model, the dynamic energy release rate reads
GIII(v) = KIII^2(v) / (2μ * γ(v))

Combining Eqs. (11) and (13) yields
GIII(v) = (Δσ^2 / (2μ)) * π * v * t * (E^2(γ(v)) / γ(v)) = (Δσ^2 / (2μ)) * π * v * t * c(v), with c(v) = E^2(γ(v)) / γ(v)

Assuming the rupture propagating over a fault with a given fracture energy rate G will reach a steady-state velocity v after extending to a size a = vt, from Eq. (14) it follows:
GIII(v) = (Δσ^2 / (2μ)) * c(v) * a

Considering a rupture governed by the slip-weakening friction, the fracture energy rate is
G = (1/2) * Dc * Δτs

where Dc is the slip-weakening distance and Δτs is the strength drop. Equating the fracture and dissipation energy rates (Eqs. 15 and 16) yields a condition for the increase of Dc with the increase of crack size a,
Dc = (π * Δσ^2 / (μ * Δτs)) * c(v) * a

Figure 9 shows the dependence of c(v) on the ratio of v/β and how it can be well approximated by a linear function for a wide range of reasonable rupture velocities.

Numerical experiments with circular and elliptical ruptures suggest that the Dc rate must be 1.3-2.2 times smaller than predicted by Eq. (17), perhaps due to the transition from 2D to 3D rupture model. Therefore, ec in Eq. (10) is 1.3-2.2 times smaller than c. Equation (17) is used to define the Dc rate for the initial rupture models. For our rupture velocity (constant in the initial model, but variable in the sustained), min and max v obtained assuming 2.2-1.3 times larger values than predicted by Eq. (17), respectively, should fall within 1-3 km s-1.

Data availability
The waveform data, locations and source metadata are provided through the website of the Engineering Strong-Motion (ESM68,69; https://esm-db.eu/) and the Italian ACcelerometric Archive (ITACA70), https://itaca.mi.ingv.it/ under IDs 20161030_0000130 (Mw 4.2) and 20160824_0000007 (Mw 4.5). The FAS dataset is available from https://shake.mi.ingv.it/central-italy/. The data and results of the dynamic inversions presented in this work are available at Zenodo71 (https://doi.org/10.5281/zenodo.17417423).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Gallovic et al. employ a novel dynamic rupture inversion technique to investigate the rupture processes of two magnitude 4 earthquakes in central Italy. Their key conclusions are that small-scale (~100 m) heterogeneities in rupture propagation are essential to fit the high-frequency radiation and that fractally distributed dynamic parameters (prestress, friction drop, characteristic slip-weakening distance, etc.) are required to explain the data. While previous kinematic inversion studies and theoretical work have yielded similar insights, this study is the first to apply dynamic rupture inversion to reach the same conclusion, as the authors correctly claim. Although the results are potentially significant and the inversion's ability to accurately resolve only the average properties of source parameters is reasonable and useful for understanding the source processes of small earthquakes, I believe additional simulations and clarifications are necessary to support the main conclusions (as stated below).\nFurthermore, this study analyzes only two small events, limiting the generalizability of the inferred source characteristics. As the authors acknowledge, further research is needed to determine whether these findings extend to other small earthquakes. However, I question the validity of the inferred fractalness, where pre-stress, strength excess, and Dc are all fractally distributed in the two analyzed events. While this might be valid for these specific cases, specific fractal distributions obtained in this study, particularly for the pre-stress and strength excess, cannot necessarily be extrapolated to either tiny or large earthquakes. If self-similarity holds, these fractal distributions would not be generalizable to other small earthquakes.\nMajor comments:\nThe authors state (Line 170) “We would like to emphasize that the smooth average dynamic model would be unable to fit the observed AS data (see Materials and Methods and cf. Fig. S2).” This is an important result and should be discussed in more detail. I recommend moving Fig. S2 into the main manuscript. More importantly, additional synthetic tests and alternative parameterizations of dynamic rupture models are needed to support this claim. Classical models by Sato and Hirasawa (1973) and Madariaga (1976), among others, incorporate a sudden rupture arrest, which results in an omega-squared high-frequency decay. However, Figure S2 suggests that the rupture gradually arrests at the edge of the source region, likely due to the assumptions made in setting up the “smooth” rupture model. For instance, Figure S2 indicates that, just outside the slip zone, the initial stress is only slightly smaller and the static minus dynamic friction coefficient is relatively large. This setup does not produce a sudden rupture arrest, which is typically required for high-frequency energy generation. If the authors were to implement a sudden arrest model, the results might differ significantly from those presented in Figures S2a and S2d. Given the importance of this conclusion, additional results and discussion are warranted.\nThe mean spectral bias between observed data and the MAP results (Figure 3a) indicates that low-frequency (<0.8 Hz) energy is systematically missing in the MAP model. Since the vertical axis is logarithmic, the discrepancy reaches approximately 30% at 0.5 Hz, which is substantial given that this is the most resolvable frequency range. Please provide an explanation for this discrepancy and justify the difference.\nThe authors state (Line 223) “We fit Eq. (1) using L2 minimization in logarithms to the synthetic MR spectra from our ensemble, considering scalar moment M0 as known.” Please include the fitting results, as they are crucial for assessing how well Eq. (1) fits across a broad frequency range. Additionally, are there any biases between estimated fc and the goodness of fit?\nThe authors state (Line 255) “This seems to contradict the findings of Yoshida (2019), who showed that stress drop estimates can be almost twice as low when the asymmetry of the rupture propagation is neglected.” More context is needed here, as the line of reasoning is unclear. First, what is meant by “the theoretical model of Yoshida (2019)” (Line 259)? Yoshida (2019) used the model of Dong and Papageorgiou (2003) rather than developing a new theoretical model. Furthermore, Kaneko and Shearer (2015) showed that their asymmetrical rupture model is set up to produce a longer source duration than the symmetrical model, reflecting the smaller k values in the asymmetrical model. Since the present study infers similar source durations for the two events, it is inappropriate to interpret the results using the Kaneko and Shearer (2015) models. Please provide more background before comparing the results to previous theoretical models.\nThe authors state (in Line 293) that the inferred k^-2 decay of PSD corresponds to k^-1 decay of the amplitude spectrum, and such a decay was considered for the initial stress in previous theoretical studies. However, the previous studies (e.g., Ripperger et al, 2007) assumed that stress heterogeneity is described by a power law decay of its wave number spectrum. I might be missing something here, but Figure 8 shows increasing amplitude with wavelength, and hence, the trend is opposite from those assumed in the previous studies.\nMinor comments:\nLine 19: The statement “Heterogeneity in the rupture propagation down to the smallest scales is necessary to fit the high-frequency radiation” could be misleading, as it does not consider the debate on whether high-frequency energy arises from elastic heterogeneities near the source rather than rupture propagation. Since this study does not account for elastic heterogeneities, the statement should be revised to avoid confusion.\nLine 35: Yoshida (2024) should be corrected to Yoshida (2019).\nLine 45: The phrase “can attain different than” is unclear. Please clarify.\nLine 73: The claim that this study performs “an inversion of AS free of path and site effects” seems overstated. Complete removal of path and site effects is generally impossible, particularly due to elastic heterogeneities near the source (e.g., fault damage zones within a few hundred meters). Consider revising this statement.\nLine 90: The phrase “described in Text S1” is inappropriate for the Introduction section, as it references results before they are presented. Move this sentence to the Results section.\nLine 145: “due to the prescribed weak nucleation”. This is not explained anywhere.\nLine 458: The authors should discuss numerical resolution in terms of cohesive zone sizes (e.g., Day et al., 2005), as is common in dynamic rupture simulation studies.", "reviewer_2": "I have reviewed the code at a high level. It seems in good order and documented. However, I did not have tome to actually run and test it. For a code at this complexity, that would take a time commitment of many hours.", "reviewer_3": "I read the paper “Physics-based broadband characterisation of weak earthquakes” by Gallic et al. with considerable interest. It is a fine piece of work and appears to be of high quality and completeness. It shows how spatial variations in earthquake amplitude spectra can be used to determine important details about the rupture properties of even relatively small earthquakes, in this case of magnitudes 4.2 and 4.5. It seems remarkable to see rupture directivity resolved for an earthquake as small as magnitude 4.2 Moreover it shows how failing to account for these details can sometimes lead to very significant bias in the conventional methods for estimating important parameters like stress drop. I think the paper is an important step forward in estimating earthquake source parameters.\nHowever, while I myself found the paper very interesting I would have to say that I think it appeals mainly to a narrow audience of ground motion seismologists. If Nature is fine with accepting a paper that may not have broad appeal, then the paper can be published with only minor revision. If instead it is Nature’s goal to publish only papers with broad appeal, then this paper may not be such a good fit. In order to make it appeal to a wider audience I think it would require major revision, such as perhaps considering a much larger dataset - difficult, given that these two events each took 30 days of run time.\n1) My main substantive comment is that in the section on “GIT derived apparent spectra”, some important assumptions are mentioned, such as the value of the high-frequency attenuation coefficient kappa (=0.015, Line 411), and the use of random perturbations in the strike, dip, and rake to mimic random perturbations to crustal velocity structure. These all may be valid assumptions, but without further evidence I’m left wondering whether there could be trade-offs between some of the source properties and the kappa value. Random crustal velocity perturbations will not only obscure the S-wave radiation pattern but can also produce frequency dependent apparent attenuation due to scattering. Could such effects be mapped into the source properties if they are not accounted for? In any case, I feel more discussion of these is warranted.\n2) As with all Bayesian inversions, an important question is: Has it converged? With the sampling of 250k models over 30 days run time, I could well believe the solution has converged. However, some indication of how the authors assessed this, such as examining the character of the chain histories of likelihood values, should be discussed.\nLine 519: Shouldn’t the total number of models visited by 5k steps x 8 chains x 6 jobs = 240k, not 250k?\nFigure 1: It can be difficult to make out the stars that denote the epicentres. Could they be made bigger, or a different colour so they’d stand out more? The lines connecting the beachballs to the respective epicenters are also very thin.\nFigures 6 & 7: These are a little confusing because the usnits of the x-axis at the bottom of the figure are actually given as part of the figure title, at the top of the figure. I recommend moving the units to the bottom of the figure, just below the x-axis, as has been Doen with other figures."}, "author_response": "Dear Editor and Reviewers,\nWe thank you for your time and efforts in carefully reading our manuscript and for the\nprovided comments. We have addressed all of them as we describe below, and we believe this has\nimproved our manuscript considerably.\nReviewer #1:\nGallovic et al. employ a novel dynamic rupture inversion technique to investigate the rupture\nprocesses of two magnitude 4 earthquakes in central Italy. Their key conclusions are that small-\nscale (~100 m) heterogeneities in rupture propagation are essential to fit the high-frequency\nradiation and that fractally distributed dynamic parameters (prestress, friction drop, characteristic\nslip-weakening distance, etc.) are required to explain the data. While previous kinematic inversion\nstudies and theoretical work have yielded similar insights, this study is the first to apply dynamic\nrupture inversion to reach the same conclusion, as the authors correctly claim. Although the results\nare potentially significant and the inversion's ability to accurately resolve only the average\nproperties of source parameters is reasonable and useful for understanding the source processes of\nsmall earthquakes, I believe additional simulations and clarifications are necessary to support the\nmain conclusions (as stated below).\nFurthermore, this study analyzes only two small events, limiting the generalizability of the\ninferred source characteristics. As the authors acknowledge, further research is needed to\ndetermine whether these findings extend to other small earthquakes. However, I question the\nvalidity of the inferred fractalness, where pre-stress, strength excess, and Dc are all fractally\ndistributed in the two analyzed events. While this might be valid for these specific cases, specific\nfractal distributions obtained in this study, particularly for the pre-stress and strength excess,\ncannot necessarily be extrapolated to either tiny or large earthquakes. If self-similarity holds, these\nfractal distributions would not be generalizable to other small earthquakes.\nWe agree with the reviewer that the inferred fractal characteristics of the dynamic parameters\ncannot be immediately generalized to other events and magnitudes without further studies. For\nnow, we have found the fractal properties only of the two events, but more importantly, with\nsignificantly different shapes of their apparent spectra. We point out that the perturbations are\ndescribed “on top of” the mean spatial spectra of dynamic parameters, and do not describe the\nmean spectra themselves. Based on the theoretical description of von Karman power spectrum and\na few preliminary tests that we performed, the fractal properties can be rescaled both to larger and\nsmaller events and generate various spectral behaviors. However, it warrants a devoted, in-depth\nstudy. We have added such a note in Conclusions: “Future studies are needed to test whether the\nfractal variations in prestress, strength excess, and Dc, inferred for our two events with distinct\nAS shapes, apply to other events and magnitudes.”\nThe authors state (Line 170) “We would like to emphasize that the smooth average dynamic\nmodel would be unable to fit the observed AS data (see Materials and Methods and cf. Fig. S2).”\nThis is an important result and should be discussed in more detail. I recommend moving Fig. S2\ninto the main manuscript. More importantly, additional synthetic tests and alternative\nparameterizations of dynamic rupture models are needed to support this claim. Classical models\nby Sato and Hirasawa (1973) and Madariaga (1976), among others, incorporate a sudden rupture\narrest, which results in an omega-squared high-frequency decay. However, Figure S2 suggests that\nthe rupture gradually arrests at the edge of the source region, likely due to the assumptions made\nin setting up the “smooth” rupture model. For instance, Figure S2 indicates that, just outside the\nslip zone, the initial stress is only slightly smaller and the static minus dynamic friction coefficient\nis relatively large. This setup does not produce a sudden rupture arrest, which is typically required\nfor high-frequency energy generation. If the authors were to implement a sudden arrest model, the\nresults might differ significantly from those presented in Figures S2a and S2d. Given the\nimportance of this conclusion, additional results and discussion are warranted.\nWe agree that the statement was too strong and thus confusing. Strong stopping phases would\nindeed result in omega-squared radiation in some directions, as analyzed in detail by, e.g., Kaneko\nand Shearer (2015). The topic of whether the rupture is necessarily heterogeneous is still open and\nvery interesting, as it is difficult to extract from real observations.\nBased on the reviewer’s suggestion, we performed additional synthetic tests, in particular a\nsimple parametric study for the non-directive Mw4.5 event. In this case, the rupture geometry is\nsimple—circular with nucleation in the center of the fault. The initial stress and strength were\nprescribed to be homogeneous inside the circle (except for the increased initial stress in the small\nnucleation patch), and the barrier around the circle was prescribed with zero initial stress and 15\ntimes larger strength.\nThe parametric study was performed for two main parameters: stress drop and rupture\nvelocity. The stress drop Δσ (determined mainly by the initial stress) controls the size of the rupture\nwith radius 𝑟 based on the scaling relationship for a given seismic moment 𝑀 : 𝑀 = 16 𝛥𝜎𝑟3. We\n0 0\n7\ntested 3 values of Δσ: 2 MPa, 4 MPa, and 8 MPa. The second parameter, the rupture velocity, is\ndetermined by the prescribed Dc rate, describing the assumed linear increase of Dc from nucleation\naccording to Eq. (17) in the main text. We assumed several values of rupture velocity ranging\nbetween 0.4 and 1.0 of β.\nThis simple setting results in models with almost constant stress drop and rupture velocity\nalong the fault, i.e., 2-parameter space. For these barrier dynamic models, we compare the\nsynthetic apparent spectra with the observed ones. The result of this simple parametric study with\nmodels color-coded by misfit (left) or spectral bias (right) is shown in this figure:\non the resulting dynamic simulation and thus slightly differ from the original (prescribed) value.\nThe best model has a misfit of ~52. Based on this value, it might have been accepted in our\nensemble if it had been present in the MCMC sampling. We note that the best model of this\nparametric test should not be compared with the best model of the Bayesian inversion, as in the\nformer case, the parameters were not optimized. Nevertheless, based on this simple parametric\nstudy, we can infer some general characteristics of the barrier model with apparent spectra similar\nto the observations.\nProperties of the best barrier model:\nhave higher misfit (>70). Thus, the obtained stress drop value lies within the marginal posterior\ndistribution, as found by the Bayesian inversion (cf Fig. 6).\nContrary to the heterogeneous model, the barrier model requires a higher rupture speed (2.9\nkm/s, i.e., Vr/β~0.8), likely to generate high-frequency radiation spectra comparable with the\nobserved spectra, i.e., the rupture hits the barrier with high enough energy. For slower ruptures,\nthe spectra are weaker than the observed ones (negative bias), while the faster ruptures result in\ntoo strong spectral radiation (positive bias). In addition, there is a small (negative) trade-off\nbetween the Vr and the stress drop. However, in the heterogeneous model found by the Bayesian\ninversion, the high mean rupture speed was not necessary as the rupture includes small-scale\naccelerations and decelerations over the whole evolution and thus permits smaller mean rupture\nvelocities.\nWe note that the barrier and heterogeneous models differ in radiation efficiency (~0.2 and\n~0.3 for the heterogeneous and barrier models, respectively). Unlike the rupture velocity, this\n4\ndifference is not as distinctive, as the value 0.3 lies in the upper bound uncertainty of the marginal\nposterior PDF of the Bayesian inversion.\nAs noted by the reviewer, the smooth models presented in the supplement do not exhibit\nenough stopping phases due to the weak barrier and the consequent gradual rupture arrest. Based\non this simple parametric study, this can be improved only by models with ruptures propagating\nfast enough. On the other hand, generally heterogeneous source models have the capability of\nincorporating multiple deceleration and acceleration episodes at various scales without the mean\nrupture velocity being so high.\nThe barrier model can thus be considered as a complementary approach to our\nheterogeneous-up-to-small-scale parameterization, as we cannot discern the true rupture velocity\nbased on the amplitude spectral data. Although not a primary point, we believe that our study can\ncontribute to the discussion on whether small-scale heterogeneity is necessary. However,\naddressing this point properly requires an in-depth study including an exhaustive parametric study\nof the barrier model and additional data sensitive to rupture velocity, and thus, we keep it for future\nwork.\nWe keep the gradually arresting smooth model examples in the Supplement as example\nmodels corresponding to the ensemble-averaged smooth model. The new comment to Fig. S2 now\nreads: “In Fig. S2, we show two representatives of smooth models of the directive event, one with\nsymmetric and one with asymmetric rupture propagation. The rupture gradually arrests at the edges\nof the elliptical source region due to the prescribed weak barrier represented by decreased initial\nstress and increased strength (both 2 times with respect to the conditions within the rupture). Both\nmodels underestimate the AS at frequencies >~5 Hz. Our symmetrical model does not capture the\nazimuthal dependence of the AS. Interestingly, even the asymmetrical model exhibits unexplained\nazimuthal dependence, confirming the sensitivity of the synthetic AS to the rupture evolution. Note\nthat prescribing a stronger barrier with fast rupture propagation would result in omega-squared\nspectra due to the radiation from the sudden stop of the rupture in some directions, as shown in the\nclassical models of Sato and Hirasawa (1973) and Madariaga (1976) and extensively studied by\nKaneko and Shearer (2015). Although such parameterization warrants further testing against real\ndata, it would require additional data, e.g., data sensitive to rupture velocity, which is beyond the\nscope of this work. Here, we employ parameterization with generally heterogeneous dynamic\nparameters, which naturally incorporates multiple deceleration and acceleration episodes at\nvarious scales (see Results).”\nWe have also changed the section title to “Model parameterization and Bayesian inversion”.\nMoreover, we rewrote several sentences in the text. In Abstract, sentence “We find that the\nheterogeneity in the rupture propagation down to the smallest scales (~100 m) is necessary to fit\nthe high-frequency radiation.” has been updated to “In our models, the heterogeneity down to the\nsmallest scales (100m) is required to fit the observed high-frequency radiation.” In Conclusions,\nwe have changed text “A significant result is that the approximately omega-squared apparent\nsource spectral falloff can be explained only by heterogenous ruptures with fractal properties, in\ncontradiction to smooth models with faster spectral decay.”, for “The approximately omega-\n5\nsquared apparent source spectral falloff observed empirically can be explained by heterogeneous\nruptures with fractal properties.”\nThe mean spectral bias between observed data and the MAP results (Figure 3a) indicates that\nlow-frequency (<0.8 Hz) energy is systematically missing in the MAP model. Since the vertical\naxis is logarithmic, the discrepancy reaches approximately 30% at 0.5 Hz, which is substantial\ngiven that this is the most resolvable frequency range. Please provide an explanation for this\ndiscrepancy and justify the difference.\nThe underestimation is below the corner frequency and thus related to the seismic moment\nM0. Below, we show the bias for a model from our ensemble with the highest M0 (3.5e15 Nm vs.\n2.0e15 Nm of the MAP model shown in Fig. 3a), demonstrating the consequent overestimation at\n0.8 Hz by a similar amount (~30%). The corresponding Mw values are 4.3 and 4.13 for the largest\nand the best model, respectively. Recalling that the prior constraint was 4.2±0.1, the models fit the\nprior constraint similarly well. This apparent discrepancy (similarly to other frequencies) is\ncovered by the estimated uncertainty of the inverted models. We have added a note to the end of\nthe Best rupture model section: “Note that the ~30% underestimation at the lowest frequencies is\nrelated to the seismic moment M_0; models with higher M_0 overestimate the data by a similar\namount in that frequency range.”\nThe authors state (Line 223) “We fit Eq. (1) using L2 minimization in logarithms to the\nsynthetic MR spectra from our ensemble, considering scalar moment M0 as known.” Please\ninclude the fitting results, as they are crucial for assessing how well Eq. (1) fits across a broad\nfrequency range. Additionally, are there any biases between estimated fc and the goodness of fit?\n(Fig. 4a) and the fitted Brune spectrum:\nThick black curve corresponds to the MAP model. Thick and thin green curves are the\nensemble average and one standard deviation. We have added the plots as new Fig. S3.\nNote that while the directive Mw4.2 event has an omega-squared source spectrum, the bias\nis close to zero, suggesting the omega-squared spectrum is a suitable approximation. Conversely,\nthe spectrum of the nondirective Mw4.5 event decays faster, and thus, the omega-squared spectrum\noverestimates the MR spectrum up to ~2 Hz, whereas it underestimates it at higher frequencies.\nThis is indeed reflected in the corner-frequency vs. goodness-of-fit plot shown below, where the\nmisfit is generally higher for the non-directive Mw4.5 event. Despite this, the corner frequencies\nand the Madariaga stress drop estimates result in correct estimates.\nThe following figures show the estimated corner frequency vs M0, shaded by misfit (L2\nnorm in ln):\nExcept for the possible trade-off between fc and M0, there is no clear trade-off between fc\nand misfit. However, the lighter misfit colors clearly demonstrate the generally larger misfit values\nof the nondirective Mw4.5 event.\nslightly worse for the nondirective Mw4.5 event but only above ~3 Hz (i.e., well above the corner\nfrequency) because its spectral decay is faster than the omega-squared.”\nThe authors state (Line 255) “This seems to contradict the findings of Yoshida (2019), who\nshowed that stress drop estimates can be almost twice as low when the asymmetry of the rupture\npropagation is neglected.” More context is needed here, as the line of reasoning is unclear. First,\nwhat is meant by “the theoretical model of Yoshida (2019)” (Line 259)? Yoshida (2019) used the\nmodel of Dong and Papageorgiou (2003) rather than developing a new theoretical model.\nFurthermore, Kaneko and Shearer (2015) showed that their asymmetrical rupture model is set up\nto produce a longer source duration than the symmetrical model, reflecting the smaller k values in\nthe asymmetrical model. Since the present study infers similar source durations for the two events,\nit is inappropriate to interpret the results using the Kaneko and Shearer (2015) models. Please\nprovide more background before comparing the results to previous theoretical models.\nThe paragraph mentioned by the reviewer is primarily about the determination of the corner\nfrequency and the stress drop inferred from it, which depends on the value of k. The preceding\nparagraph compares stress drops estimated from the corner frequency inferred from the moment\nrate spectra of the dynamic models. In that case, the theoretical models (Madariaga) provide\ncorrect stress drops. However, when the corner frequency is estimated directly from the apparent\nsource spectra (AS), it is biased due to the averaging over the stations. Indeed, the apparent source\ntime function (and thus AS) corresponds theoretically to the moment rate function only for a station\nlocated perpendicular to the fault at an infinite distance. For normal faults, such a condition is\ndifficult to meet in a real station configuration. To further demonstrate this, we show in the figure\nbelow the comparison of the apparent source time functions (left) and spectra (right, in violet) at\nall our stations (gray curves) with the moment-rate characteristics for a symmetrical barrier model\n(as shown in response to Q2 above).\nAs there are no stations that would be exactly perpendicular to the fault, the average observed\nspectrum (thick green) is very different from the moment rate spectrum (black). Therefore, it is\ntrue that our examples should not be compared with the results of Kaneko and Shearer (2015),\n8\nbecause they performed station averaging over the full focal sphere, which might reduce the\ndiscrepancy of the stress drop estimates. In our manuscript, we tried to emphasize that the problem\narises only if the corner frequency is estimated directly from the azimuthally-limited data (AS)\ninstead of the moment rate spectra. However, the moment rate can be acquired from the data only\nby the final-fault inversion. Note also that although the source duration of the two events is similar,\nthe events differ in magnitude, and thus, the relatively longer duration of the smaller event can be\nassigned to the rupture asymmetry, consistent with Kaneko and Shearer.\nWe are sorry for our confusing discussion that mixes the various corner frequency estimates\nand their applicability based on the symmetry of the rupture. We have simplified some parts of the\nsection and the Fig. 7 caption to distinguish direct and indirect corner frequency (and stress drop)\nestimates. As suggested by the reviewer, we have also revised the comparison with the published\nworks. Hopefully, the text is now clearer and more precise:\n“We also calculate the Madariaga stress drop using the corner frequencies inferred directly\nby fitting the GIT AS to the omega-square model (shown as arrows in Fig. 7). The stress drop is\noverestimated significantly for the directive Mw4.2 event, resulting from the overestimated corner\nfrequency. We attribute this to the directivity of the rupture model and the insufficient averaging\nover the focal sphere (Kaneko and Shearer, 2015). This suggests that the simple spectral stress\ndrop estimates from averaged AS (provided physics-based models, such as Madariaga, are\nutilized) may be biased in practice, especially for events with strong azimuthal variability of the\nAS. Therefore, the practitioners should first examine the variability of the AS shapes as their large\nspread may indicate a directive event for which the standardly inferred source properties might be\nsignificantly biased (Fig. 4a). A potential improvement may be in such cases represented by fitting\napparent (station-dependent) corner frequencies estimated from the AS using asymmetrical rupture\nmodels (e.g., Yoshida, 2019). However, the most reliable approach, regardless of the event’s\ndirectivity, is to derive directly the moment rate functions by, e.g., dynamic or kinematic finite-\nfault inversion (Gallovič et al., 2025).”\nThe authors state (in Line 293) that the inferred k^-2 decay of PSD corresponds to k^-1 decay\nof the amplitude spectrum, and such a decay was considered for the initial stress in previous\ntheoretical studies. However, the previous studies (e.g., Ripperger et al, 2007) assumed that stress\nheterogeneity is described by a power law decay of its wave number spectrum. I might be missing\nsomething here, but Figure 8 shows increasing amplitude with wavelength, and hence, the trend is\nopposite from those assumed in the previous studies.\nThe inferred spectral decay is indeed k^-2, i.e., the spectral amplitudes decrease with\nincreasing wavenumber (=decreasing wavelength). The confusion probably stems from plotting\nthe x-axis (wavelengths) in Fig. 8 in a reverse manner. To avoid the confusion, we have added a\nnote that the wavelengths are shown “plotted as a function of wavelength (the inverse of\nwavenumber k), displayed from long to short wavelengths (i.e., decreasing from left to right)”. We\nhave also added the term “power law decay” to the Fig. 8 caption and text where appropriate.\nLine 19: The statement “Heterogeneity in the rupture propagation down to the smallest scales\nis necessary to fit the high-frequency radiation” could be misleading, as it does not consider the\ndebate on whether high-frequency energy arises from elastic heterogeneities near the source rather\nthan rupture propagation. Since this study does not account for elastic heterogeneities, the\nstatement should be revised to avoid confusion.\nWe have corrected the statement to “rupture heterogeneity” to avoid confusion.\nLine 35: Yoshida (2024) should be corrected to Yoshida (2019).\nWe are sorry for the wrong reference. The correct one is Yoshida (2024) in SRL. We have\ncorrected it in the reference list.\nLine 45: The phrase “can attain different than” is unclear. Please clarify.\nWe have modified the sentence to “the empirical spectra exhibit decay at high frequencies\nthat is different from the omega-squared decay”.\nLine 73: The claim that this study performs “an inversion of AS free of path and site effects”\nseems overstated. Complete removal of path and site effects is generally impossible, particularly\ndue to elastic heterogeneities near the source (e.g., fault damage zones within a few hundred\nmeters). Consider revising this statement.\nWe agree and have corrected the expression to be more appropriate: “AS corrected for path\nand site effects.”\nLine 90: The phrase “described in Text S1” is inappropriate for the Introduction section, as\nit references results before they are presented. Move this sentence to the Results section.\nWe have removed the phrase from the Introduction, keeping the first reference to Text S1 in\nResults.\nLine 145: “due to the prescribed weak nucleation”. This is not explained anywhere.\nWe agree that the nucleation specified only in Tab. 1 was insufficient. Therefore, we have\nadded a text specifying the nucleation conditions and reference to Table 1 to the Methods section:\n“We prescribe a weak nucleation with a radius of 0.2 km and a maximum mean overstress of 1\nMPa (Tab. 1) to avoid a strong effect on the synthetic AS at high frequencies.”\n(e.g., Day et al., 2005), as is common in dynamic rupture simulation studies.\nSince the models are strongly heterogeneous, cohesive zone size varies, and their estimates\ndo not have to be precise. Nevertheless, we have tested the modeling by densifying the FD spatial\n(and temporal) grid steps to half (5 m, 0.2 ms). For the MAP models, we provide a comparison of\nthe MR spectra and modeling bias below:\nMoment jerk spectra Original sampling Finer sampling\nMw4.2:\nMw4.5:\nThe misfit value changes only by 2-4% for the finer sampling, which is negligible\nconsidering the accepted models’ misfits vary by up to 100%. The dynamic simulations are thus\nnumerically robust. We have added the following text to Methods (section Dynamic rupture\nsimulation and synthetic AS calculation): “We tested the numerical convergence of the dynamic\nrupture simulation by decreasing the spatial and temporal grid size to half (5 m, 0.2 ms) for the\nMAP models. The models differ by 2-4% in their AS misfit, which is negligible considering the\nfrequency range up to 25 Hz.”\nThis manuscript considers dynamic rupture inversion for earthquakes of moment magnitude\n4. Dynamic rupture models are essentially a research frontier in understanding the physical\nprocesses of earthquake rupture. These models permit more detailed insight into physical processes\nthan kinematic models, spectral models, or point source models. For example, the physical\nconstraints/correlations build into these models often ensure more physically meaningful results.\nTo my knowledge, only few works exist that characterize sources with kinematic models at these\nmagnitudes. After the recent success of dynamic rupture studies, extension to lower magnitudes is\nquite desirable.\nThe authors choose a probabilistic framework for their new method which considers\napparent-spectra data. In my view, this is attractive since (1) it permits intrinsic\nuncertainty quantification, (2) permits general parametrizations since the method is nonlinear/numerical in\nnature, and (3) in some ways future proofs the method since it can be extended because of (1) and\n(2). There is significant new knowledge to be gained if the methods presented in this manuscript\nwere to be applied more broadly.\nI think this is an excellent manuscript and should be of interest to the entire earthquake\ncommunity and likely beyond. I recommend considering publishing it after addressing the\ncomments below.\nSpecific comments:\n1) The method employs AS for data, and these are “corrected” for certain effects by an\ninversion process (termed GIT). Since Bayesian methods generally rely on some form of data error\nestimation or assumption. However, regardless of the assumptions made, the errors need to be\npreserved in some meaningful way by the pre-processing, here including GIT. Is my understanding\ncorrect? Could the authors comment on this or explain in a bit more detail which impacts the GIT\nhas on the overall results?\nGIT employs an L2 minimization of the smoothed Fourier Amplitude Spectra into individual\nterms (site, path, and source). It does so under the assumption of uncorrelated, homoscedastic\nerrors, but without explicitly accounting for measurement uncertainty. In principle, it could use an\nestimate of the data covariance matrix to obtain the posterior uncertainty of individual terms.\nHowever, to the best of our knowledge, the current implementation of GIT, including its Python-\nbased version (GITpy; available at https://gitlab.rm.ingv.it/inversion/gitpy), does not incorporate\na data covariance matrix in the inversion process. Therefore, it does not provide formal posterior\nuncertainty estimates for the individual inverted terms, including the source spectral term. In\naddition, our data also include the GIT residual term, which contains all (non-systematic) effects\nof source (e.g., directivity), path (e.g., 3D structure), or site (e.g., site effect anisotropy). Although\nthe source effects seem to dominate this term (as seen for the directive event), it is difficult to\nproperly assess the uncertainty related to the other remaining parts. Considering this, we made an\nassumption on the apparent source spectra uncertainty and tried to explain it in Text S2.\nfrequency part, where the spectra should theoretically converge to a point source estimate. We\nhave added a new explanation of our choice of the data error based on the data standard deviation\nat low frequencies, where all the AS should collapse, to Text S2. We have added a new Fig. S12\nshowing the standard deviation (i.e., measure of variability with respect to the average spectrum)\nas a function of frequency, see below:\n2) Figure 1: The figure is generally presented well. A small criticism: I found that the color\ncoding for station azimuth is not clear, in particular w.r.t. azimuthal gaps. Could a rose diagram or\nsimilar be included to see the azimuthal distribution more clearly? Alternatively, one could state\nthe max azimuthal gap (I could not find it in the manuscript).\nWe are sorry for the confusion in Fig. 1 that resulted from using a similar color palette in the\na and b panels. The stations in a) are color coded by their AS variations and not by their azimuth,\nwhile in b) the spectra were color coded by the corresponding station azimuth. To avoid further\nconfusion, we have modified the color palette for Fig. 1a. In addition, we have added rose\nhistograms to Fig. 1a as insets to show the azimuthal coverage and have modified the caption\naccordingly. The maximal azimuthal gaps for directive and nondirective events are 16° and 35°,\nrespectively. We include the comment on this in the text:\n“The azimuthal coverage is good with 16° and 35° gaps for the directive and nondirective\nevents, respectively, yet with a prevalence of stations in the NW-SE direction (see the rose diagram\nin Fig. 1a), following the orientation of the Apennines.”\n3) The rupture is parametrized on a 125-m grid (and there is a second 10-m for FD\ncalculations that is irrelevant for my question). How is this grid size chosen? Does the result\ndepend on this choice? Does the dynamic model require any form of regularization? I have worked\non the effects of discretization choice (with trans-dimensional parametrizations and with\nregularization) and generally found that model selection is quite important for the results. It would\nbe useful to briefly discuss the parametrization choice and whether regularization or some other\nform of model selection is useful here.\nVs=3.5 km/s, it corresponds to fmax=28 Hz) and on our previous experience: Valentová et al.\n(2021) generated a dynamic scenario database and found that the moment rate spectra decayed\nfaster above 2 Hz, which corresponds to the wavelength of their employed model grid (1.2-1.4\nkm). Therefore, for coarser grids, the rupture heterogeneities might not be sufficient to generate\nomega-squared radiation in the frequency range (up to 25 Hz) considered. On the other hand,\nparameters on a finer grid would be unconstrained by the AS data, and thus would require AS\nextended above 25 Hz to infer the possible k^-2 power law signature of fractal properties down to\nsuch lower scales. In addition, a denser grid would significantly increase the number of model\nparameters, making the posterior sampling by the MCMC insufficient. The choice of the grid thus\nserves also as a regularization of the inversion. Another implicit regularization of the kinematic\nrupture parameters is represented by the dynamic modeling itself (elastodynamic equation and\nfriction law). We have added the following text on this issue: “The parameters are defined on a\ngrid of control points 125 m apart, from which they are bilinearly interpolated onto the finer FD\ngrid. The grid size was chosen so that the reliably resolved maximal simulated frequency\ncorresponds to the maximum frequency of the employed data, 25 Hz (compare with, e.g.,\nValentová et al., 2021). For more details, see, e.g., Gallovič et al. (2019a,b).”\n4) Further on the use of GIT: The goal to remove path and station-side effects. This seems a\nbit too good to be true, in particular at frequencies up to 25 Hz. Are the residual errors mentioned\non line 419 in some way accounted for in data error covariance estimation?\nWe discuss the data error estimation in relation to the residual error in Text S2. Moreover,\nwe have added a reference to Text S2 into this part of the manuscript. See also our reply to question\n1.\n5) The methods section (L428 ff) mentions a combination of GIT and BayesISOLA\nsoftware to constrain inputs for the subsequent inversion. Is there any type of uncertainty propagation in\nthis workflow? It was unclear to me if centroid moment tensor results are used in a probabilistic\nsense or as point estimates. A similar mention is on L 481; is Mw constrained according to the\nposterior from BayesIsola or some arbitrary symmetric distribution (there is a mention of\nplus/minus 0.05)?\nIndeed, the centroid as well as hypocenter location uncertainties estimated from BayesIsola\nare not part of the Bayesian dynamic source inversion, i.e., they are not propagated to the\nuncertainty of the source properties. The only exception is M0 estimated from BayesIsola, which\nserves as a prior mean value. We assume a Gaussian PDF in Mw, which corresponds to a log-\nnormal PDF in M0 (see Fig. 6), with a conservative error estimate of 0.05, as the reviewer\n14\nmentioned. As shown in Fig. 6, the dynamic source inversion provided a smaller posterior\nuncertainty in this parameter.\nConsidering the fault plane orientation (strike-dip-rake), we also performed a simple\nsynthetic test, where we assumed an incorrect fault orientation by 7 degrees in Kagan angle, one\nof the BayesIsola results. The effect on the apparent spectra is shown as a spectral bias below and\ndemonstrates that the effect is rather random and more significant at frequencies above the corner\nfrequency. Furthermore, we assume that there are various other effects at these frequencies, such\nas fault plane nonplanarity, contributing to the variability of apparent spectra that are difficult to\ndetermine. Therefore, we neglected the uncertainty in fault plane orientation in the dynamic source\ninversion for now.\nWe thank the reviewer for the suggestion. We admit that neglecting the centroid variability\nand location uncertainty is a simplification that can potentially be improved in further applications.\nWe mention it in the newly added Text S3 about the method limitations (see further).\nMinor comments:\n1) L45: Wording unclear w.r.t. “attain”\nWe have modified the sentence to “Empirical source spectra often exhibit high-frequency\ndecay that deviates from the classical omega-squared model.”\n2) L71: functions (plural)\nCorrected.\n3) S2: 5,600 and not 5.600\nWe have corrected it.\n15\n4) Overall, some editing for English is required. I think I understood the meaning of the text\nmost of the time but both manuscript and supplement will require editing for many small issues as\nin the previous comments.\nWe used Grammarly to check the English grammar in the original version. In the revised\nmanuscript, we have employed AI to improve the English of several sentences in the main text and\nSupplement.\n5) The setup of PT (L515 ff) appears odd to me. Randomly assigning temperatures between\n1 and 100 seems unlikely to result in good acceptance rates for swaps. Were acceptance rates for\ntemperature swaps considered in tuning the algorithm? I my experience this can make an enormous\ndifference and I have never seen an approach similar to the one described here.\nWe consider our choice reasonable because T=100 corresponds to a 10x higher data error.\nSince the data error selection is ambiguous (see Text S2 in the Supplement), especially due to data\ninter-frequency correlation, we think the assumed data error range is rather conservative.\nThe acceptance rates in the chains of our dynamic inversion are mainly affected by the\nparameters of the dynamic models themselves. For example, a perturbation close to the nucleation\nregion generates models that are discarded (not accepted) due to differences in magnitudes. We\nadmit that we tuned the algorithm mainly considering these perturbations and paid less attention\nto the temperature swapping, which should be improved in the future for better effectiveness of\nthe posterior PDF sampling.\n6) L521: The chain thinning applied here should be noted to not positively affect the posterior\nensemble. Rather, it is merely applied to manage storage size. The full posterior (all samples,\nregardless of correlation) always produce better results.\nWe thank the reviewer for this suggestion. We admit that it was introduced to reduce the\npseudoreplication of the samples because the correlation between the accepted dynamic models is\nusually high. We can avoid chain thinning in our future applications. As suggested, we have\ncorrected the sentence to read: “...keeping only every 10th accepted model to manage the storage\nsize of the output models.”\nI hope you find these comments useful and constructive.\nBest wishes, Jan Dettmer (Yukon Geologic Survey / UCalgary)\nI read the paper “Physics-based broadband characterisation of weak earthquakes” by Gallic\net al. with considerable interest. It is a fine piece of work and appears to be of high quality and\ncompleteness. It shows how spatial variations in earthquake amplitude spectra can be used to\ndetermine important details about the rupture properties of even relatively small earthquakes, in\nthis case of magnitudes 4.2 and 4.5. It seems remarkable to see rupture directivity resolved for an\nearthquake as small as magnitude 4.2. Moreover it shows how failing to account for these details\ncan sometimes lead to very significant bias in the conventional methods for estimating important\nparameters like stress drop. I think the paper is an important step forward in estimating earthquake\nsource parameters.\nHowever, while I myself found the paper very interesting I would have to say that I think it\nappeals mainly to a narrow audience of ground motion seismologists. If Nature is fine with\naccepting a paper that may not have broad appeal, then the paper can be published with only minor\nrevision. If instead it is Nature’s goal to publish only papers with broad appeal, then this paper\nmay not be such a good fit. In order to make it appeal to a wider audience I think it would require\nmajor revision, such as perhaps considering a much larger dataset - difficult, given that these two\nevents each took 30 days of run time.\n1) My main substantive comment is that in the section on “GIT derived apparent spectra”,\nsome important assumptions are mentioned, such as the value of the high-frequency attenuation\ncoefficient kappa (=0.015, Line 411), and the use of random perturbations in the strike, dip,\nand rake to mimic random perturbations to crustal velocity structure. These all may be valid\nassumptions, but without further evidence. I’m left wondering whether there could be trade-offs\nbetween some of the source properties and the kappa value. Random crustal velocity perturbations\nwill not only obscure the S-wave radiation pattern but can also produce frequency dependent\napparent attenuation due to scattering. Could such effects be mapped into the source properties if\nthey are not accounted for? In any case, I feel more discussion of these is warranted.\nWe understand the concerns of the reviewer. Accounting for all these effects introduces a\nsystematic change in all spectra over the whole frequency range. These changes can affect the\nestimation of magnitude at low frequencies or small-scale source characteristics at high\nfrequencies. For example, the effect of various values of kappa, corresponding to high frequency\nattenuation of the carefully selected reference sites, is shown in the figure below:\npronounced small-scale rupture characteristics. Conversely, the higher kappa value might have the\nopposite effect.\nDiscussing with other practitioners who use these types of spectral data to derive source\nproperties, we understand that there is no certain evidence on which parameters are the most\nappropriate. Therefore, we tried to account for the site amplification by combining all the effects\nto be consistent with other studies using spectral data. For example, our value of kappa was\nselected to be consistent with rock site conditions in Central Italy (Lanzano et al., 2022; Morasca\net al., 2023). Further research might help to improve the characterization of the residual spectral\namplification and thus improve the source parameter estimates. We expand the discussion on the\npossible trade-offs that can be introduced by these effects in new Text S3.\n2) As with all Bayesian inversions, an important question is: Has it converged? With the\nsampling of 250k models over 30 days run time, I could well believe the solution has converged.\nHowever, some indication of how the authors assessed this, such as examining the character of the\nchain histories of likelihood values, should be discussed.\nWe calculated the Gelman-Rubin criterion (Gelman and Rubin, 1992) for several mean\nparameters (slip, stress drop, or energy estimates). We did not use any particular model parameter\nto track the convergence, because the dynamic parameters may significantly differ between the\nindividual models due to the missing phase information in the data. Histograms for the mean stress\ndrop over individual MPI jobs (pseudochains, in colors) for cold chains (T=1) are shown below,\ndemonstrating the differences between the pseudo-chains:\nwith values <1.1. Only the values of the mean stress drop were a little higher, 1.14 for the directive\nevent. For the non-directive event, we found higher Gelman-Rubin values for stress drop 1.12 and\nrupture area 1.19.\nLine 519: Shouldn’t the total number of models visited by 5k steps x 8 chains x 6 jobs =\n240k, not 250k?\nOur calculation is only approximate. Since the number of steps was larger than 5k, we have\nmodified the sentence as follows: “...each chain performed more than 5,000 steps, i.e., about 250k\nmodels were visited.”\nFigure 1: It can be difficult to make out the stars that denote the epicentres. Could they be\nmade bigger, or a different colour so they’d stand out more? The lines connecting the beachballs\nto the respective epicenters are also very thin.\nAs suggested, we have increased the size of the star (epicenter) and made the line thicker.\nHowever, we have kept the epicenter's black color because we believe it is now sufficiently visible.\nFigures 6 & 7: These are a little confusing because the units of the x-axis at the bottom of\nthe figure are actually given as part of the figure title, at the top of the figure. I recommend moving\nthe units to the bottom of the figure, just below the x-axis, as has been done with other figures.\nWe understand the potential confusion, but including the units in the titles was chosen to\nsave space in the plots, i.e., to keep the individual panels sufficiently large. Note that we also use\nthe same style in other figures (2, 5, 7) and in our previous publications. If acceptable, we would\nrather retain the present style.\nbelow. Note that we have also made some minor improvements to the text."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I appreciate the authors' efforts to revise the manuscript and provide detailed responses to all comments. The study continues to offer a promising approach to dynamic rupture inversion for small earthquakes, and I acknowledge the added value of the new parametric study shown in the rebuttal. However, I am not satisfied with the revision in its current form.\nMy primary concern remains the same as in my initial review. I feel that the authors have not adequately addressed the first major comment of Reviewer #1 regarding the role of abrupt rupture arrest models. While I acknowledge that additional simulations were conducted and shared in the rebuttal letter, the results of this important analysis are not included in the revised manuscript or its Supplementary Information. Presenting such key results only in the rebuttal is inappropriate, as it prevents readers of the published paper from critically evaluating the comparison between abrupt arrest and heterogeneous rupture models.\nThe revised manuscript continues to emphasize the necessity of small-scale heterogeneity to fit the observed high-frequency spectra. This is a central conclusion of the paper and, if correct, would represent a significant advance. However, without properly documenting the performance of sudden arrest models—beyond the gradual-arrest smooth models shown in the SI—the claim that heterogeneity is required remains insufficiently supported. The authors themselves acknowledge that barrier models with fast rupture can produce omega-squared falloff and fit the observed spectra reasonably well under certain conditions. Yet this important nuance is hidden in the rebuttal and omitted from the actual manuscript.\nIf the barrier models performed significantly worse than the heterogeneous models, then the comparison should be documented and discussed explicitly in the paper. If the differences are more subtle—perhaps depending on trade-offs with rupture speed or stress drop—then that too should be part of the discussion. In either case, the reader must be given access to these results in order to independently evaluate the authors' interpretation.\nI also note that the revised version makes several definitive statements implying that heterogeneity is required (e.g., Abstract, Lines 21–22: \"In our models, the heterogeneity down to the smallest scales (~100 m) is required to fit the observed high-frequency radiation.”). In light of the alternative modeling presented in the rebuttal, such statements are too strong without presenting the full scope of tested models.\nTherefore, I encourage the authors to revise the manuscript by:\n(i) Incorporating the abrupt arrest (barrier) model results into the manuscript or SI in a clearly labeled section, including key figures and discussion of misfits, parameter trade-offs, and spectral comparisons.\n(ii) Clarifying the framing of their conclusions regarding the necessity of heterogeneity in light of the alternative models.\n(iii) Moderating claims in the Abstract, Results, and Conclusions to reflect the fact that other models (with abrupt arrest) can also explain aspects of the high-frequency radiation, albeit with limitations or different trade-offs.\nThese additions would substantially improve the rigor, transparency, and impact of the paper. Until these revisions are made, I would not consider the current manuscript suitable for publication in Nature Communications.\nMinor comment:\nLine 62: The authors mention that “Yoshida (2019) demonstrated that asymmetrical rupture propagation may result in biased stress drop estimates if the directivity effect on the source spectra is neglected.” It could be helpful to also acknowledge that methods such as the second moment approach (e.g., McGuire et al., 2017; McGuire and Kaneko, 2018) offer an alternative pathway to improved stress drop estimates by better resolving source geometry and rupture propagation direction.", "reviewer_2": "I would like to thank the authors for an exceptionally thorough respoinse to my comments. I have no further concerns with this manuscript and recommend it be considered for publication. It's excellent work in my opinion."}, "author_response": "We agree that some of our statements regarding the rupture heterogeneity were confusing.\nIndeed, in our study, the heterogeneous distribution of parameters is an assumption, and what we\ninfer is the power-law properties of the parameters. We are unable to disprove that the other types\nof models are incorrect, but we consider our parameterization more realistic when fitting real data\nthan the smooth ruptures, which are undoubtedly important from the theoretical point of view. We\nalso agree that performing an appropriate inversion with the smooth barrier models would be\ninteresting, but it cannot be fully incorporated in our manuscript as it requires a devoted in-depth\nstudy.\nheterogeneity is required (e.g., Abstract, Lines 21–22: \"In our models, the heterogeneity down to the\nsmallest scales (~100 m) is required to fit the observed high-frequency radiation.”). In light of\nthe alternative modeling presented in the rebuttal, such statements are too strong without\npresenting the full scope of tested models.\nWe have modified the Abstract to clearly distinguish between our parameterization\n(heterogeneous) and the properties required to fit the data (power-law spectra).\nTherefore, I encourage the authors to revise the manuscript by:\n(i) Incorporating the abrupt arrest (barrier) model results into the manuscript or SI in a clearly\nlabeled section, including key figures and discussion of misfits, parameter trade-offs, and spectral\ncomparisons.\nWe were very cautious to present the parametric study of the barrier models within the paper\n(and keeping it public only as a rebuttal), as it provides a very limited number of manually set\nsimulations compared to the McMC inversion that visited 250k models. Therefore, it cannot be\nconfronted with the results of the Bayesian inversion in the main text. In addition, the spectral\nproperties of such models have already been properly studied, e.g., Kaneko and Shearer (2015),\nwhich we have referenced. Based on the suggestion, we have added the tests presented in the\nprevious rebuttal letter to the SI as a new Text S4. We refer to this text when clarifying our\nparameterization in the Materials and Methods.\n(ii) Clarifying the framing of their conclusions regarding the necessity of heterogeneity in\nlight of the alternative models.\nWe have modified the Abstract and other statements to clearly express that our inversion\ninfers the power law for the assumed heterogeneity.\n(iii) Moderating claims in the Abstract, Results, and Conclusions to reflect the fact that other\nmodels (with abrupt arrest) can also explain aspects of the high-frequency radiation, albeit with\nlimitations or different trade-offs.\nIn all relevant parts, we modify the text to clearly distinguish our assumption (heterogeneity)\nfrom the inferred properties (power-law). We also mention the smooth barrier models as an\nalternative, but we prefer heterogeneity as a more realistic description of the underlying complex\nfaulting process, which we support by adding references.\nUntil these revisions are made, I would not consider the current manuscript suitable for publication\nin Nature Communications.\nMinor comment:\nLine 62: The authors mention that “Yoshida (2019) demonstrated that asymmetrical rupture\npropagation may result in biased stress drop estimates if the directivity effect on the source spectra\nis neglected.” It could be helpful to also acknowledge that methods such as the second moment\napproach (e.g., McGuire et al., 2017; McGuire and Kaneko, 2018) offer an alternative pathway to\nimproved stress drop estimates by better resolving source geometry and rupture propagation\ndirection.\nWe agree that another way of improving the stress drop estimates is by constraining the fault\nsize instead of the corner frequency, following also the simulation analysis in Gallovič and\nValentová (2020). However, such an estimate by, e.g., the second moment approach will also be\nsubject to uncertainty, which will be translated into uncertainty of the stress drop. We have added\na note together with new references to the Discussion part, where we discuss the potential bias in\nstress drop estimation in detail.\nReviewer #2 (Remarks to the Author):\nI would like to thank the authors for an exceptionally thorough response to my comments. I\nhave no further concerns with this manuscript and recommend it be considered for publication. It's\nexcellent work in my opinion.\n3"}, {"round": 3, "reviewer_comments": {"reviewer_1": "I appreciate the careful and constructive responses to the previous review comments. The authors have appropriately revised the manuscript text to ensure consistency with the presented results and have added a detailed description of the abrupt-rupture-arrest model simulations in the Supplementary Information (Text S4). Although they did not conduct a full inversion or quantitative comparison for abrupt-rupture-arrest models as I originally requested, the revised version now clearly articulates the underlying assumptions, limitations, and rationale of their analysis. (I still consider such an analysis to be an important next step for future work.) In any case, the central findings are now presented in a balanced and transparent manner. Given these revisions, I consider the manuscript suitable for publication in Nature Communications in its current form."}, "author_response": ""}]

Decision letter