Analysis of the Szarvas, Hungary, Earthquake Swarm of 19 August 2023: A Case Study in a Low‐Seismicity Sedimentary Basin Open Access

Between 19 August and 26 September 2023, an earthquake sequence took place near the town of Szarvas, southeastern Hungary (Fig. 1). During this period, 105 earthquakes with local magnitude $ML$ between 0.9 and 4.1 were recorded in the region (Süle et al., 2024). The sequence was initiated by two events of $ML$ 4 and 4.1, occurring just 3 min apart, and causing minor damage in the epicentral area.

The temporal and magnitude distribution of the earthquakes did not follow a typical pattern (Fig. 2). We cannot distinguish a single mainshock followed by several aftershocks with decreasing magnitude. Instead, there were seven $ML≥3.5$ earthquakes in the first 24 hr, including three $ML≥4.0$ events. Then, on the fourth day, another $ML$ 3.7 earthquake could be observed. Overall, 89 earthquakes occurred in the first 24 hr and 95 in the first 48 hr, representing 85% and 90% of the total number of recorded events, respectively. After 25 August, only two events occurred. Consequently, rather than a typical series of earthquakes, we classify these events as an earthquake swarm that persisted for several days. Earthquake swarms of this type are unusual in Hungary. Therefore, it has attracted a lot of attention not only from researchers but also from the public and the media.

The central region of the Pannonian basin, encompassing Hungary, exhibits moderate seismic activity (Tóth et al., 2008). The distribution of epicenters spans the entire country; nonetheless, there are localized areas with increased seismicity for which destructive earthquakes have been documented in the last centuries. As a consequence of the earlier extension, the Pannonian basin features a thin lithosphere, a shallow Moho discontinuity, and a high surface heat flow (Lenkey et al., 2002; Grad and Tiira, 2009; Kalmár et al., 2021; Timkó et al., 2024). Therefore, earthquake activity is limited to the upper crust; most seismic events occur at depths ranging from 5 to 15 km (Zsíros, 1996; Czecze et al., 2023).

The epicentral area of the 2023 Szarvas earthquake swarm does not fall within the regions of highest seismic activity in the country. A total of 27 events are known to have occurred within 30 km of the current epicenters, but only two of them had a magnitude greater than 3. Both occurred in 2008, on the same day at Kondoros, approximately 10–15 km from the location of the present earthquakes (Fig. 3). Their local magnitude were 3.7 and 3.5. However, at the same time, we should mention the Békés earthquake of 22 June 1978, which occurred about 40 km from the current epicenters (Zsíros, 2000). It was the most powerful seismic event in the region with a magnitude of $ML$ 4.8 and was followed by several felt aftershocks. Apart from the relatively distant Békés earthquake, most of the events in the area were weak. In contrast, during the 2023 Szarvas earthquake swarm, three earthquakes of approximately the same magnitude 4 occurred within the first 24 hr, which has never been experienced before in the region or even in the country.

In 2022, the Kövesligethy Radó Seismological Observatory (KRSO) of the Institute of Earth Physics and Space Science (EPSS) joined the international AdriaArray project (Kolínský and Meier, 2023). The AdriaArray project addresses fundamental questions related to the geodynamics and the deformation and stress field of the Adriatic plate in southeastern Europe. It will lead to a significant improvement in our understanding of the geodynamic causes of plate deformation and associated geohazards. In the framework of the AdriaArray project, 32 temporary broadband stations have been deployed in Hungary. As a result, together with the 15 permanent broadband stations, there are now 47 stations operating in the country. Such a dense network of stations has never operated in Hungary before. Therefore, the 2023 Szarvas earthquake swarm is the first one in the region that was recorded by a significant number of high‐quality broadband seismic stations, offering an excellent opportunity to investigate earthquake sources in an area where fault parameters and small‐scale tectonic structures are poorly understood.

Taking advantage of the dense temporary seismic network operating in the country, in this article, we thoroughly analyze the source parameters of the 2023 Szarvas earthquake swarm. After a brief review of the tectonic setting, macroseismic observations, and data sources used in this study, we recalculate the hypocenter locations of the events applying the probabilistic BayesLoc multiple‐event location algorithm (Myers et al., 2007, 2009). By focusing on a limited spatial and temporal domain, we demonstrate the effectiveness of the method in a dense, local network setting. Although previous research has illustrated the advantages of BayesLoc in large‐scale applications (Myers et al., 2015; Hicks and Myers, 2017; Hofstetter and Shamir, 2021; Bondár et al., 2023; Czecze, 2024), our study provides valuable insights into how it performs in a small‐scale local seismic sequence. We also estimate the full moment tensor (MT) of the eight largest events using the probabilistic nonlinear Monte Carlo moment tensor (MCMT) waveform inversion method (Wéber, 2006, 2009, 2016a). The probabilistic approach allows for estimating the uncertainty of the focal mechanisms, which is a crucial prerequisite for reliable geodynamic interpretation.

This work demonstrates that the state‐of‐the‐art probabilistic BayesLoc and MCMT methods are powerful tools that can be successfully applied for detailed analysis of local event clusters, encouraging their further use in other small‐scale seismicity studies.

Hungary is located in southeast Central Europe, in the central part of the Pannonian basin. The Pannonian basin is a sediment‐filled back‐arc basin for which formation started in the Early Miocene as a result of the convergence between the European and African plates (Horváth and Cloetingh, 1996; Tari et al., 1999; Matenco and Radivojević, 2012; Horváth et al., 2015; Balázs et al., 2016). The subduction of the European foreland and the subsequent rollback of the subducted slab induced lithospheric extension and caused asthenospheric updoming during the Middle Miocene, leading to the formation of distinct basins due to the extension (Schmid et al., 2008; Handy et al., 2015; Matenco and Radivojević, 2012; Handy et al., 2021). By the Late Miocene, the cooling of the asthenospheric dome resulted in the subsidence of the entire basin (Horváth and Cloetingh, 1996; Horváth et al., 2006). In the latest Pliocene and Quaternary periods, the movement of the Adriatic microplate halted the subsidence and initiated the inversion of the Pannonian basin (Horváth and Cloetingh, 1996; Fodor et al., 2005).

The recent tectonic activity of the region is controlled primarily by continuous northeast compression and counter‐clockwise rotation of the Adriatic microplate relative to Europe around a pole in the Western Alps (Bada et al., 1999; Bada, Horváth, et al., 2007). Horizontal crustal velocities determined by Global Navigation Satellite Systems measurements relative to the stable Eurasian plate show an up to 3 mm/yr, north‐northeast‐directed motion of the Dinarides and an up to 1.5 mm/yr, south‐southwest‐directed motion of the Eastern and Southern Carpathians (Fig. 4). The Pannonian basin and the Western Carpathians are characterized by very low velocities of 0.1–0.5 mm/yr, with slightly higher values in the southwest part of the Pannonian basin (Porkoláb et al., 2023). As a result, compressional stress with complex stress pattern dominates in the basin (Bada, Horváth, et al., 2007; Békési et al., 2023; Porkoláb et al., 2024), and it is considered responsible for horizontal shortening and the reactivation of older faults, which is manifested in the intraplate seismicity of the region. Because of the compressional stress, strike‐slip to compressive faulting can be observed inside the basin (Békési et al., 2023; Koroknai et al., 2023; Porkoláb et al., 2024).

The 2023 Szarvas earthquake swarm occurred in the southeast part of Hungary. In this area, the Conrad surface is situated at a depth of 16 km, whereas the Moho discontinuity is located at 22 km depth, based on P‐to‐S receiver function analysis conducted by Kalmár et al. (2021). This means that the thinnest crust within the Pannonian basin is located here, with the lower crust measuring only 6 km in thickness. The earthquakes occurred at the northwest edge of the Békés basin (Fig. 1), which is one of the deepest (more than 6 km deep) and largest Neogene sub‐basins inside the Pannonian basin. The thickness of the sediments in the epicentral area is about 3.4–3.7 km (Balázs et al., 2018; Kalmár et al., 2021). The surface of the pre‐Neogene basement is fragmented; its geological structure is complex.

As many as 13 events of the 2023 Szarvas earthquake swarm were felt by the population, 10 of them in the first 24 hr. There was minor damage to buildings (tiles and plaster falling, chimney damage, and wall cracks) in the epicentral area, but no people were injured.

More than 800 macroseismic questionnaires were collected via the online form of the KRSO. In the case of the first two $ML≥4.0$ events, the evaluation of the questionnaires was difficult because there was only a 3‐minute time difference between them. It was not always clear which earthquake’s effects were reported in a given questionnaire. Nevertheless, the maximum epicentral intensity for these two events, and the whole swarm, could be estimated as V–VI on the European Macroseismic Scale 1998 (EMS‐98) (Grünthal et al., 1998; Süle et al., 2024). The maximum intensities of the felt earthquakes are listed in Table 1. More details on the collected questionnaires and the intensity maps for the 13 felt events can be found in Süle et al. (2024).

We have also created ShakeMaps (Worden et al., 2020; Wald et al., 2022) for the intensity distribution of the three largest $ML≥4.0$ earthquakes. For their preparation, we used both instrumental data and the intensities derived from the submitted questionnaires. The obtained isoseismals are symmetric around the epicenters, that is, no anomalous intensities due to directionality or local effects can be observed.

Waveform data used in the present study were mainly recorded by four broadband seismic networks. Besides the permanent stations of the Hungarian National Seismological Network (HNSN, network code HU) and the Romanian Seismic Network (network code RO), stations of two temporary networks (network codes Z6 and Y8) belonging to the AdriaArray initiative (Kolínský and Meier, 2023) played an important role in detecting the investigated earthquake swarm. The achieved good azimuthal coverage and short epicentral distances would not have been possible without the temporary AdriaArray stations (Fig. 1).

At the time of the 2023 Szarvas earthquake swarm, there were 15 permanent and 32 temporary broadband seismic stations in Hungary, forming a network with a density of about 40–50 km. Moreover, on 21 August 2023, the KRSO installed two additional short‐period temporary stations (STMP1 and STMP2, Fig. 1) in the close vicinity of the epicentral region to better record even the smallest events. Both were equipped with Le3D/5s seismometers. Permanent and temporary broadband stations in Romania also played a crucial role in achieving the exceptional azimuthal coverage. Additional data sources are detailed in the Data and Resources. All broadband stations used in this study were equipped with three‐component seismometers with a long‐period corner of 120 s.

Reliable determination of earthquake locations and magnitudes is crucial for understanding tectonic processes, identifying active fault zones, and conducting seismic hazard assessments. At the KRSO, hypocenter locations are routinely calculated by the iLoc single‐event location algorithm (Bondár and Storchak, 2011). Bondár et al. (2018) demonstrated that the iLoc algorithm produces high‐quality locations in the Pannonian basin when using travel‐time predictions from the global 3D upper‐mantle regional seismic travel time (RSTT) model (Myers et al., 2010). Therefore, earthquake localization using the iLoc algorithm constitutes a fundamental part of our observatory’s routine processing (Czecze et al., 2023).

Although iLoc provides reliable single‐event locations, multiple‐event location techniques allow for the simultaneous adjustment of multiple seismic events within a cluster, leveraging their relative positions and common travel paths. To further improve the location accuracy of the 2023 Szarvas earthquake swarm, we used the BayesLoc multiple‐event location algorithm in this study (Myers et al., 2007, 2009). BayesLoc employs a Bayesian approach to jointly locate multiple events by statistically accounting for uncertainties in both travel times and event locations. This method optimally uses the spatial and temporal correlations between events in a cluster, significantly improving the relative accuracy of event locations compared to single‐event methods. BayesLoc reduces location uncertainties, particularly in regions with complex geological structures or sparse station coverage (Wéber, 2016b; Schardong et al., 2021; Visnovitz et al., 2021; Bondár et al., 2023; Czecze, 2024). The method is also able to identify and remove outliers, which can otherwise significantly bias the results of single‐event algorithms.

BayesLoc has been extensively applied in regional and global‐scale seismic studies, as demonstrated by previous works on the Pannonian basin (Czecze, 2024), the Arctic plate boundary (Hicks and Myers, 2017), and the Caucasus region (Bondár et al., 2023). Other applications include global seismicity studies (Myers et al., 2015) and revisions of national earthquake catalogs (Hofstetter and Shamir, 2021). All these studies focus on large‐scale datasets of regional or teleseismic events. In contrast, in this study, BayesLoc is applied to a small, localized earthquake swarm, demonstrating its usability in high‐resolution local seismicity studies.

In this section, we first provide a brief description of the BayesLoc method and how it was configured in this study. Then, we present and analyze the obtained hypocentral distribution of the investigated earthquake swarm.

BayesLoc (Myers et al., 2007, 2009) applies a statistical model that decomposes the multiple‐event problem into three model components: the travel‐time model, the arrival‐time model, and the prior model. Each model component is represented as a conditional probability and integrated into a joint posterior probability distribution (PPD) through the application of Bayes’ theorem. BayesLoc employs Markov chain Monte Carlo (MCMC) sampling to generate realizations of hypocentral parameters from the joint PPD.

The travel‐time model allows for corrections to model‐based predictions to mitigate location biases and improve the accuracy of the applied earth model by achieving a closer fit to the observed data. The corrections consist of multiple components: adjustments to the slope of the travel‐time curve and static corrections for each station, event, phase, event‐phase pairs, and station‐phase pairs. The prior distribution for these correction parameters is Gaussian with zero mean. In this study, we applied phase‐specific and event‐specific static corrections, as well as a scaling factor that adjusts the slope of the travel‐time curve.

Errors in the measured arrival times are addressed by the arrival‐time model once travel‐time corrections have been implemented. The a priori errors in arrival time are assumed to follow a Gaussian distribution, characterized by zero mean and variance derived from the pick precision (1/variance). Pick precision is decomposed into three factors: phase, event, and station precisions. This decomposition accounts for common data quality issues, such as the impact of lower‐quality stations or weak signal amplitudes in the case of low‐magnitude events. Precision factors act to assign greater weight to low‐variance data and lesser weight to high‐variance data. The prior distribution for precision follows a Gamma distribution. The parameters of the Gamma distribution are treated as random variables and are iteratively adjusted during the BayesLoc process. In this study, we included all three precision factors in the calculations.

Probabilistic seismic phase labels are also integrated into BayesLoc to handle potential phase misidentifications (Myers et al., 2009). Finally, the origin prior model allows for the specification of constraints on parameters such as the epicenter, depth, or origin time for known events. In this study, we used the high‐quality iLoc hypocenters as initial locations for the MCMC sampling but applied noninformative priors to the parameters of event origins (epicenters, depths, and origin times).

We used vertical‐component seismograms to manually pick P‐wave arrivals and horizontal‐component waveforms to identify S‐wave arrivals where possible. To ensure that only high‐quality events were subjected to the BayesLoc analysis, we performed a preliminary event selection. Those events that were detected by fewer than five stations and had less than seven phase readings were excluded from the analysis. After applying these selection criteria, we found 83 $ML≥0.9$ earthquakes with a sufficient number of high‐quality arrival‐time observations. We collected 2093 arrivals in total, consisting of 655 Pg, 567 Pn, 633 Lg, and 238 Sn phases.

The geometry of a station network recording a seismic event is best described by the secondary azimuthal gap, which represents the maximum azimuthal gap when one station is removed from the recording stations. The preliminary event selection procedure described earlier ensured that the largest secondary azimuthal gap remained below 236° for the whole dataset, indicating a well‐functioning station geometry. More specifically, for the largest events (⁠$ML≥3.5$⁠) the secondary azimuthal gap does not exceed 65° and is below 150° for 52 of the relocated 83 earthquakes. Only for the smallest earthquakes does it increase above 150°, probably due to the limited detection threshold of the station network.

To map the joint PPD of the 2023 Szarvas earthquake swarm, we applied eight independent MCMC chains with chain mixing to ensure robust convergence. For each chain, the last 5000 of 10,000 MCMC samples were used to produce the results presented here, with the initial 50% of iterations serving as a burn‐in period to minimize the effect of the poorly constrained initial locations. In fact, numerical tests demonstrated that in our case the burn‐in period is no longer than 2500–3000 iterations.

The seismic stations that provide data for the Szarvas earthquake swarm are situated in very different geological settings with very different local velocity structures (Fig. 1). Therefore, it is not possible to define a single 1D velocity model that is appropriate for each seismic station. Instead, we used the global ak135 earth model (Kennett et al., 1995) with ellipticity corrections to predict travel times, and allowed BayesLoc to adjust the model according to the observed data through its statistical travel‐time model.

BayesLoc provides a probabilistic assessment of the accuracy of the input arrival time picks by evaluating how well they fit the model. If a pick is identified as inconsistent with the model, BayesLoc either reassigns the phase label or flags it as an outlier (Myers et al., 2009). According to our results, over 96% of the input phase labels are in agreement with the posterior ones. A total of 70 arrivals, representing 3.3% of the dataset, have a most probable posterior phase name that differs from the input phase label. More specifically, 6.7% of the input Pn arrivals were relabeled as Pg, 0.2% of the Lg arrivals as Sn, and 13.8% of the Sn arrivals as Lg. No Pg arrivals were relabeled. Only one input Pn phase was marked as an outlier. Overall, it can be concluded that the manually picked input phase labels are of high quality.

The maximum‐likelihood hypocentral parameters obtained for the 2023 Szarvas earthquake swarm are summarized in Table 1 and Figures 5 and 6. The earthquakes are concentrated near Csabacsűd in a narrow region with a radius of a few kilometers. This result is in good agreement with the observations reported in the collected macroseismic questionnaires. The standard deviation of the epicentral coordinates is mostly below 2 km, and even below 1 km for the larger events (Table 1). Figure 5 demonstrates that the posterior BayesLoc epicenters are more spatially concentrated than the prior locations calculated by routine observatory data processing (iLoc).

The depth distribution of the swarm does not show any particular structure (Fig. 6). The earthquakes essentially occurred at depths between 1 and 10 km with a standard deviation of mostly around 2–4 km. According to the histogram in Figure 6c, more than one‐third of the events took place at shallow depths in the topmost 3–4‐km thick depth range.

Figure 7 presents 2D marginal PPDs illustrating the joint horizontal and vertical distributions of the hypocenters. The joint PPD reveals where the events of the swarm are most likely to have occurred. Again, the epicentral PPD shows that the earthquakes are concentrated within a limited geographical area, indicating a localized seismic activity. The two vertical PPDs also support our earlier finding that the earthquakes originated in the top 10 km of the upper crust. The events took place with high probability near the surface, in the uppermost 3–4 km depth range, as we have already seen in Figure 6c. On the other hand, the local maximum at depths around 8 km in the 2D PPDs is not so pronounced in the depth histogram. Without additional independent data, it is now difficult to decide whether or not this deep local maximum reflects a true small‐scale fault structure.

The distribution of depth uncertainties (standard deviations) for the relocated events is depicted in Figure 8a. Most events have depth uncertainties concentrated between 2 and 6 km, with a peak around 3 km, indicating relatively well‐constrained depths for the majority of events. Standard deviations extend up to around 10 km, but only a few events exhibit such large uncertainties. In contrast, the standard deviation for the iLoc depths typically reaches values of as large as 6–10 km.

Figure 8b compares the size (area) of the formal error ellipses for the prior and posterior hypocenters, calculated at the 90% confidence level. The error ellipses reflect the horizontal uncertainty in the hypocenter locations. The prior hypocenters have a broader range of ellipse sizes, with a significant portion of events exhibiting high uncertainty, including some reaching an area of $400km2$⁠. In contrast, the posterior hypocenters show a significant reduction in ellipse size. Most relocated events have error ellipses smaller than $50km2$⁠. To estimate a more realistic location error instead of a formal one, we would need the presence of ground‐truth events (Bondár and McLaughlin, 2009). However, no anthropogenic events with known locations are available within a practically relevant distance to the swarm.

Summing up, the BayesLoc multiple‐event algorithm resulted in substantial improvements in both depth and horizontal uncertainties after the relocation process, as illustrated in Figure 8. Depth uncertainties are concentrated at lower values, and the size of the error ellipses has significantly decreased for most events relative to the initial iLoc locations.

We also investigated the relation between hypocentral depth and local magnitude. Local magnitudes were calculated using the iLoc software (Bondár and Storchak, 2011), which follows the SeisComp3 (SC3) methodology (Helmholtz‐Centre Potsdam‐GFZ and GEMPA GmbH, 2008). Figure 9 illustrates the magnitude distribution and hypocentral depth, highlighting the depth–magnitude relationship. According to the magnitude histogram (Fig. 9a), about 70% of the events fall within the $1.5≤ML≤2.5$ range, demonstrating that the 2023 Szarvas earthquake swarm is dominated by minor earthquakes. These minor events predominantly occurred at shallow depths (above 6 km). In contrast, earthquakes with $ML≥3.0$ consistently originated at depths greater than 6 km (Fig. 9b).

In this section, we determine the source mechanisms (full MTs) for eight $ML≥3.5$ events from the 2023 Szarvas earthquake swarm. We begin by briefly describing the probabilistic nonlinear waveform inversion methodology applied and then examine the selected events in detail, discussing the obtained MT solutions.

In this study, we applied the probabilistic nonlinear MCMT inversion method (Wéber, 2006, 2009) to retrieve the seismic MT for the selected earthquakes. It has already been demonstrated that this inversion technique can produce reliable focal mechanism solutions for both local and near‐regional events (Wéber and Süle, 2014; Wéber, 2016a,b; Wéber et al., 2020; Békési et al., 2023). In the next few paragraphs, we provide a brief description of the method.

We describe a general seismic point source by a symmetric MT. For a given hypocenter and source time function (STF), the MT is estimated by deconvolving the Green’s functions (GFs) from the observed seismograms. We utilize a Monte Carlo simulation technique (e.g., Rubinstein and Kroese, 2008) to assess the overall uncertainties of the derived MT solutions. The Monte Carlo simulation method enables the determination of how variability in the input data propagates to the uncertainty of the output. Applying this principle to our waveform inversion problem, we generate many new realizations of input datasets by randomly generating new hypocenters and waveforms according to their respective distributions (Wéber, 2016a; Wéber et al., 2020). Then each generated input dataset is inverted for MT (output). The distribution of the obtained set of MT solutions approximates well the PPD of the MT. The final estimate for the best MT is defined by the maximum‐likelihood point.

Given the samples from the PPD, several statistical properties of the MT can be determined. First of all, we calculate the principal axes for each MT solution obtained by the aforementioned Monte Carlo simulation. We adopt the convention of Sipkin (1993) that the P and T axes always point upward and the principal axes form a right‐handed coordinate system. Then, we construct the 2D histograms of the principal axes on the focal sphere and determine the confidence zones for the desired confidence levels. The confidence contours of the P and T principal axes are then plotted on top of the focal mechanism plot of the maximum‐likelihood mechanism (Wéber et al., 2020). We also decompose each MT solution into a double‐couple (DC), a compensated linear vector dipole (CLVD), and an isotropic (ISO) component (Jost and Herrmann, 1989). To assess the relative amounts of these components in an MT, we calculate their percentages as well.

For full details of the MCMT inversion method, see Wéber (2006, 2016a, b).

Stations within about 80–100 km (or even more) from the epicentral region are located on top of a several kilometer thick sedimentary basin, the Great Hungarian Plain (Fig. 1). In such a geological environment, earthquake signals may be contaminated by significant amounts of microseismic and anthropogenic noise, especially in daytime. Therefore, it is essential to perform rigorous noise analysis when selecting the events to be processed and the waveforms to be used for waveform inversion. Therefore, before inversion, we first carefully analyzed the frequency‐dependent signal‐to‐noise ratio (SNR) for each component of each observing seismic station to select the waveforms suitable for inversion, together with their optimal frequency band. Because the high‐frequency content of the waveforms may be strongly affected by small‐scale heterogeneities of the medium not included in our simple 1D velocity model, it is necessary to use as low frequencies as possible with good SNR to reliably estimate earthquake focal mechanisms by waveform inversion.

After selecting the waveforms of sufficient quality, we applied a Butterworth band‐pass filter to the displacement seismograms according to the previously selected frequency band. The displacement GFs were subjected to the same filtering process. Because the inversion frequency band was set below the events’ corner frequency, we assumed the STF to be an impulse. In this study, we used a simple triangular STF before band‐pass filtering. The length of the processed time window was chosen according to the epicentral distance and the applied frequency band. The synthetic GFs were windowed in the same way as the observed seismograms. Because a simple 1D earth model is not able to accurately predict the arrival times of the different seismic phases, we applied a time shift between the observed waveforms and the synthetics to obtain the maximum correlation between them.

After carefully evaluating the frequency‐dependent SNR of the recorded seismograms, we finally selected eight earthquakes with $ML≥3.5$ for detailed source analysis. The most usable frequency band ranges between 0.05 and 0.1 Hz with SNR greater than 3. At these low frequencies, waveforms are unlikely to be significantly affected by small‐scale heterogeneities of the medium not incorporated in our 1D earth model. However, it should be noted that for many stations it was not possible to use all three components in the inversion due to the low SNR at almost all frequencies.

In this study, we employed a frequency–wavenumber integration method (Wang and Herrmann, 1980; Herrmann and Wang, 1985; Herrmann, 2013) for calculating the synthetic GFs. Data from seismic stations with epicentral distances greater than about 150 km were not included in the waveform inversion procedure, so we used the local, horizontally layered earth model of Gráczer and Wéber (2012) to calculate the synthetics. To estimate the uncertainty of the resulting MT, we performed 10,000 Monte Carlo simulations and thus generated 10,000 MT solutions from the posterior distribution. The final MT solution is given by the maximum‐likelihood point. Because we believe that the previously estimated BayesLoc epicenters are accurate enough for waveform inversion, we searched only for hypocentral depths in the MCMT inversion.

The resulting focal mechanisms are listed in Table 2 and plotted in Figure 10. The figure shows the focal mechanism plot of the deviatoric part of the full MT solutions together with the 50%, 68%, 90%, and 95% confidence contours for the P and T principal axes. For verification, P‐wave polarities are also displayed on top of the focal mechanism plots. However, it should be noted that polarities were not included in the inversion. The scalar seismic moments and principal axes of the maximum‐likelihood mechanisms are listed in Table 3 together with their 95% confidence regions.

According to Figure 10 and Tables 2 and 3, the resulting MTs all represent strike‐slip events with very similar orientations. The azimuth of the sub‐horizontal P‐axes varies within a narrow range between 70° and 77°. The mechanisms are generally in good agreement with the available first‐arrival P‐wave polarities. Both the P and T principal axes are strongly clustered around well‐defined directions with small uncertainty. Because of the fairly large number (36–58) of the inverted waveforms (Table 2) and the good azimuthal coverage, the derived MT solutions are reliable and of high enough quality to facilitate further tectonic investigations.

Table 2 demonstrates that the DC component of the obtained MTs ranges from 82% to 95%, whereas the CLVD and ISO components do not exceed 10%. The negligible amount of the ISO component suggests that the studied events are of tectonic nature. The centroid depths achieved by the waveform inversion fall within the 95% confidence intervals of the hypocenter depths (compare with Tables 1 and 2).

The quality of waveform fitting is illustrated in Figure 11, which compares the observed seismograms for event 37 (⁠$ML$ 3.8) with synthetic waveforms derived from the best (maximum likelihood) MT. Two quantities are given for each seismogram: the normalized correlation coefficient c and the variance reduction vr. Because of the good SNR in the inversion frequency band and the azimuthally well‐distributed recording stations, the synthetics and the data usually correlate satisfactorily. The variance reduction values, on the other hand, highlight that the observed waveform amplitudes are not always well reproduced; nonetheless, the waveform matching achieved is acceptable.

In the online catalog of the German Research Centre for Geosciences (GFZ), regional MT solutions have been published for the first two $ML≥4.0$ events of the investigated earthquake swarm. These regional MTs show very good agreement with the source mechanisms retrieved in this study (Fig. 12, Table 2).

To quantitatively compare our solutions to those published by the GFZ, we follow Pasyanos et al. (1996) and use their MT difference parameter $μ$⁠. It is defined as the root mean square of the differences of the MT elements normalized by their respective scalar moment. According to the authors, $μ<0.25$ corresponds to nearly identical focal mechanisms, whereas divergence becomes evident for $0.25<μ<0.50$⁠, and $μ>0.50$ indicates significantly differing mechanisms.

As seen in Table 2, the $μ$ values for events 1 and 2 are 0.19 and 0.23, respectively. Both values are below 0.25 confirming that the agency solutions are practically identical to ours. Indeed, the differences between the P axes are negligible, only the T axes differ slightly (Fig. 12). We consider our solutions to be better because they contain much smaller non‐DC components and we have used more local stations in the inversion.

In this article, we have thoroughly analyzed the source parameters of the 2023 Szarvas earthquake swarm, an exceptional seismic event cluster beneath a low‐seismicity sedimentary basin in Hungary. The unusual temporal and magnitude distribution of the swarm poses a challenge to the conventional interpretation of earthquake sequences in the region.

Utilizing the BayesLoc multiple‐event location algorithm, we have recalculated the hypocentral locations of 83 earthquakes with $ML$ between 0.9 and 4.1. The application of BayesLoc significantly reduced the location uncertainties in both horizontal and vertical directions, compared to the routinely used iLoc algorithm. The events are concentrated in the close vicinity of Csabacsűd, within a small area with a radius of a few kilometers. The epicentral region is situated on top of the Dévaványa–Gyomaendrőd (DG) fault zone (Fig. 5). The 2023 Szarvas earthquake swarm provides the first clear evidence that this fault zone is active.

The earthquakes occurred in a confined volume of a few kilometers horizontally and about 10 km vertically. At the same time, the 90% confidence zone of the obtained hypocenters falls within the range of ± (1−5) km horizontally and around ± 5 km vertically (Fig. 8, Table 1). This resolution does not allow the resulting hypocentral distribution to reveal any small‐scale tectonic features within such a small source volume, which may explain why Figures 5–7 do not display any well‐defined patterns.

In the epicentral area of the 2023 Szarvas earthquake swarm, the basement is covered by about 3.4–3.7 km thick sediments (Balázs et al., 2018; Kalmár et al., 2021). Figures 6c and 9b show that more than one‐third of the events, all with $ML<3.0$ magnitude, occurred above the basement in the sedimentary layers. The larger $ML≥3.0$ events, on the other hand, originated at depths between 6 and 10 km in the upper crust, well below the sediments.

By employing a probabilistic waveform inversion methodology, we have estimated the full MTs for eight earthquakes with $ML≥3.5$⁠. The non‐DC components within the retrieved MT solutions are statistically insignificant. The negligible ISO component indicates that the events are of tectonic nature. The centroid depths obtained by the waveform inversion fall within the 95% confidence intervals of the hypocenter depths.

The MT solutions displayed in Figure 10 reveal that all of the investigated earthquakes have strike‐slip mechanisms with either right‐lateral slip on an approximately north‐northeast–south‐southwest striking, or left‐lateral movement on a roughly west‐northwest–east‐southeast‐striking nodal plane. The azimuths of neither nodal plane coincide with the trend of the DG fault zone, so it is most likely that one of the sub‐faults was activated during the 2023 Szarvas earthquake swarm. To identify the small‐scale tectonic structure of the source region, however, requires further research using various geophysical and geological data.

Recently, Békési et al. (2023) and Porkoláb et al. (2024) have updated the stress database in the Pannonian region and reconstructed the neotectonic stress field in the area. In the epicentral region of the 2023 Szarvas earthquake swarm, the authors suggest an approximately east‐northeast–west‐southwest‐striking maximum‐horizontal stress direction. The azimuth of the P‐axes obtained in this study varies between 70° and 77°, which is consistent with the present‐day stress orientation reported by Békési et al. (2023) and Porkoláb et al. (2024).

We hope that the results presented in this article will help to better understand the otherwise poorly understood small‐scale tectonic structure of the source region. We also believe that beyond regional interest, this work contributes to the broader seismological community by illustrating that the probabilistic BayesLoc and MCMT methods are versatile tools that can be successfully used for detailed local event analysis. Our results may encourage further applications of the methods in other small‐scale seismicity studies worldwide.

This article used waveform data from the following seismic networks: CR—Croatian Seismograph Network (doi: 10.7914/SN/CR), University of Zagreb, Croatia; CZ—Czech Regional Seismic Network (doi: 10.7914/SN/CZ), Institute of Geophysics, Academy of Sciences of the Czech Republic, Czech Republic; HU—Hungarian National Seismological Network (doi: 10.14470/UH028726), Institute of Earth Physics and Space Science (EPSS) Kövesligethy Radó Seismological Observatory (KRSO), Hungary; OE—Austrian Seismic Network (doi: 10.7914/SN/OE), Central Institute for Meteorology and Geodynamics (ZAMG), Austria; RO—Romanian Seismic Network (doi: 10.7914/SN/RO), National Institute for Earth Physics, Romania; SJ—Serbian Seismological Network (doi: 10.7914/SN/SJ), Seismological Survey of Serbia, Serbia; SK—Slovak National Seismic Network (doi: 10.14470/FX099882), Geophysical Institute, Slovak Academy of Sciences, Slovakia; SL—Seismic Network of the Republic of Slovenia (doi: 10.7914/SN/SL), Slovenian Environment Agency, Slovenia; UA—Carpathian seismological Network of Ukraine, Carpathian Branch of the Institute of Geophysics, National Academy of Sciences of Ukraine, Ukraine; Z6—AdriaArray Temporary Network: Albania, Austria, Czech Rep., Germany, Hungary, Kosovo, Montenegro, Slovakia (doi: 10.7914/2cat-tq59); Y8—AdriaArray Temporary Network: Bulgaria, Moldova, Poland, Romania, Ukraine (doi: 10.7914/b1sc-0n71). The data have been obtained through the European Integrated Data Archive (EIDA) available at https://www.orfeus-eu.org/data/eida/. Data access to the Z6 and Y8 AdriaArray temporary networks is restricted. Information about AdriaArray is available at https://orfeus.readthedocs.io/en/latest/adria_array_main.html. For historical earthquake data, we used the EPSS Earthquake Catalogue. Regional moment tensor (MT) solutions were taken from the public catalog of the German Research Centre for Geosciences (GFZ; https://geofon.gfz-potsdam.de). The synthetic Green’s functions (GFs) were computed using the computer program hspec96, version 3.3, developed by Robert Herrmann, Department of Earth and Atmospheric Sciences, Saint Louis University. This program is contained in Computer Programs in Seismology (Herrmann, 2013), a software package available at https://www.eas.slu.edu/People/RBHerrmann/CPS330.html. Earthquake locations were estimated using the BayesLoc software package (https://gs.llnl.gov/nuclear-threat-reduction/nuclear-explosion-monitoring/bayesloc). Figures were prepared using the Generic Mapping Tools (GMT) software, version 4.5.8 (Wessel and Smith, 1998). All websites were last accessed in December 2024.

The authors acknowledge that there are no conflicts of interest recorded.

The reported investigation was financially supported by the National Research, Development and Innovation Fund, Hungary (Research Grant Hungary, RGH_L 151351). The authors are grateful to the operators of the permanent seismic networks that make their data available through the European Integrated Data Archive (EIDA). The authors would like to thank all members of the AdriaArray Field Team, especially T. Czifra and B. Süle at the Kövesligethy Radó Seismological Observatory (KRSO) for the installation and maintenance of the temporary AdriaArray stations used in this study. The invaluable work of the AdriaArray Seismology Group is also highly appreciated. In addition, the authors would like to thank the two anonymous reviewers and Associate Editor Cezar Trifu for their insightful comments and suggestions, which have significantly helped them to improve this article.