Relocation of earthquakes in the southern and eastern Alps (Austria, Italy) recorded by the dense, temporary SWATH-D network using a Markov chain Monte Carlo inversion

In this study, we analyzed a large seismological dataset from temporary and permanent networks in the southern and eastern Alps to establish high-precision hypocenters and 1-D V_P and $V_{\mathrm{P}}/V_{\mathrm{S}}$ models. The waveform data of a subset of local earthquakes with magnitudes in the range of 1–4.2 M_L were recorded by the dense, temporary SWATH-D network and selected stations of the AlpArray network between September 2017 and the end of 2018. The first arrival times of P and S waves of earthquakes are determined by a semi-automatic procedure. We applied a Markov chain Monte Carlo inversion method to simultaneously calculate robust hypocenters, a 1-D velocity model, and station corrections without prior assumptions, such as initial velocity models or earthquake locations. A further advantage of this method is the derivation of the model parameter uncertainties and noise levels of the data. The precision estimates of the localization procedure is checked by inverting a synthetic travel time dataset from a complex 3-D velocity model and by using the real stations and earthquakes geometry. The location accuracy is further investigated by a quarry blast test. The average uncertainties of the locations of the earthquakes are below 500 m in their epicenter and ∼ 1.7 km in depth. The earthquake distribution reveals seismicity in the upper crust (0–20 km), which is characterized by pronounced clusters along the Alpine frontal thrust, e.g., the Friuli-Venetia (FV) region, the Giudicarie–Lessini (GL) and Schio-Vicenza domains, the Austroalpine nappes, and the Inntal area. Some seismicity also occurs along the Periadriatic Fault. The general pattern of seismicity reflects head-on convergence of the Adriatic indenter with the Alpine orogenic crust. The seismicity in the FV and GL regions is deeper than the modeled frontal thrusts, which we interpret as indication for southward propagation of the southern Alpine deformation front (blind thrusts).

1 Introduction

The Alpine orogen resulted from the collision of the European and Adriatic plates (Dewey et al., 1973; Schmid et al., 2004; Handy et al., 2010, 2015). This collision led to localized deformation patterns within the plates and along the suture, often accompanied by seismic activity. The seismicity of the Alps has been investigated both on the scale of the Alpine orogen (e.g., Diehl, 2008) and on the scale of parts of the chain (e.g., the southern Alps; Amato et al., 1976; Cagnetti and Pasquale., 1979; Galadini et al., 2005; Anselmi et al., 2011; Scafidi et al., 2015; Bressan et al., 2016; Slejko, 2018; Viganò et al., 2008, 2013, 2015). The southern and eastern Alps have heterogeneous earthquake distribution, with regions of high activity, e.g., the eastern southern Alps (ESA); rare activity, e.g., the eastern Alps; or inactivity, e.g., parts of western southern Alps (WSA) (Fig. 1).

The dense temporary seismic SWATH-D network was deployed from 2017 to 2019 in the southern and central Alps (Heit et al., 2017, 2021) with a total number of 161 stations. Given the average station spacing of only 15 km, this network was designed to capture the shallow crustal seismicity (especially the depths of the earthquakes) with higher precision than the coarser permanent networks, including the AlpArray network (AASN, with 52 km station spacing; Hetényi et al., 2018), and to obtain an increased spatial resolution of images (e.g., local earthquake tomography, receiver functions), even in a region with low or moderate local seismicity. The SWATH-D network is located in a key part of the Adriatic indenter, for which a switch in subduction polarity was proposed at the transition from the central to eastern Alps (e.g., Lippitsch et al., 2003).

In this study, we use a subset of 335 local earthquakes recorded by the SWATH-D network and complemented by a selection of 112 AlpArray network stations nearby to precisely relocate the hypocenters of the local earthquakes. The dense, high-quality travel time picks created in this study potentially lead to constrained hypocenter solutions with high internal consistency. This will enable us to identify the general pattern of seismicity on the surface and at depth throughout the region and contribute to the understanding of active tectonic processes. A further aim of the study is to derive a high-quality dataset suitable to be used in local earthquake tomography (LET).

Although locating earthquakes has been a routine task in seismology for decades, there are several challenges related to obtaining a precise location. One of them is the trade-off between the hypocenters and the velocity structure (so-called hypocentre–velocity coupling; Kissling, 1988; Thurber, 1992; Kissling et al., 1994). To yield accurate locations (especially depths), either the velocity model should be well known in advance or, particularly in the case of earthquakes occurring within a network (local earthquakes), it has to be simultaneously inverted for (Kissling, 1988; Kissling et al., 1994; Thurber, 1983; Husen et al., 1999). This is conventionally being done by employing iterative inversion strategies based on damped least squares. These methods are quite robust and have been successfully applied for years. However, because they use a linearization and the damped least-squares approach, they depend not only on proper choices of initial values for hypocenter coordinates (close enough to the true location) and the velocity model(s) but also on model parametrization (i.e., layers) or regularization (i.e., damping factors), which all have to be carefully checked and selected.

Recently, a transdimensional, hierarchical Bayesian approach utilizing a Markov chain Monte Carlo (McMC) algorithm was implemented for simultaneous inversion of hypocenters, 1-D velocity structure, and station corrections, specifically for the local earthquake case (Ryberg and Haberland, 2019). This approach has the advantage of being largely independent of prior knowledge of model parameters (starting velocity model, starting locations of earthquakes), and the inversion results are data-driven. Moreover, the results can be statistically analyzed, and thus errors and ambiguities can be estimated. The method extends the probabilistic relocation approaches (Lomax et al., 2000) by inverting for a set of velocity models well explaining the data.

2 Tectonic setting and local seismicity

The study region (Fig. 1) is part of the collisional zone between the Adriatic and European tectonic plates. The convergence of these plates began no later than 84 Ma, with the onset of collision at ca. 35 Ma (Dewey et al., 1973; Schmid et al., 2004; Handy et al., 2010) and culminating since ca. 24 Ma with north-northwestward indentation of the Adriatic continental lithosphere (Handy et al., 2015). Indentation is still ongoing today (Cheloni et al., 2014; Aldersons, 2004; Serpelloni et al., 2016; Sánchez et al., 2018).

Deformation in the target area (red box in Fig. 1) occurred within a corner of the Adriatic indenter, which is delimited to the north by the Periadriatic Fault (PAF) and to the west by the Guidicarie Fault (GF). The PAF is a late orogenic fault active in Oligo–Miocene time, which was sinistrally offset by the GF in Miocene time. Just north of the northwestern corner of the indenter, where the GF and PAF meet, the Tauern Window (TW) exposes remnants of the European lower plate at the surface (Scharf et al., 2013; Schmid et al., 2013; Favaro et al., 2017; Rosenberg et al., 2018). Therefore, shortening and exhumation within the TW are attributed to faulting along the Adriatic indenter (Scharf et al., 2013; Favaro et al., 2017; Reiter et al., 2018). Counterclockwise rotation of the Adriatic plate with respect to the Eurasian plate (Le Breton et al., 2017) is interpreted to have induced seismic activity in the southern Alps (Anderson and Jackson, 1987; Mantovani et al., 1996; Bressan et al., 2016).

The thrust front in the eastern southern Alps (ESA) accommodates most of the Adria–Eurasia convergence with pronounced N-oriented present-day motion (Serpelloni et al., 2016; Sánchez et al., 2018). A N-oriented horizontal deformation of ∼2 mm a⁻¹ is observed in the ESA front (northern part of the Venetian-Friuli Basin; Sánchez et al., 2018; Métois et al., 2015; Cheloni et al., 2014) with highest recorded geodetic strain rate along the Venetian Front (within the Montello region, around 45.8^∘ N, 12.25^∘ E) (Serpelloni et al., 2016; Verwater et al., 2021). From the Venetian part of ESA toward the Po Basin (PoB), the motion vectors have a slight westward rotation with decreasing magnitudes. In contrast, a progressive eastward rotation of motion vectors is observed from the Venetian area toward the Pannonian Basin. The entire Alpine chain is still undergoing uplift, with a maximum in the inner Alps, e.g., along the border between Switzerland, Austria, and Italy (near 46.6^∘ N, 11^∘ E) with values of ≥ 2 mm a⁻¹, whereas uplift rates decrease toward the foreland (Sánchez et al., 2018).

Two sedimentary basins, the Molasse and the Po basins (MoB and PoB in Fig. 1) dominate the shallow structure at the northern and southern borders of the target area. The MoB, running along the northern front of the Alps, thickens going across strike from the exposed Variscan basement (0 km) to the Alpine Front (5 km; see, e.g., depth-to-basement contours in the map of Bigi et al., 1989). The Quaternary and Tertiary sediments of the PoB thicken from 0 km at its northern margin to about 6 km in the south below the Apenninic front (Bigi et al., 1989; Waldhauser et al., 2002).

Figure 1 Geological map of the Alps showing major deformation fronts, faults, geographical subdivisions, tectonic units, and seismicity (seismic events between January 2000 and December 2012 with a magnitude larger than 3 from the European-Mediterranean Seismological Centre). The light green zone indicates the outline of the Adriatic plate. The tectonic units of the study area (red box) are shown on the upper-right inset map. The maps are modified from Handy et al. (2010) and Schmid et al. (2004, 2008). Abbreviations are as follows: WA stands for western Alps, CA stands for central Alps, SA stands for southern Alps, EA stands for eastern Alps, ADR stands for Adriatic plate, IVR stands for Ivrea body, PB stands for Pannonian Basin, VFB stands for Venetian-Friuli Basin, PoB stands for Po Basin, MoB stands for Molasse Basin, TW stands for Tauern Window, PAF stands for Periadriatic Fault, NGF stands for northern Giudicarie Fault, SGF stands for southern Giudicarie Fault, ESA stands for eastern southern Alps, WSA stands for western southern Alps, MT stands for Montello, GL stands for Giudicarie–Lessini region, FV stands for Friuli-Venetia region, and AN stands for Austroalpine nappes.

The seismicity in the study area and surrounding regions has been monitored routinely for many years with stations operated by different national and regional networks ( National Institute of Geophysics and Volcanology, INGV, and National Institute of Oceanography and Experimental Geophysics, OGS, in Italy; Central Institution for Meteorology and Geodynamics in Austria, ZAMG; and Swiss Seismological Service in Switzerland, SED). Besides routine catalogs, the seismicity has also been investigated in a few seismological studies through various time periods and regions (Blundell et al., 1992; Solarino et al., 1997; Galadini et al., 2005; Diehl, 2008; Viganò et al., 2013; Reiter et al., 2018, among others). Diehl (2008) compiled the local earthquake data from 14 seismic networks in the greater Alpine region from 1996 to 2007 and created a uniform and consistent relocated event catalog. The SHARE European Earthquake Catalog (Grünthal et al., 2013) comprised homogeneous earthquakes from 1000 to 2006.

These studies indicate that seismic activity is clustered in the Friuli-Venetia (FV), Giudicarie–Lessini (GL), and the Inntal regions (Fig. 1). The FV region is located along the thrust-and-fold belt marking the active Adria–Europe plate boundary and was struck by several earthquakes of M_W≥6 over the last several centuries (Slejko, 2018). The energy was released on a system of E–W south-verging thrusts (some of which are blind or unknown), as well as on backthrusts and oblique-slip faults (Nussbaum, 2000; Galadini et al., 2005; Slejko, 2018; Romano et al., 2019). The GL region is located in the deformation zone along the western margin of Adriatic indentation with two major fault and fold systems characterized by high seismic activity: the Giudicarie Belt and the Schio-Vicenza strike-slip fault (Viganò et al., 2015). In contrast, the area east of the northern Guidicarie Fault (NGF) has obviously been much less active since at least 1994 (Viganò et al., 2015).

3 Data acquisition and processing

3.1 Network data

For this study, we used mainly data of the SWATH-D network, which was temporarily deployed in a roughly rectangular region in northern Italy and southeastern Austria (black box in Fig. 2). The installation started in Summer 2017, and the network ran for almost 2 years (Heit et al., 2017, 2021). The network consisted of 151 broadband seismometers with an average inter-station spacing of 15 km. A total of 78 stations equipped with Güralp-3ESPC seismometers (https://www.guralp.com/products/surface, last access: 1 November 2020) and Earth Data EDR-210 recorders (http://www.earthdata.co.uk/edr-210.html, last access: 1 November 2020) transmitted data in near real-time via the cellular network (online). The other 73 stations, with data accessed during station services (offline), were equipped with two different kinds of seismometers and digitizers: 58 of them with Nanometrics portable Trillium Compact seismometers (https://www.nanometrics.ca/products/seismometers/trillium-compact, last access: 1 November 2020) and CUBE recorders (https://www.gfz-potsdam.de/en/section/geophysical-deep-sounding/infrastructure/geophysical-instrument-pool-potsdam-gipp/pool-components/seismic-pool/recorder-dss-cube3/, last access: 1 November 2020), and the rest with Güralp-3ESPC seismometers and Earth Data PR6-24 recorders (http://www.earthdata.co.uk/pr6-24.html, last access: 1 November 2020). In Autumn 2018, the network was complemented with an additional 10 online stations in the northeastern edge by LMU-Munich. All SWATH-D seismic stations recorded the data continuously at 100 samples per second (sps).

In order to enlarge the azimuthal coverage for earthquakes occurring in the periphery of the SWATH-D framework, we additionally selected 112 stations of the larger-scale AlpArray Seismic Network (AASN; Hetényi et al., 2018) to include their waveform data in the analysis.

Figure 2 Distribution of seismic stations of the SWATH-D network (red, blue, purple) and selected stations of the AASN (white). The black box indicates the periphery of the SWATH-D network. The faults are from Schmid et al. (2004, 2008) and Handy et al. (2010).

3.2 Seismic events, arrival time picks, and phase-type identification

The national seismological agencies of Italy (INGV and OGS), Austria (ZAMG), and Switzerland (SED) provide comprehensive earthquake catalogs (minimum magnitude −0.8 by ZAMG) in our study region. Therefore, the process of event detection was skipped and an integrated event list from the national agencies, after removal of common events, formed a proper list of 2639 local events (see Appendix A for more information) for arrival time picking. Considering this large number of events and 273 selected permanent and temporary stations, we applied a modified version of the automated multi-stage workflow from Sippl et al. (2013) to the waveform data for picking the arrival times of P and S waves. This workflow was originally implemented for producing a complete catalog of well-located earthquakes and reliable arrival times using continuous passive seismic data in the Pamir–Hindu Kush region (Sippl et al., 2013). We, subsequently, assessed the performance of the picking algorithms by quality-checking of time and phase-type uncertainties of the picks (Appendix B).

As described in Sect. B2, the number of mispicks and picks ignored by the automatic procedure was quite large. Therefore, a dataset with relatively highly constrained events was selected for visual (manual) inspection. Since the number of picks in this stage (after the automatic procedure) is not the most reliable information on the location precision (see Sect. B2) and also considering that this data is aimed for high location accuracy of the occurring earthquakes, we selected events with azimuthal gap less than 200^∘ and root mean square (rms) less than 1 s, which resulted in 384 event (18 390 P and 7762 S picks). This visual (manual) inspection of this dataset was carried out by removing or modifying obvious mispicks (mainly at large epicentral distances) and also re-picking the missed arrivals (mostly S picks at small epicentral distances).

This manual (visual) inspection was mainly performed on individual station data. The direct Pg (Sg), Moho-reflected PmP (SmS), and Moho-refracted Pn (Sn) phases arrive closely spaced in time, especially in the triplication zone, potentially leading to phase-type misidentification. The PmP (SmS) is always a secondary arrival, but its amplitude can dominate the first arrival and thus be easily mispicked. On the other hand, for epicentral distances larger than the triplication zone, the Pn (Sn) amplitude can be so weak that the Pg (Sg) or PmP (SmS) be wrongly identified as the first arrival. In the local earthquake case, and particularly for regions with varying Moho topography and significant lateral variations in the crustal structure, something to be expected for the Alps, the phase-type identification is even more challenging (Diehl et al., 2009b). Since the inversion is based on the first arrival times we looked at the waveforms in epicentral distance plots for larger events for which we expect the later phases, inspected the picks, and corrected or adjusted where needed. However, we would like to point out that even with these procedures we cannot rule out a certain amount of misidentified phases.

An inspection for identifying potential quarry blasts (and other anthropogenic sources) was simultaneously performed during the visual inspection. Criteria for potential quarry blasts were (1) relatively small number of S-picks, (2) relatively small $\mathrm{S}/\mathrm{P}$ amplitude ratios (see, e.g., Walter et al., 2018), and (3) large surface waves (observed dispersive waveform characteristics). Based on our assessment, we classified the events into 344 earthquakes, 15 potential blasts, and 25 unclear events.

A Wadati diagram (Wadati, 1933; Kisslinger and Engdahl, 1973; Diehl, 2008) was used to calculate the $V_{\mathrm{P}}/V_{\mathrm{S}}$ ratio of the earthquakes, which is 1.72 (calculated by linear least squares; Fig. 3). Furthermore, S picks with S–P travel time differences larger than 4 s off the main trend were removed (however, this was only 0.3 % of all picks). After this removal, only 2 % of the whole observations have $t_{\mathrm{s}}-t_{\mathrm{p}}>\mathrm{0.72}{t}_{\mathrm{p}}\pm \mathrm{3}$ s. We notice that at P travel times larger than ∼ 25 s the observations tend to slightly larger S–P travel time differences, potentially indicating a higher $V_{\mathrm{P}}/V_{\mathrm{S}}$ ratio at larger depth, i.e., in the upper mantle.

Figure 3 Modified Wadati diagram based on P and S picks after visual (manual) analysis (the color indicates number of picks; please see the color palette table on the right). The origin times (t₀) of the earthquakes are provided by the single-event localization using HYPO71 hypocenter determination program (Lee and Lahr, 1975). The solid black line shows the fit to the data by linear least squares, and its slope indicates the $V_{\mathrm{P}}/V_{\mathrm{S}}$ ratio of 1.72. The dashed lines show values with $t_{\mathrm{s}}-t_{\mathrm{p}}=\mathrm{0.72}{t}_{\mathrm{p}}\pm \mathrm{3}$ s.

For the simultaneous inversion of hypocenters, 1-D velocity models, and station corrections, we used the events with at least 10 P picks and 5 S picks and an azimuthal gap less than 180^∘, which comprises 301 local earthquakes. The arrival time dataset consists of 11 084 P picks and 6496 S picks (quality classes of 0 to 3) with average picking errors of 0.12 and 0.21 s for P and S observations, respectively (Table 1). The average number of picks per event is 32 for P and 18 for S, with maximum of 157 and 118, respectively (Fig. 4).

Table 1 Statistical parameters of high-quality P and S pick dataset of 301 well-constrained local earthquakes used for simultaneous inversion of hypocenters, 1-D velocity model, and station corrections.

Figure 4 Histograms of high-quality P and S picks per event for 301 well-constrained local earthquakes used for simultaneous inversion of hypocenters, 1-D velocity models, and station corrections.

In order to investigate the seismicity pattern, we use the travel time data of an additional 43 local earthquakes with slightly fewer picks (minimum of eight P picks and two S picks) to the main dataset and relocate them using the V_P, $V_{\mathrm{P}}/V_{\mathrm{S}}$, station corrections, and noise levels from the simultaneous inversion (Sect. 6.3). This updated dataset contains 344 earthquakes, 12 420 P picks, and 7192 S picks.

4 Method

4.1 Probabilistic Bayesian inversion

Considering that hypocenter locations and velocity model are inherently linked to each other (coupled hypocenter-velocity problem, e.g., Kissling, 1988; Thurber, 1992; Kissling et al., 1994), especially in the local earthquake case, a simultaneous inversion for hypocenters and velocity structure (and/or station corrections) is needed. In this study, different to the conventional approach of damped least squares, we use a Bayesian approach (Bayes, 1763). Bayesian approaches have been applied in a number of geophysical studies (Tarantola et al., 1982; Duijndam, 1988a, b; Mosegaard and Tarantola, 1995; Gallagher et al., 2009; Bodin et al., 2012a, b; Ryberg and Haberland, 2018). Ryberg and Haberland (2019) recently implemented a hierarchical, transdimensional Markov chain Monte Carlo (McMC) approach for the joint inversion of hypocenters, 1-D velocity structure, and station corrections for the local earthquake case. While with classical inversion techniques a discrete best-fitting model m (i.e., locations, velocity model, etc.) is derived, in the Bayesian inference all the parameters of the model m are represented probabilistically. For the theoretical background the reader is referred to Bodin and Sambridge (2009) and Bodin et al. (2012a, b). We define a uniform and wide range of values for each individual model parameter (before knowing the observed data; prior) and after performing the inversion for random combinations of the parameters, we derive a probability distribution individually for each model parameter (after combining the prior information with the observed data; posterior). Thus, a large suite of models is generated, all of them fitting the travel time observations. The choice of a uniform and extensively wide range of model parameters (Table 2) guarantees that the final model is dominated by the data rather than by the prior information (i.e., starting model, choice of parameters).

Table 2 Prior distribution (model space), starting model, and width of Gaussian distribution for each model parameter.

4.2 Model parametrization and forward problem

Hypocenters are described by spatial (x₀,y₀,z₀) and temporal (t₀) coordinates and the 1-D velocity structure is described by a variable number of horizontal layers with individual values of V_P and $V_{\mathrm{P}}/V_{\mathrm{S}}$. Moreover, the model m comprises station corrections for P and S waves (τ^P and τ^S), which account for travel time effects (delayed or earlier arrivals) due to deviations of the 1-D model from the real 3-D velocity structure in the shallow subsurface beneath the stations. The inversion also derives the quality of the data expressed as noise levels for P and S picks (σ_P and σ_S), which are also unknowns in the inversion. The noise includes actual picking errors, systematic errors of observations or measurements, and any approximation or errors of forward travel time calculation (see below) and is assumed to be normally distributed and uncorrelated. Therefore, the complete model to be inverted for is defined as follows:

$$\begin{array}{}\text{(1)}& m=[K,{V_{\mathrm{P}}}_i,V_{\mathrm{P}}/{V_{\mathrm{S}}}_i,h_i,x_j,y_j,z_j,t_j,{{\mathit{\tau }}^{\mathrm{P}}}_k,{{\mathit{\tau }}^{\mathrm{S}}}_k,{\mathit{\sigma }}_{\mathrm{P}},{\mathit{\sigma }}_{\mathrm{S}}],\end{array}$$

where K is the number of layers; h is layer depth; and i, j, and k refer to layer, earthquake, and station index, respectively. It should be noted that the number of layers is also unknown. In other words, the model space dimension is not fixed in advance and hence the posterior is a transdimensional function.

In order to sample models from the prior probability distribution, we apply McMC approach to create a huge number of different models m by random selection of model parameters (see Eq. 1 and Sect. 4.3). Accordingly, to compute travel times the forward problem has to be solved a very large number of times. We use the 2-D finite differences (FD) Eikonal solver (Podvin and Lecomte, 1991) for the forward calculation. This algorithm requires a regular mesh. Therefore, the irregular velocity model is converted to a fine and uniform mesh by setting the velocity at each mesh point to the value of the nearest point from the irregular model (V_Pi or $V_{\mathrm{P}}/V_{\mathrm{S}}$_i). The fine mesh used by the Eikonal solver has a cell spacing of 1 km vertically and horizontally.

Once the travel times are calculated, a misfit function, particularly for each model, is defined as the summed squared differences between the observed (d) and calculated travel times. The misfit function is then used to build the Gaussian likelihood function and the posterior values of the proposed model (for detailed information one could refer to Ryberg and Haberland, 2019).

4.3 Markov chain Monte Carlo (McMC) algorithm

Given the Bayesian approach described above, the McMC algorithm generates an ensemble of models with parameters within the prior distribution. We mainly follow the hierarchical, transdimensional procedure proposed by Ryberg and Haberland (2018, 2019) that supports the calculation of both model parameters and model dimensionality. The evolution of a model along the Markov chain consists of four main steps: (1) choose a random initial model m (Table 2). (2) Generate a new model from the prior distribution by perturbing the current model parameters (changing one of the velocity parameters of a random node, shifting the position of a random earthquake, adding a new layer, removing a random layer, altering one of the noise parameters, changing a randomly selected station correction. The changes in the V_P, $V_{\mathrm{P}}/V_{\mathrm{S}}$, cell position (i.e., layer), and noise levels must be according to Gaussian probability distribution centered at the current value. The values of the new model are randomly selected within very wide bounds (Table 2) and thus do not truncate the posterior probability distribution. (3) The forward calculation (Sect. 4.2) is then performed on the new model. The newly estimated travel time data are compared with the observed data d and then the misfit function, data likelihood, and the posterior probability are determined. (4) The newly proposed model is accepted or rejected based on the criteria of Bodin and Sambridge (2009). If the new model is rejected, then the current or old model is retained by reiterating step 2. However, there is still a probability of acceptance even if the fit of the new model is worse than that of the old model. If the model is accepted (i.e., when it is better than the previous model) then it acts as a starting model for another iteration (step 1).

By reiteration of steps 1 to 4, a chain of models is produced, which is in fact the Markov chain. This chain is continued until the misfit is no longer significantly decreasing (burn-in phase). Thereafter, a stationary model space sampling is achieved. If these sequences are repeated long enough, a chain provides an approximation of the posterior distribution for the model parameters. To accelerate the model space sampling, up to 1000 separate and independent chains are investigated in parallel.

The main idea of the McMC method is to keep prior assumptions of the model parameters (conventionally used as initial values) at a minimum and thus minimize unwanted artificial effects on the final model. Therefore, the McMC method only uses the travel times and does not directly depend on initial hypocenters, origin times, velocity models, or even the model parametrization (e.g., grid node spacing). Nevertheless, the selection of picks in the (semi-automatic) picking procedure involves the comparison of the observed travel times with those based on preliminary hypocenters (which depend on initial velocity models). Therefore, the results of the inversion can sensu stricto depend to some degree on initial values.

5 Uncertainty estimates using synthetic tests

As a variable Moho topography and complex lithosphere structure are expected for our study region, the 1-D model derived by and used in our inversion can only be a rough abstraction of the real conditions. Nevertheless, this simplification could potentially have influence on earthquake location correctness and accuracy. To assess the performance of our inversion strategy in this respect, using a 3-D velocity model and earthquake and station distribution of the real data, we calculated the synthetic travel times, added synthetic noise and finally inverted them in the same way as the real data (see below). Comparing input event locations (synthetic) and the inverted (output) ones allows us to study the recovery of the hypocenters, location consistency, and potential systematic errors related to the use of a 1-D model, which we can generally expect for the derived real hypocenters. For example, it can be studied whether events at the periphery or in certain parts of the model have systematically larger uncertainties (e.g., due to their location and/or spatial distribution of picks). Furthermore, we can study how the (output) 1-D model looks in comparison to the (input) 3-D model, how large the derived noise is in relation to the synthetic input noise, and how the pattern of station corrections corresponds to the shallow velocity anomalies. Similar tests are standard in structural studies (i.e., LET) to study the recovery of certain features.

In our synthetic 3-D background model, the course of the crust–mantle boundary is inspired by the Moho topography of the European and Adriatic plates from Spada et al. (2013); however, we adopted it in a simplified way without features such as crustal stacks (overlapping Mohos of the two involved plates) or the proposed Moho hole. Our main aim of this simplified adaption is studying the influence of a 1-D model class on the derived hypocenters after the inversion and the goal is not recovering the exact velocity structure. We assigned typical crustal velocities to the upper layer, starting with 6 km s⁻¹ and increasing to 7 km s⁻¹ at the base of this crustal layer. Thereafter, the velocity increases from 7 to 8 km s⁻¹ within a thin zone of 3 km with a strong gradient. Following this, the velocity increases again gradually until it reaches 8.05 km s⁻¹ at 90 km depth (Kennett et al., 1995; Christensen and Mooney, 1995). The Moho depth in the target area changes from 22 km along the flank of the Ivrea body to more than 55 km along the Adria–Europe plate boundary (Spada et al., 2013). We considered the model's surface to be −5 km to avoid wave propagation through the air.

In order to simulate realistic 3-D effects in the context of expected crustal structure, we superimposed shallow high- and low-velocity anomalies onto the background 3-D model presented above. We assumed that the hard rocks of the European basement exhumed in the TW have higher velocities than the surrounding Austroalpine nappes and thus form a high-velocity anomaly (+5 %) within a rough outline of the TW and from the surface to the depth of 10 km (for more information on the TW, see, e.g., Bertrand et al., 2015). Furthermore, we designed lower-velocity anomalies for the PoB and MoB sedimentary basins (the coarse outlines of the basins are taken from Waldhauser et al., 2002). An average V_P of 4.35 km s⁻¹ is considered for the entire MoB from the model's surface down to 2 km depth. For the PoB, a layer of 6 km thickness starting from the model's surface and with V_P of 2 km s⁻¹ less than the background model is constructed. The shallow anomalies are then inserted in the background model. Taking the shallow anomalies into account not only makes the model more realistic but also allows us to check whether the station corrections reflect travel time deviations due to velocity variations in the shallow subsurface. Figure 5 displays the 3-D synthetic V_P model from different views and also the 1-D model at two different viewpoints. The $V_{\mathrm{P}}/V_{\mathrm{S}}$ ratio in the entire model is set to $\sqrt{\mathrm{3}}$.

The P and S arrival times are calculated using the 3-D FD Eikonal solver (Podvin and Lecomte, 1991; Tryggvason and Bergman, 2006). We generated an FD grid with an increment of 1 km horizontally and vertically, resulting in a grid dimension of 601 × 321 × 96 with a total number of 18 520 416 nodes. For simulating a realistic travel time dataset, we adopted the geometry of stations and hypocenters of the real data. Using 344 earthquakes, in total 12 534 P and 7258 S synthetic picks resulted from the forward calculation. Thereafter, according to the manual arrival time uncertainties in Table B2, we added random noise (0.05, 0.1, 0.2, 0.3, 0.1, 0.2, 0.3, and 0.4 for P-0, P-1, P-2, P-3, S-0, S-1, S-2, and S-3, respectively) to the arrival times, and this dataset was then used for the simultaneous inversion by McMC.

Figure 5 The 3-D synthetic V_P velocity model. The $V_{\mathrm{P}}/V_{\mathrm{S}}$ ratio was fixed to $\sqrt{\mathrm{3}}$ in the entire model. (a) and (b) Contours of the Moho and shallow anomalies of the PoB, MoB, and TW from different view directions. (c) Modified Moho depth map based on Spada et al. (2013). (d) A 1-D velocity model of two different locations with shallower and deeper Moho retrieved from the 3-D synthetic velocity model (the locations are shown with red and black crosses in (a) and (c), respectively). The faults are the same as those in Fig. 2.

To explore the model parameters, 100 Markov chains each with 1000 random initial models are used for simultaneous inversion of McMC. Following the strategy in Ryberg and Haberland (2019), we derived ∼ 15 000 of the final models from the Markov chains (more information is given in Sect. 6) to define model parameters based on the average μ and standard deviation σ. For the earthquake epicenters (x and y) and station corrections (τ^P and τ^S) the classical average of the Gaussian distribution is used. However, the depth distribution of shallow earthquakes is truncated by the upper model boundary. In order to accommodate this, we used an algorithm based on truncated Gaussian distributions (Ryberg and Haberland, 2019) to derive true depth averages and uncertainties. The V_P and $V_{\mathrm{P}}/V_{\mathrm{S}}$ values are defined based on the modified average (Ryberg and Haberland, 2018) and standard deviation.

Figure 6 represents the recovery of the earthquakes. Figure 6a shows that the epicenters are recovered very well (the blue dots are located on top of red circles). Further assessment of the earthquake recovery (Fig. 6b) demonstrates that the average misfit between synthetic and recovered hypocenters is quite close to zero with uncertainties of ∼ 400 m in epicenter, 1.31 km in depth, and 0.12 s in origin time. Nevertheless, we note that although the epicenters of most of the earthquakes throughout the study area are recovered well, the depths of the events at the edges of the network are less well recovered and that their uncertainties are also systematically larger (see Fig. 6). The recovered noise levels, which are representative of the unresolved part of the travel times, are calculated for individual pick types and quality classes separately (Fig. 6c). They are close to the random noise, which was added to the synthetic arrival times (manual arrival time uncertainties in Table B2); however, they are slightly higher. This deviation could be explained by the forward modeling errors and the fact that the inversion derives a 1-D velocity model from data of a 3-D input model.

Figure 6 (a) Recovery of the hypocenters after synthetic test. Red circles are recovered earthquakes with a 1σ error bar (error bars of latitude and longitude are invisible at this scale), and blue dots are the synthetic locations. The gray triangles show the stations. (b) Histograms of the misfit between recovered hypocenters and original ones. The hypocenters are recovered very well, with an average misfit of 20 m (±380 m) in longitude, 210 m (±400 m) in latitude, 110 m (±1.131 km) in depth, and 0.12 s (±0.08 s) in origin time. (c) The recovered noise levels after McMC for individual pick types and quality classes are slightly higher than the random noise, which was added to the synthetic travel time data.

Figure 7a shows the derived V_P and $V_{\mathrm{P}}/V_{\mathrm{S}}$ models as heat maps (posterior distribution of all the models). In addition, it shows the average value (red line), standard deviation of all the inferred models (yellow zone), and a reference model with maximum likelihood (blue line), i.e., maximum posterior probability (similar to Ryberg and Haberland, 2019). As seen, no clear single Moho velocity jump is recovered, but it is modeled by a gradual increase of the velocities at depths from 30 to 50 km. We attribute this to the variable Moho topography, which cannot be modeled by the 1-D velocity model per se, and potentially lower resolution due to less ray coverage in this depth range. Because of a dense ray coverage from the surface down to a depth of ∼ 20 km, the V_P uncertainty is almost zero in this depth range. In the depth interval with strong Moho topography, the uncertainty varies between 0.1 and 0.66 km s⁻1. The uncertainty values are compatible with the standard deviation of the lateral heterogeneity in the synthetic 3-D velocity model (∼ 0.8 km s⁻1 above 2 km, ∼ 0.05 km s⁻1 between 2 and 20 km, and between 0.1 and 0.5 km s⁻1 below 20 km depth). Thus, we think that the increased standard deviations beneath ∼ 30 km depth indicating reduced and fading resolution are due to both variable Moho topography and reduced ray coverage; however, their exact contribution is hard to discern.

The $V_{\mathrm{P}}/V_{\mathrm{S}}$ ratio was fixed to the square root of 3 for the whole region in the forward modeling, and the same value is recovered down to ∼ 35 km depth with almost zero uncertainty (less than 0.02). The uncertainty of $V_{\mathrm{P}}/V_{\mathrm{S}}$ increases below 35 km depth to a maximum value of 0.046, indicating that this deeper part of the model is not well resolved.

The station corrections (Fig. 7b) correspond to a large extent to the shallow anomalies of the 3-D synthetic model (Fig. 5). The shallow high-velocity anomaly in the TW is expressed by negative corrections (large circles) corresponding to earlier arrivals, while the shallow low-velocity anomalies of the PoB and MoB correspond to positive corrections (large crosses) reflecting later arrivals.

Figure 7 (a) Recovery of the V_P and $V_{\mathrm{P}}/V_{\mathrm{S}}$ models after the synthetic test. The results are shown with heat maps (probability histogram vs. depth; warm and cold colors correspond to higher and lower probabilities, respectively) with gray areas showing zones with no model achievement. The modified average model (solid red line), most probable model (solid blue line), and corresponding uncertainty of 1σ (yellow zone) for both V_P and $V_{\mathrm{P}}/V_{\mathrm{S}}$ models are shown as well. The $V_{\mathrm{P}}/V_{\mathrm{S}}$ ratio of $\sqrt{\mathrm{3}}$ is displayed with the dashed green line. (b) Recovered P wave station corrections of all the stations that were involved in the McMC inversion. Blue crosses show positive values reflecting lower velocities, and red circles display negative values indicating higher velocities than expected. Regions indicated in gray have shallow high- or low-velocity anomalies (see Fig. 5 for more information). Symbol size corresponds to correction amplitude. The faults are same as Fig. 2.

6 Results and discussion

For the simultaneous inversion of the hypocenters, the 1-D velocity model, and the station corrections, we used the travel time dataset of the earthquakes (Sect. 3.2). In the simultaneous inversion, for the first 300 000 models of the Markov chains, only the hypocenters are allowed to be changed. Thereafter, V_P, $V_{\mathrm{P}}/V_{\mathrm{S}}$ ratio, station corrections, and noise levels additionally are permitted to change, adding 300 000 more models to the Markov chain. After inversion of 400 000 models in each Markov chain, the model space reached the stationarity (end of burn-in) phase. Thereafter, badly converging chains are automatically excluded, and every 1000th model is stored from the remaining models after the burn-in phase. The stored models of each chain are combined (in total ∼ 15 000 models) to calculate the posterior probability. The final models have stationary rms misfits of 0.36 s (Fig. 8).

Figure 8 On the left histogram, the rms misfit along Markov chains during the inversion evolution is shown (color indicates number of models; please see color palette table on the left). Before the red line (300 000 models), only the hypocenters are perturbed and beyond that the V_P, $V_{\mathrm{P}}/V_{\mathrm{S}}$ ratio, station corrections, and noise levels are allowed to change as well. The green line indicates the end of the burn-in phase after ∼ 400 000 models. On the right histogram, the rms misfit for all the models after burn-in phase is shown. The best-fitting models are characterized by an average 0.36 s rms misfit.

6.1 The 1-D velocity model and station corrections

It turns out that the algorithm needs roughly five layers to best fit the data (Fig. 9a). Similar to the synthetic test (Fig. 7), the final velocity models are shown with the average of all remaining models and the model with maximum posterior probability. In general, for both V_P and $V_{\mathrm{P}}/V_{\mathrm{S}}$ models (Fig. 9b) the average value (red line) varies gently, whereas the model with maximum posterior probability (blue line) is rather coarse and contains discontinuities.

The derived V_P model starts with a rather high value of 5.94 km s⁻¹ at the surface and increases only very gently down to a depth of around 20 km. After a zone of higher velocity gradient in the average model around (and a step-like increase of the reference model at) ∼28 km depth, velocities reach ∼6.8 km s⁻¹ at ca. 30 km (down to ca. 45 km depth). Down to this depth the uncertainties (1σ) indicating the resolution increase from 0.01 km s⁻¹ (very good) to 0.28 km s⁻¹ (fair). Our derived V_P model in the upper crust – except for the uppermost layer shallower than 3 km – is similar to the Diehl et al. (2009b) model, an increase for values above 6.5 km s⁻¹ is observed deeper than in the latter model.

Below ca. 45 km depth, velocities increase further up to 7.85 km s⁻¹ at 60 km depth; however, standard deviations (1σ) between 0.3 and 0.42 km s⁻¹ indicate only poor resolution in this depth range. A single step-like increase from “crustal” to “upper mantle” velocities (as at 35 km depth in the Diehl et al. (2009b) model) cannot be observed, which we attribute either to (1) the expected strong depth variability of the Moho in the study area (see the discussion of the synthetic model in Sect. 5) and a gradual course of the 1-D model due to lateral averaging in this depth range, (2) the relatively small number of Pn arrivals, and/or (3) limited maximum observation distances. Influences by wrongly picked late-arriving Pg arrivals at large distances (instead of Pn phase) cannot totally be ruled out.

As a 1-D model cannot be representative of the 3-D structure, especially in a region with expected irregular Moho and complex crustal structure, a geologically meaningful interpretation of the derived V_P model is hardly possible.

The $V_{\mathrm{P}}/V_{\mathrm{S}}$ model starts with high values at the surface (to ∼ 5 km depth) and shows reduced values of around 1.70 down to 22 km depth before reaching 1.77 at greater depths (with uncertainty between 0 and 0.02). However, 1σ values of 0.02 to 0.045 indicate only poor resolution at depth larger than around 30 km, which we associate to only a few Sn arrivals. This was basically expected from the Wadati diagram (Fig. 3) and is in agreement with values previously derived, e.g., by Viganò et al. (2013). Both V_P and $V_{\mathrm{P}}/V_{\mathrm{S}}$ models are available in the Supplement S1.

The station corrections derived from the McMC simultaneous inversion potentially indicate local, shallow 3-D velocity anomalies in the subsurface, which cannot be accounted for by the 1-D model. The McMC inversion assumes that P and S station corrections have an average of zero. Negative corrections indicate earlier wave arrival and thus higher velocities in the (shallow) subsurface, whereas positive correction implies delayed arrival that is indicative of lower velocities.

The pattern of corrections (Fig. 9c) shows coherent negative corrections associated with the EA, ESA, and CA. The large negative values in the eastern part (east of 15^∘ E) might be related to proximity to the edge of the network and thus be dominated by mantle phases (faster arrivals). Surprisingly, in between the negative values in the Alpine Chain, a pattern of slight positive corrections in the WSA is observed. Besides, extreme positive corrections are seen in the PoB and MoB as expected for sedimentary basins, which is also consistent with the results of the synthetic test (see Sect. 5). A pronounced pattern specifically related to the TW as seen in our synthetic test is not observed so clearly in the real data, suggesting a different structure in the shallow subsurface. The pattern of stations in the ESA and a few stations in the WSA and CA agree very well with results by Diehl et al. (2009b).

A detailed interpretation of the pattern of the station corrections is highly ambiguous because they contain a superposition of site effects and/or effects from 3-D structural variations, as well as a lot of smearing. Nevertheless, the general concept of simultaneously inverting local earthquake datasets for a (simplified) 1-D velocity model, hypocenter positions, and origin times and reactions proved to be very powerful for accurately localizing earthquakes (see e.g., Kissling, 1988). Moreover, the general pattern of corrections can indicate the consistency of phase data.

Figure 9 (a) On the left histogram, the probability of the number of layers of the random models introduced by Markov chains during the evolution are shown (color indicates number of models). On the right histogram, the number of layers of the models after burn-in phase is given. (b) V_P, $V_{\mathrm{P}}/V_{\mathrm{S}}$, and V_S models after McMC simultaneous inversion. Figure characteristics are similar to Fig. 7. The velocity models are recovered down to ∼ 63 km depth. (c) P wave station corrections of all the stations that were involved in the inversion corresponding to the 1-D velocity model in (b). Negative corrections (red circles) indicate earlier arrivals (indicative for higher velocities in the shallow subsurface) and positive corrections (blue crosses) indicate delayed arrivals (representative for lower velocities underneath the station). The faults are same as in Fig. 2. Abbreviations are as follows: CA stands for central Alps, EA stands for eastern Alps, ESA stands for eastern southern Alps, WSA stands for Western Southern Alps, MoB stands for Molasse Basin, and PoB stands for Po Basin.

6.2 Estimation of hypocenter accuracy based on relocation of quarry blasts

To validate the localization procedure, the detected 15 quarry blasts (based on visual inspections; see Sect. 3.2) were independently relocated by McMC using the previously derived 1-D velocity model and station corrections from the simultaneous inversion (Sect. 6.1). Figure 10 focuses on the blast distribution associated with two quarry areas close to the villages of Albiano (Italy) and Gummern (Villach, Austria). After the relocation of the blasts, we see that the epicenters are located within the quarry areas and that the depths are in the range of the quarry topography (considering the average uncertainty of 1σ), some of them are offset by a maximum of hundreds of meters. Based on the average mislocation of the blasts (relative to the centers of the quarry areas) and also their location uncertainties, we estimate the absolute location errors in the range of 1 km horizontally and 500 m vertically. This indicates that although the number of picks is generally lower for blasts, the McMC routine provides high-precision hypocentral solutions. It is expected that the errors of earthquakes (which are potentially deeper than blasts) are smaller than those estimated from this test because they are less affected by the heterogeneous shallow structure, which was poorly accounted for by the model (Kissling, 1988; Husen et al., 1999).

The hypocentral solutions of the blasts, in comparison with those obtained by INGV/ZAMG, have an average difference of ∼ 160 m in longitude, ∼ 1 km in latitude, and ∼ 4.6 km in depth, and thus the depths are considerably better resolved and recovered in our study due to the availability of the denser SWATH-D network.

Figure 10 (a) Map view showing the location of quarry blasts relocated using single-event mode of McMC, with solutions of the 1-D velocity model and station corrections from the simultaneous inversion of earthquakes. (b) and (c) Epicenter and depth distribution of the blasts associated with two quarries close to Albiano in Italy and Gummern in Austria, respectively. The light gray band in the cross sections shows the surface elevation variation within the map view boundary. The dark gray zones show the location and surface elevations of the quarries. The depth histogram is also shown for each quarry (the bar in the depth histograms displays the surface elevation within the quarry area). Please note that elevations above the sea level are shown with negative values.

6.3 Seismicity distribution

The distribution of seismicity in the southern and eastern Alps is shown in map view and cross sections in Figs. 11 and 12. This includes the same 301 events used for the simultaneous inversion (Sect. 6.1), as well as 43 additional earthquakes with slightly lower quality (see Sect. 3.2). For this relocation, we used the modified average V_P and $V_{\mathrm{P}}/V_{\mathrm{S}}$ models, station corrections, and noise levels that resulted from the simultaneous McMC inversion. After relocation, the depth probability distribution of nine earthquakes was not Gaussian, and thus (due to their high depth uncertainty) we decided to remove these earthquakes from our catalog. Based on statistical analyses of the McMC inversion, the average epicentral and depth uncertainties (1σ) are ∼ 500 m and ∼ 1.7 km, which are compatible with precision estimates by the synthetic test (Sect. 5). However, the absolute location errors estimated by quarry blasts test (Sect. 6.2) are ∼ 1 km for epicenter and ∼ 500 km for depth.

Since the national agencies are (probably) using much less data for the location (smaller number of stations used, larger inter-station distances), a significant difference between their hypocenter solutions and those obtained in this study is expected (average of 2.4 km in epicenter and 3.7 km in depth). The earthquake depths calculated by McMC are systematically shallower than those by national agencies (by an average of 1.1 km). The maximum and minimum differences in the epicenters and depths (between McMC and national catalogs) are seen for the earthquakes from the INGV and SED, respectively.

The derived hypocenters in this study are not a representative seismicity catalog for the region (the national catalogs contain also many small, poorly constrained events in a much longer period) but form excellent data for further seismological studies, e.g., local earthquake tomography (LET). Moreover, these highly precise hypocentral data allow further tectonic inferences. The hypocenter solutions of constrained earthquakes derived in this study are available in the Supplement S2.

The seismicity is clustered in the same areas as in previous seismic studies of the region (e.g., Reiter et al., 2018). The long-term seismicity by, e.g., merged national catalogs and the SHARE catalog (Grünthal et al., 2013) also indicates the similar seismicity pattern. Most activity is seen within the orogenic retro-wedge in the FV region, with somewhat lesser activity in the SW part of the Giudicarie Belt (GL region). Further activity is located in the Austroalpine nappes north of the PAF and the region around Innsbruck (upper and lower Inntal and the Stubai Alps). The earthquakes reach a maximum depth of ∼ 20 km (Fig. 11), most of them occurring in the depth range of 5 to 15 km. Regions of little or no seismicity are observed at the northeastern corner of the eastern Adriatic or Dolomites indenter, southeast of the NGF (e.g., Reiter et al., 2018), and north of the PAF (Fig. 11).

The overall pattern of seismicity reflects the head-on convergence of the Adriatic indenter with the Alpine orogenic crust, accommodated along thrust faults and folds in the FV region and segmented in the western part of the ESA by strike-slip faults of the Giudicarie Belt.

Figure 11 Distribution of the 335 well-located earthquakes in the target area after McMC inversion (single-event mode). Circles are earthquakes that were incorporated in the simultaneous inversion, and pentagons are 43 additional earthquakes with fewer picks. The color and size of the circles indicate the depth and magnitude of the earthquakes, respectively (magnitude is taken from national catalogs). The dashed black box indicates the periphery of the SWATH-D network. Purple stars show the location of Friuli earthquake of 6 May 1976 (M_L 6.4) and its major aftershock on 15 September 1976 (M_L 6.1) (Slejko, 2018), and the blue star is the Salò earthquake on 24 November 2004 (M_L 5.2; Viganò et al., 2015). Black lines indicate cross sections shown in Fig. 12. Faults are compiled from geological maps of the Italian geological survey, including Avanzini et al. (2010); Bartolomei et al. (1967); Bosellini et al. (1967); Braga et al. (1971); Castellarin et al. (1968, 2005); and Dal Piaz et al. (2007) and complemented by fault traces published in Handy et al. (2010) and Schmid et al. (2004, 2008). Abbreviations are as follows: ESA stands for eastern southern Alps, WSA stands for western southern Alps, FV stands for Friuli-Venetia region, GL stands for Giudicarie–Lessini region, AN stands for Austroalpine nappes, PAF stands for Periadriatic Fault, NGF stands for northern Giudicarie Fault, SGF stands for southern Giudicarie Fault. The histogram on the right side indicates the depth distribution of the earthquakes in the crust.

Figure 12 Seismicity distribution in depth along three profiles of A, B, and C in Fig. 11. For better clarity, the depth and length scales of cross section A are magnified by a factor of 1.5 (without any distortion) compared to cross sections B and C. Red dots indicate hypocenters, projected from a swath extending to 15 km on either side of the profiles. Purple stars show the location of Friuli earthquake of 6 May 1976 (M_L 6.4) and its major aftershock on 15 September 1976 (M_L 6.1) (Slejko, 2018), and the blue star is the Salò earthquake on 24 November 2004 (M_L 5.2; Viganò et al., 2015). Geological structures of the FV region along cross section A are modified following Nussbaum (2000) with the blind Trasaghis thrust taken from Galadini et al. (2005). Geological structures NW and SE of the NGF in cross section B are modified after Rosenberg and Kissling (2013) and Verwater et al. (2021), respectively. Abbreviations are as follows: NGF stands for North Giudicarie Fault, CA stands for Campofontana Fault, PT stands for Priabona–Trambileno Fault, and CAL stands for Calisio Fault.

6.3.1 Friuli-Venetia (FV) region

In accordance with the previous studies (and references therein Bressan et al., 2016; Reiter et al., 2018) and national catalogs, most of the seismicity in our dataset occurs in the ESA, i.e., within the FV region, coinciding with the eastern part of the deformation front between the Adriatic plate and the PAF. This was also the location of the destructive Friuli earthquake of 6 May 1976 (M_L 6.4) and its major aftershock on 15 September 1976 (M_L 6.1) (purple stars in Figs. 11 and 12a and c), which are associated with the Susans–Tricesimo Fault and buried or “blind” Trasaghis thrust, respectively (Poli et al., 2002; Galadini et al., 2005; Slejko, 2018). We find most earthquakes around 13^∘ E (close to the villages of Tolmezzo and Gemona) down to a depth of around 17 km (Fig. 12a and c – around km 220, corresponding to 13^∘ E). However, the seismicity is distributed over a wider area (see also Fig. 12c eastern part; between km 140 and 280). As stated above, a direct connection of individual earthquakes in our dataset to known faults near the surface is difficult. Nevertheless, the clustering of seismicity in the FV region between the Alpine frontal thrusts and the Fella–Sava fault in the north suggest that several frontal thrusts and backthrusts are active. Moreover, Petersen et al. (2021) shows a dominant E–W striking thrust mechanism for this region.

According to the fault distribution (and naming) of Galadini et al. (2005), the seismicity along the cross section A (Fig. 12a) indicates that the Susans–Tricesimo and Trasaghis faults (and potentially the Maniago thrust) are the most active faults in this region. Most earthquakes are located below the Maniago thrust, one of the main Alpine frontal thrusts in the tectonic model of Nussbaum (2000). This suggests that there is an active fault at greater depth (Fig. 12a). We interpret this to indicate a blind (i.e., that has not reached the surface) southward-propagating thrust front. In the profiles of Galadini et al. (2005) and Merlini et al. (2002) further east, this deeper blind thrust (named Trasaghis thrust in Galadini et al., 2005) reaches the surface, where it offsets Plio–Pleistocene sediments. More near-surface activity down to a depth of 6 km can generally be related to the Friuli thrust system (see Nussbaum, 2000, and Fig. 12a).

6.3.2 Giudicarie-Lessini (GL) region

The seismicity in the GL region correlates with the NNE–SSW transpressive fault system of Giudicarie Belt and the NW–SE-trending Schio–Vicenza strike-slip fault system. Seismicity of the Giudicarie fault system is concentrated mainly in the south (southwest of Lake Garda at depths ranging from a few kilometers down to 15 km; Fig. 11) and is in agreement with Viganò et al. (2015). For this part of the Giudicarie Belt, Viganò et al. (2008) suggested the kinematic regime of thrust with a strike-slip component (see also Petersen et al., 2021). The big Salò earthquake on 24 November 2004 with M_L=5.2 and a thrust mechanism (blue star in Figs. 11 and 12c; Viganò et al., 2015) was suggested to have occurred on a low-angle thrust connected to the steep Ballino–Garda Fault, which detaches the sedimentary cover (hanging wall) from the underlying crystalline basement in its footwall (Viganò et al., 2008, 2013, 2015).

The NNW–SSE trending system of vertical strike-slip faults making up the Schio–Vicenza Fault is located between the active NW-dipping thrusts of the Giudicarie Belt to the west and the NNW-dipping thrusts of the Valsugana thrust system to the east. The short observational period of our network and the low seismicity rate limits the unambiguous association of the earthquakes to specific faults. However, some of the captured seismicity seems to take place along the Campofontana (CA), Priabona–Trambileno (PT), and Calisio (CAL) faults (Fig. 12c – western section). The Schio–Vicenza fault system is the most active structure in the GL region and was shown with predominant right-lateral strike-slip mechanism (Viganò et al., 2008).

Figure 12b shows the distribution of seismicity superposed on the main geological structures, which indicates that seismicity is focused in its southern part of the GL region. This suggests that the most southern and deepest faults (Fig. 12b) remain active, while the more internal fault systems have become seismically inactive (Verwater et al., 2021). Moreover, we observe that the seismicity is deeper than the frontal thrust modeled from balanced geological cross sections (Fig. 12b; Verwater et al., 2021), similar to observations within the FV region (as mentioned above). We interpret this to reflect southward propagation of the southern Alpine deformation front (blind thrusts) towards the Po Plain.

6.3.3 Lateral variations in clustering of seismicity from the WSA to ESA

A cross section running orogen-parallel from the GL to the FV region (Fig. 12c) indicates a similar depth of earthquakes for both regions. The FV region shows relatively high seismic activity, located at the junction with the Dinaric Front (Doglioni and Bosellini, 1987). However, in the central part of cross section C (around km 150 in Fig. 12c), sparse seismicity may indicate a seismic gap (Anselmi et al., 2011; Burrato et al., 2008). Alternatively, the sparse seismicity (especially at shallow depths) could indicate that deformation within this area is occurring aseismically, as proposed by Barba et al. (2013) and Romano et al. (2019) for the Montello thrust (Fig. 11). Strain rates within the Montello region are among the highest within the ESA (Serpelloni et al., 2016), which combined with its sparse seismicity indicates that the majority of deformation within this area is most likely occurring aseismically (Barba et al., 2013).

6.3.4 Engadine and Austroalpine nappes

The seismicity in the Austroalpine nappes (close to Italy–Switzerland border, Fig. 11) reaches a rather homogeneous depth around 11 km and does not seem to coincide with known geological structures in the area (e.g., Engadine Fault, Fig. 12b, left part). Although earthquakes are situated at a uniform depth, the distribution of these events in map view is quite diffuse (Fig. 11). Note that there are no seismicity clusters around the Engadine Fault in cross section B (Fig. 12b, left part), although the map view shows that seismicity coincides with the fault trace of the Engadine Fault NE and SW of cross section B. This could be related to block rotation along the Engadine Line (Schmid and Froitzheim, 1993).

6.3.5 Inntal region

A cluster of seismicity is found in the Inntal region at a depth range of 5 to 18 km (Fig. 11). This region is known for earthquakes with strike-slip and oblique-slip mechanisms that have been attributed to the activity of the Inntal Fault (Reiter et al., 2018; Petersen et al., 2021) and the Alpine basal thrust (Reiter et al., 2003). Historical earthquakes (e.g., the earthquake of 17 July 1670 in Hall in Tyrol; Hammerl, 2017) were also reported. Current seismicity (Fig. 11) extends further south into the Brenner region, where it has been associated with activity of the Brenner normal fault (Reiter et al., 2005).

6.3.6 Other structures

Several earthquakes in a depth range of 5 to 14 km are clearly aligned along the PAF between 12 and 12.5^∘ E, possibly indicating ongoing activity of this fault (see also Reiter et al., 2018). The ENE–WSW Valsugana fault system does not seem to be seismically active, similar to the conclusion that Viganò et al. (2015) reached based on his analysis in the time frame of 1994–2007.

7 Conclusions

In this study, a semi-automatized picking procedure was applied to the waveform data recorded by SWATH-D and selected AlpArray network stations of earthquakes listed in national seismic catalogs. During the visual inspection of the picks, several quarry blasts were identified and separated from the earthquakes. In order to establish a very high-quality, consistent earthquake catalog in the study region, the arrival time data of the most constrained events (in terms of azimuthal gap, rms, and number of picks) were inverted simultaneously for hypocenters, 1-D velocity models, and station corrections using a McMC approach. A synthetic recovery test, using a close-to-realistic 3-D velocity model, validated the reliability of the method for yielding highly precise hypocenters even in complex 3-D velocity conditions. A quarry blast test demonstrated that the blast hypocenters are located very close to the quarry areas with average horizontal and vertical mislocations (relative to the centers of the quarry areas) of 740 and 376 m, respectively. The average uncertainties calculated by McMC for the earthquake locations are ∼ 500 m in epicenter and ∼ 1.7 km in depth. The derived V_P and $V_{\mathrm{P}}/V_{\mathrm{S}}$ models are well resolved down to a depth of 30 km with an uncertainty below 0.45 km s⁻1 and 0.05, respectively. At greater depth the resolution decreases. Below 30 km depth, velocities gradually increase to mantle velocity, which is reached at around 50 km depth. The station corrections demonstrate coherent negative values associated with the EA, ESA, and CA and slight positive values in the WSA. Moreover, the PoB and MoB represent extreme positive corrections, as is expected for sedimentary basins.

The study is based on a very dense and large station network; however, it is limited to 16 months. Nevertheless, our precise hypocenter solutions with high internal consistency allow us to have insight into the active tectonic situation throughout the study region. The seismicity pattern confirms previous seismotectonic studies based on long-term records or local networks. Seismicity is found in diffuse clusters in the upper crust (0–20 km) and within broader zones of the Alpine frontal thrust, e.g., the FV region, along the GL and Schio-Vicenza domains, and in the internal Alpine and Inntal areas. The northeastern corner of the Adriatic or Dolomites indenter is obviously inactive or low active. We think that the overall pattern of seismicity reflects the head-on convergence of the Adriatic indenter with the Alpine orogenic crust, accommodated in the FV region and along the Giudicarie Belt. The seismicity in the FV and GL regions shows activities deeper than the modeled frontal thrusts, which we interpret as an indication for southward propagation of the Southern Alpine deformation front (blind thrusts). The seismicity in FV and GL indicates similar earthquake depths for both regions, while the sparse seismicity between them may indicate that the deformation within this area is most likely occurring aseismically.

The 1-D V_P and $V_{\mathrm{P}}/V_{\mathrm{S}}$ models and the reassessed precise hypocenters will be utilized for investigating 3-D seismic velocities in the region using LET. The earthquake database is provided in Supplement S1 and S2.

Appendix A: Integration of the seismic event catalogs

The process of integrating events from the national seismological agencies is summarized in Fig. A1. It starts with getting events between September 2017 and the end of 2018 and confined to our study area with M_L ≥ 1 from catalogs of INGV, OGS, ZAMG, and SED. Thereafter, three individual event lists of A for Italy, Slovenia, and Croatia from INGV; B for Germany and Austria from ZAMG; and C for Switzerland from SED are selected. The lists will be subsequently updated by the other three catalogs (one at a time, as described in Fig. A1). Eventually, the three lists are merged into one list, which totals 2639 local events in the study region.

Figure A1 Flow chart showing the process of merging seismic catalogs of INGV, OGS, ZAMG, and SED in order to derive an integrated seismic event list in the study region.

Appendix B: Automatized arrival time picking and earthquake relocation

B1 Workflow

Our modified workflow consists of three main stages: (1) the NonLinLoc ray tracer (Lomax et al., 2000) for predicting initial P arrival times based on preliminary origin time and hypocenters (integrated earthquake catalog from national agencies) and the minimum 1-D velocity model of Diehl (2008); (2) the MPX algorithm (Aldersons, 2004; DiStefano et al., 2006) for improving the first alerts of P picks and classifying them into four quality classes (0 to 3, Table B2) based on Fisher statistics (Diehl, 2008, Appendix D: Users guide for MPX picking system); and (3) the Spicker algorithm (Diehl et al., 2009a) for picking S onsets based on P picks, event location, and a preliminary S pick (determined by NonLinLoc ray-tracing; Lomax et al., 2000). Spicker also rates picks into four quality classes (Table B2).

After each stage, in order to remove the effect of low-quality picks on hypocenter locations, a multiple relocation procedure using HYPO71 (Lee and Lahr, 1975) and various subsets of picks is performed. The procedure starts with localization of the events using the most reliable picks (higher-quality classes and closer distance to the event), and this first location is considered to be the reference for testing further picks. Thereafter, the picks with lower quality and those from distant stations are incorporated one-by-one. The newly added pick is kept if its residual is less than a distance-dependent, user-defined value and also if the root mean square (rms) residual of the event location is still acceptable (less than 1 s). Picks leading to unacceptable rms residuals are removed from the dataset. Moreover, in this relocation procedure, not only low-quality picks but also events with a small number of high-quality picks are removed. This procedure results in the location not being dominated by low-quality picks or not having sorted out better picks. Since MPX and Spicker depend on the correctness of the predicted arrival times (the predicted arrival time itself depends on the event location), picks and location are highly inter-dependent. This can be controlled by iteratively improving the event location, re-picking arrival times based on the updated location, and adding or removing dubious picks in each step.

The implementation of this automatic workflow for the event list by national catalogs (Appendix A) resulted in a dataset with 2639 events, 68 099 P picks, and 17 151 S picks. The statistics of the picking quality classes are given in Table B1.

Table B1 Distribution of the picks and their quality classes after implementation of the automatic picking algorithms for the event list created by national catalogs.

B2 Uncertainty assessment of arrival time picks

To assess the performance of the picking algorithms, a dataset consisting of 14 earthquakes with various magnitude ranges that are randomly scattered in the region are manually analyzed. Based on the arrival time uncertainty, the P and S picks are qualified into five classes (0, 1 ,2, 3, and 4, meaning 100 %, 75 %, 50 %, 25 %, and 0 % contribution into each location, respectively), and they are considered as reference for calibration (Table B2).

Table B2 Statistics of the automatic and reference picks associated with the quality classes for a dataset of 14 earthquakes.

* The mean of misfits between automatic and reference picks.
** The standard deviation of misfits between automatic and reference picks.

The misfits between reference and automatic arrival times are calculated (individually for each quality class of the P and S picks) and the mean and standard deviation values are shown in Table B2. The picks with higher qualities have smaller mean values, which means their misfits are lower (i.e., automatic picks are closer to reference picks). Moreover, except for S-3, the standard deviation of each class is very comparable with manual arrival time uncertainty.

Another way of assessing the performance of the automatic picking algorithms is the so-called confusion matrix. The confusion matrix, in general, shows the performance of predicted vs. the actual data of different classifications in machine learning (Kohavi and Provost, 1998). For our dataset, it shows automatic picks vs. reference picks (similar to DiStefano et al., 2006; Diehl et al., 2009a; Sippl et al., 2013). In the confusion matrices displayed in Table B3, three values represent the performance of the automatic pickers in each quality class: (1) number of picks, (2) hit rate (fraction of reference picks of a certain class that are assigned to a specific class by the automatic picking) that is summed up to 100 % in each row, and (3) the standard deviation of the residuals between automatic and reference picks. The red and orange colors refer to picks that are upgraded to higher qualities by the automatic algorithms (unfavorable), whereas the green sectors show the picks that are underestimated by the automatic picker (more eligible). Having low hit rates in the orange and red cells and high hit rates in the green fields point toward a conservative assessment of the picking algorithms. The white fields in the confusion matrices illustrate the picks that are ideally qualified based on the reference dataset.

Table B3 Confusion matrices of P and S picks, comparing the automatic picks with reference picks. Columns are quality classes of the automatic P picks from MPX (a) and S picks from Spicker (b), and rows represent the pick weighting classes provided by the human analyst for reference data. The parameters in each cell of the matrices are explained in the text.

The confusion matrices of both MPX and Spicker show relatively high values of hit rate (alternatively N) in some orange and red cells that we interpret as mispicks. For instance, the MPX (Spicker) picked 157 P picks (166 S picks) that are either ignored or qualified as class 4 by a human analyst. As the human analyst adopted class 4 for far stations with unclear arrivals, far automatic picks are unreliable. On the other hand, it seems that Spicker ignored a large number of picks (high hit rates in the most right column of the Spicker matrix). This is eminent for picks with manual quality class 0 that are potentially associated with the closest stations based on human analysis.