Introduction

In July–August 2019, the area of Lake Constance (Southern Germany, near the borders with Switzerland and Austria) was hit by a series of moderate magnitude earthquakes. The strongest earthquake, felt by the population (I0 = IV), occurred at 23:17 UTC of 29 July 2019, with magnitude Ml = 3.7. Other three events exceeded Ml = 3.0.

The geological complexity of the Alpine region affects the hypocentral location of these events, because of significant anomalies of travel times observed at regional distances (i.e. 150–400 km) (see e.g. Piromallo and Morelli 2003; Wagner et al. 2013; Spada et al. 2013; Mayor et al. 2016). Anomalies in wave amplitude attenuation were also observed in this region (Babuška et al. 1990).

The purpose of this paper is the study of the effects of lateral crustal heterogeneities on the quality of earthquake locations, and the ways these effects can be moderated. Here we focus on the 2019 Constance Lake earthquake sequence as a good example of these problems.

We show that, when using arrival times picked from regional stations, and adopting a simple 1-D velocity, the effect of velocity anomalies produces relevant mislocations. These mislocations can be significantly reduced by the application of various alternative methods, such as 3-D velocity models, station specific models or, more simply, source-station specific corrections. Relative location methods, as well as the double-difference joint hypocentral determination (DDJHD) method, can be suitably used for obtaining high-resolution spatial images of the seismogenic structure in an arbitrary reference system.

Geological setting of the crust in the alpine region

The construction of the Alps is the result of a very complex geological history, which began in the late Cretaceous with the collision between the African and European plates (Argand 1924; Capitanio and Goes 2006; Handy 2010) that is still ongoing today. The NNW-SSE convergence of Europe below the Adria microplate (Piromallo and Faccenna 2004) led to the progressive closure of the Ligurian-Piedmontese Ocean and the lithospheric substratum of these oceanic relicts has been subducted appearing today as positive P-wave velocity anomalies in different tomographies (e.g. Spakman et al. 1993; Worthel and Spakman 2000; Piromallo and Morelli 2003). As a result of the collision, the Earth's crust thickened from an average of 30 km up to 60 km (Fitzsimons and Heinz 2001), producing a complex lithological and structural patterns associated with the development of a series of overthrusted nappes. Massive shortening phenomena led to the superimposition of rocky masses, originally occupying specific paleogeographic positions, moved even for several hundred km, to form a double vergence orogen with northward thrusting in the northern part and southward thrusting in the southern one. The boundary between these two areas is an important tectonic line, known as the "Insubric Line'' (Fig. 1), roughly divides the Alps into two parts, separating the Southalpine to the south from the Pennidic and the Australpine to the north. Accommodation of collisional shortening varies along strike of the alpine chain (Rosenberg and Kissling 2013). While the backbone of the European vergence domains is given by metamorphic rocks, the rocks of the southern Alps are essentially sedimentary, with a prevalence of carbonatic rocks (limestone and dolomite). Intrusive rocks make up the Adamello Massif, while volcanic rocks characterize the Atesino volcanic complex and the Euganean Hills.

Fig. 1
figure 1

Map of the seismic stations used in this study. Blue dots are stations receiving Pg waves as first arrivals, adopting the intermediate Moho depth velocity model and a hypocentral depth of 5 km. Red dots denote stations receiving Pn phases as first arrivals. The epicenter of the mainshock of 29 July 2019 is shown by a red cross. The red lines show the main thrusts and the tectonic lineaments. EH: Euganean Hills; AM: Adamello Massif; AVC: Atesino Volcanic Complex

Full size image

The complexity of the alpine crustal structure just described is reflected on the velocity model in this area. In this respect, Wagner et al. (2013) developed a 3-D P-wave velocity model by the inversion of regional secondary phases arrival times. In their analysis they adopted an initial 1-D crustal velocity model as a good representation of the normal continental crust with a near-surface crystalline basement structure, upper and lower crust, an averaged Moho depth of 32 km, and a sub-Moho mantle lithosphere velocity of 8 km s−1. To account for the varying Moho depth, the 1-D reference model was adapted in a second step to the local Moho depth (Fig. 2 in Wagner et al. 2013). In this respect, two additional cases to the normal one were proposed: (i) a shallower Moho (e.g. beneath the Upper Rhine Graben) and (ii) a deeper Moho (e.g. beneath the Alps). In our study, we have verified that the adoption of a shallower Moho or a deeper Moho velocity model for stations located in different areas achieves a better fit between the modeled and the observed arrival times.

Fig. 2
figure 2

Comparison of the epicentral locations obtained for the four considered events with a uniform velocity model (M1) and a variable velocity model (M2) for different stations (see explanation in the text). Depths are fixed at 5 km

Full size image

Data

In this study, we analyze a seismic sequence of 4 seismic events in the Constanza Lake area that occurred in July and August 2019 with magnitude Ml ranging between 3.0 and 3.7. We performed a bulletin search from the International Seismological Centre (ISC) collecting P-waves arrival times from all the nearest 58 available stations. Moreover, we downloaded numerous waveforms from other 28 stations of different European seismological networks (INGV, BGR, GFZ, ZAMG, STN, TSNI, NEISN, RSNI). For these stations we proceeded with a careful manual picking of the first arrival times, which were used for the localization of the events together with the bulletin data. The list of 86 stations with epicentral distances ranging from 0.037° (4 km) to 3.496° (388 km) used in our analysis is reported in Table 1, and their geographical positions are shown in Fig. 1.

Table 1 List of stations used in this study, distance and azimuth from the epicenter of the 19/07/2019 mainshock, number of stations analyzed for each station and average time residuals
Full size table

The geographical distribution of the stations considered in our study includes an area between 4°E and 14°E in longitude and 45°N and 50°N in latitude and the 86 stations are located within 400 km from the epicentral area. The stations are distributed along the Alpine arc in the southern area and in the European foreland to the north. About half of the stations are located at less than 150 km from the epicenter, so that the P-wave first arrivals are represented by both Pg and Pn phases. As it will be shown in the following sections, including Pn arrivals in the data set for hypocentral location may produce a contamination in the arrival times due to the heterogeneity in velocity structures. However, Pn arrivals may have an important role in constraining the determination of hypocentral depths, with respect to the use of Pg arrivals only. Moreover, the choice of considering both Pg and Pn phases is in line with the aim of our study, i.e. to analyze the effect of the Moho heterogeneities. Our main purpose is to explore the issues that arise from using a sparse network of stations to locate events in a tectonically complex region, and to compare our results with the information already available in the literature from tomographic studies of the velocity structure.

Location methods

The seismic location process is essentially the solution of an inverse problem where observed arrival times of seismic waves are compared with theoretical arrival times computed from a given velocity model. Assuming such a velocity model as known, the unknown parameters of this problem are the origin time and the three spatial coordinates of the earthquake source.

The location algorithm adopted in this study is based on a straightforward least-squares iterative best-fit procedure, implemented in a FORTRAN 77 computer code. It was originally developed at the ING (Istituto Nazionale di Geofisica) in the early 1970s with the main purpose of routine processing of seismological observations from the National Seismological Network and earthquake catalog production (see, e.g., Console and Gasparini 1976; Cagnetti and Console 1978; Console and Favali 1981). The very long and extensive use of this code, with some more recent minor improvement, has demonstrated its wide applicability and robustness. The velocity model adopted in this code is a flat 1-D layered earth. In our application, only P-arrivals have been used, as S-arrivals are generally characterized by higher uncertainty in arrival time picking. No automatic quality or distance weighting has been applied, because we don’t want to down-weight the contribution of stations characterized by systematic arrival time residuals. However, when outliers are discovered in the input data, they may be interactively removed by the operator.

It is generally recognized that, independently of the location algorithm, uncertainties in earthquake locations are dominated by two main factors (Husen and Hardebeck 2010):

  1. 1.

    Measurement errors of seismic arrival times,

  2. 2.

    Modeling errors of calculated travel times.

The location uncertainty of type 1 is obtained from the theory of linear matrix inversion, based on the hypothesis of a normal probability distribution for the observed arrival times. It leads to the definition of a probabilistic error ellipsoid, which statistically represents the volume of the space of unknowns within which the true location has a given probability of falling. The projection of the ellipsoid on the horizontal coordinates defines the well-known error ellipse of the epicentral location at a given confidence level. The shape and size of the error ellipse depend on the accuracy of measurement of the arrival times at the seismic stations, and on the azimuth distribution of the stations around the epicenter.

Modeling errors of calculated travel times (uncertainty of type 2) constitutes an additional element of concern for obtaining small location errors. It can be shown that, in real cases, the location uncertainty of type 2 has an effect that shifts the epicenter obtained from the solution of the analytical problem far from the real location by more than the size of the error ellipse. This kind of error depends systematically on the incorrectness of the velocity model applied in the location algorithm, and cannot be reduced by improving the accuracy of arrival time determinations at the seismic stations.

Several methods have been developed to minimize the systematic shifts affecting the hypocentral locations (Husen and Hardebeck 2010). These methods are all based on the fact that such shifts are quite similar for events whose hypocenters are close to each other with respect to their distance from the recording stations. Consequently, these methods can only be applied for a set of seismic events within a restricted volume.

A set of seismic events can be located relatively to a master event the location of which is independently known or just guessed (Poupinet et al. 1984; Console and Di Giovambattista 1987; Mezcua and Rueda 1994). The master-event location uses travel times residuals computed for the master event as station corrections to locate all other events. In the present study, ignoring the real location of any of the earthquakes taken into consideration, we adopted a relative location method, where the unknowns are simply the differences of origin times, horizontal coordinates and depths of two events. This method was implemented in a code, derived from the old ING FORTRAN77 code, where the input data are simply the differences between the arrival times for the two considered events at several couples of stations.

More in general, in order to reduce the systematic errors always present in the location of seismic events, different methods have been developed. Richards-Dinger and Shearer (2000) defined empirical corrections by computing station timing corrections that continuously vary as a function of source position on a local scale. Yang et al. (2001) developed Source Specific Station Corrections (SSSCs) for regional phases and demonstrated that using SSSCs improves the quality of event location. The performance of the SSSC method was confirmed later by Giuntini et al. (2013, 2017) and Materni et al. (2015, 2019). Alternatively, on a regional scale, Myers et al. (2010) have reduced systematic errors by replacing standard travel-time tables, based on 1-D seismic velocity models, by more sophisticated regional or global 3-D models. These models are commonly obtained by means of complex inversion algorithms using sources with accurate hypocenters, such as explosions or earthquakes with hypocenters that are well located by dense regional networks. However, the spatial resolution achievable with these methods is limited by the density of stations and the density of calibration events available in the investigated region. The validation of the results is possible through a comparison of the modeled travel-times with real observations of reference events, independent of those used in the inversion procedure.

For relative location of a group of two or more clustered events, regardless of the absolute location of each of them, the Joint Hypocenter Determination (JHD) method (Douglas 1967), and the Double Difference JHD (DDJHD) method (Waldhauser and Ellsworth 2000; Fisk 2002; Console and Giuntini 2006) have become popular. These methods are implemented to invert large sets of observations coming from as many stations as possible to obtain the coordinates of a set of events in a unique process. They can provide a high-resolution imaging of clustered seismicity in active areas, but they do not resolve the issue of the systematic shift of the cluster as a whole, unless the true location of at least one event in the cluster is known.

Hypocentral locations

The first trial of using the arrival times obtained both from original seismograms and from available bulletins for the four earthquakes described in Sect. 3 consisted in the application of the above mentioned least-squares single-event location algorithm. In this case, we adopted the above-mentioned 1-D average velocity model suggested by Wagner et al. (2013) for Switzerland. The results obtained by our least-squares location algorithm are reported in Table 2 for both fixed-depth (z = 5 km) and free-depth solutions.

Table 2 Event locations for fixed and free depth solutions (using a uniform velocity model)
Full size table

A map of the location obtained for the four events with their 95% confidence ellipses is shown in Fig. 2. As we can infer from both Table 2 and Fig. 2, the accuracy of our locations does not allow us to state that the locations of the four seismic events can be distinguished among each other.

As shown in Table 1, the arrival times for all our four events are characterized by large systematic time-residuals affecting mainly Pn arrivals (at distances larger than 150 km). In particular, German stations to the northeast of the epicenters and French stations to the west exhibit negative residuals larger than one second, while Italian stations located to the southeast are characterized by similarly large positive residuals, of the order of one second or more (Table 1).

Hypothesizing that the residuals shown in Table 1 are due, besides phase picking random errors, to systematic regional anomalies of travel times, we tried a test that consists in assigning, in the location algorithm, different velocity models to selected sets of stations. Therefore, following the indications contained in Table 1 we chose the shallower Moho model for some of the German and French stations (GRA1, GRC1, GRC3, GRB2, GERES, MILB, BOURR, HINF, HAU, CHMF, SFTF, LPL, WOER), and the deeper Moho for some of the Italian and other stations (BORM, CTI, APPI, PTCC, BRES, WLF, BUCH, ABTA) These models, as well as the average one, are described in Wagner et al. (2013). See Fig. 3 for a qualitative sketch of this situation. The results, reported in Fig. 2 and Table 3, show, as expected, an average reduction of about 40% in time residuals and an outstanding improvement of the location accuracy, denoted by a similar reduction of 40% in the size of the error ellipses. The epicenters of the three aftershocks appear located 1–2 km to the West of the mainshock of 29 July.

Fig. 3
figure 3

(Top) Sketch of the N-S cross section of the crustal structure in the Central Europe-Alpine region, and the ray paths from the hypocenter (denoted by the pink star) to seismic stations (denoted by blue and red triangles for Pg and Pn first arrivals respectively); the sketch shows a thinner crust to the North and a thicker crust to the South of the epicenter. (Bottom) Qualitative plot of travel times for Pg and Pn seismic waves according to the above velocity model. The colored dashed lines connect the stations in the top figure with the respective arrival times in the bottom figure. It can be noted that, assuming a uniform P velocity in the crust, Pg arrival times are not affected by the Moho depth, while Pn waves are characterized by later arrivals in the deeper Moho with respect to the shallow Moho at the same epicentral distances

Full size image
Table 3 Event locations for fixed and free depth solutions (using a variable velocity model)
Full size table

Another clear feature of these locations is the systematic shift of the epicenters obtained using a variable velocity model by about 4 km to the SSW of those obtained having fixed a unique velocity model. This is clearly caused by the shallower depth of the Moho of the ray path for the waves travelling to the North with respect to those travelling to the South, and the consequent shorter travel time of the former waves with respect to the latter (Fig. 3).

As to the effect of the choice of constraining the depths to a fixed value of 5 km, or leaving them as a free unknown of the least-squares solutions, both Tables 1 and 3 show irrelevant changes in the epicentral coordinates, while depths span a range of values of 3.8 km–7.4 km for a uniform model (Table 1) and a range of values of only 3.8–4.9 km for a variable model (Table 3).

In order to limit the effect of wave-velocity anomalies on our hypocentral locations, we first applied the relative location method by using the arrival times differences between pairs of events at each station that had recorded at least two events. The results obtained by this method are given in terms of hypocentral coordinates of the 30 July, 31 July and 29 August, having fixed the coordinates of the 29 July earthquake at its single-event solution (47.751°N, 9.110°E and 3.820 km depth). The choice of the 29 July earthquake as master event was suggested by its largest magnitude and presumably, best picking of the P-wave arrivals, but it does not have influence on the relative location obtained for each pair of events. In the application of the algorithm, we discarded from the input data, stations whose arrivals times had a residual larger than three times the standard deviation, assuming that these arrivals were affected by wrong phase picking. The results, reported in Table 4 and shown in Fig. 4, denote much closer locations and smaller size of the error ellipses, letting us infer a small size of the seismogenic structure activated in this sequence.

Table 4 Comparison of relative and JHD locations
Full size table
Fig. 4
figure 4

Relative epicentral locations (blue) and Joint epicentral locations (green) of the three aftershocks of the sequence with respect to the epicenter of the mainshock of 29 July 2019 (red)

Full size image

Lastly, still looking for relative hypocentral locations, we applied the JHD method to the four events considered in this study, using the algorithm developed by Console and Giuntini (2006), but giving as input to this algorithm the time differences obtained from the same time pickings used in the applications of the above described algorithms of single and relative locations. As for the application of the relative location method, also in this case we discarded from the input data stations whose arrivals times had a residual larger than three times the standard deviation. The results obtained by this method, shown in the last column of Table 4 and in Fig. 4, confirm the very small inter-event distances already obtained by the relative location method.

Looking at the epicentral locations reported in the last column of Table 4 and mapped in Fig. 4, we can notice that the locations of the event of 29/08/2019 by the relative and JHD location methods are quite close to each other, while the pairs of the other two aftershocks exhibit a distance of the order of 1 km. We can infer that the discrepancy in epicentral locations obtained by the two methods depends on the circumstance that for the latter two events we discarded different sets of stations. In order to test this hypothesis, we relocated the events by the JHD method discarding the same stations as we had done for the relative location method. The results of this test, as shown in the first of the two columns of Table 4 referring to the JHD method, fully confirm our hypothesis, showing that the results of the relative location and JHD methods are not significantly different, when we use the same input data sets.

Discussion and conclusions

Figure 5 shows a map of Central Europe and Alpine region, where the average time residuals of Table 2 for Pn first arrivals are displayed to show their relationship with the Moho depth (redrawn from Spada et al. 2013).

Fig. 5
figure 5

Map of the stations considered in this study, with the Moho depth of Central-Europe-Alpine region redrawn from Spada et al. 2013. The colors of the station symbols represent the observed average Pn residuals. The red areas represent crustal low Vs anomalies (redrawn from Molinari et al. 2015). The red cross indicates the 29 July 2019 epicenter

Full size image

It is noteworthy that high positive residual times (red circles in Fig. 5) correspond to high crustal thicknesses just over the Alps, while blue circles represent high negative residual times over the shallower Moho European foreland. This circumstance is easily explained by the longer ray path that Pn waves must travel from the hypocenters to the stations in a crust characterized by a deeper Moho (Fig. 3).

In order to analyze the degree of correlation shown in Fig. 5, we plot in Fig. 6 the average Pn residuals versus the Moho depth at the respective recording stations. This figure shows a qualitative but clearly visible correlation. Moreover, we also computed the Pearson correlation coefficient of the regression, obtaining a value R2 = 0.522. For a more quantitative way to assess the goodness of fit for the data set, we applied a statistical test (Kendall 1970), the result of which is a p-value = 1.12·10−6, letting us to definitely reject the null-hypothesis that these data are not correlated.

Fig. 6
figure 6

Average P-wave residuals versus the Moho depth at the respective recording stations and their linear regression (dotted line). The equation shows the linear best fit with its R-squared regression coefficient

Full size image

Both Figs. 5 and 6 show that the residual time values are not the same along the arc of the Alps, even for comparable Moho depths; we have therefore also considered other geophysical parameters. Speranza et al. (2016), calculating the Curie point depth (CPD) by spectral analysis of the aeromagnetic residuals of the Alps, noted as the shallow CPD zones in the Alps lie just above the low P and S low velocity anomalies observed beneath the chain (Di Stefano et al. 2009; Molinari et al. 2015; Lu et al. 2020). These higher time residuals in the eastern Alps correspond to the shallow CPD zones and low velocity zones, also corresponding to the deeper Moho areas.

The most outstanding outlier with respect to the trend of the regression line shown in Fig. 6 is constituted by the case of station WLS. This station lays on a relatively shallow Moho (32 km depth), but it exhibits a positive travel time residual larger than 1 s. We may argue that this residual is at least partly caused by the crustal low velocity anomaly shown in Fig. 5, just to the northwest of the hypocentral zone, which could cause a delay in the Pn seismic waves travelling from the hypocenter to the station, with a ray path crossing the crust to reach the Moho. If we removed the WLS station from the data set, the Pearson correlation coefficient of the remaining data set shown in Fig. 6 would significantly increase up to the value R2 = 0.567. The contrary situation can be recognized in the case of station MBDF, which exhibits a negative residual smaller than −0.8 s, even if it is characterized by a Moho depth of 47 km. We may note that the ray path from the hypocenter to this station crosses a deep Moho crustal structure without any low velocity anomaly.

In the western sector of the Alps, a shallow Moho (Spada et al. 2013) and a strong magnetic signature (Lanza 1975) define the Ivrea body that coincides with the location of the extinction zone of Lg-waves propagating through the Alps as recognized by Campillo et al. (1993). This high-density anomaly has been interpreted as a piece of mantle sandwiched in the lithosphere (ECORS-CROP 1989; Nicolas et al. 1990). Unfortunately, we do not have stations above the Ivrea body to map eventual peculiar velocity structures, but only one at its edges (i.e. MONC), which still register a delay in arrival time (orange circles in Fig. 5). We should, however, consider that the very shallow Moho depth observed in the Ivrea body cannot be considered as the depth of the Alpine crust at which the Pn waves travelled below the Moho to reach station MONC from the hypocenter.

Zooming on the stations that receive Pg-waves as first arrivals, we plotted in Fig. 7 a higher resolution map of the geographical area limited at a distance of the order of 150 km from the epicenters, where the average time residuals of Table 2 for the Pg first arrivals are displayed. This figure also shows the shapes of low-velocity anomalies observed in this area by Di Stefano et al. (2009) and Molinari et al. (2015). Also in this map, a clear correlation between Pg-waves travel time residuals and the crustal structure is clearly noticeable. In fact, Pg residuals at the stations located to the NNW of the epicenter are systematically larger than those at the stations located to the SSE, denoting lower velocities to the former stations with respect to the latter.

Fig. 7
figure 7

Map of the stations receiving Pg-waves as first arrivals from the hypocenters, with the contours of crustal low Vs structures, redrawn from Molinari et al. 2015. The colors of the stations symbols represent the observed average Pg first arrivals residuals. The red cross indicates the 29 July 2019 epicenter

Full size image

We may conclude that the higher accuracy of the hypocenter locations obtained by both the relative location and the JHD methods, with respect to those obtained by the single-location method, allows us to state that the sources of the four considered events are separated by distances of the order of one km both in horizontal coordinates and in depths. This is consistent with the size of a tectonic structure generating a seismic sequence whose mainshock had magnitude Ml = 3.7, with three aftershocks exceeding Ml 3.0.

The analysis of travel time residuals has shown a good correlation with the crustal structure of the studied region, both for Pg waves, at a local scale in the range of 0–150 km distance from the epicenter, and for Pn waves, at epicentral distances ranging from 150 to 400 km for the Central Europe-Alpine region considered in this study.