1 Introduction

This geophysical investigation addresses the contact area between the Adriatic microplate and European plate in the northern Dinarides, which belongs to a very complex tectonic system in the southern part of the European lithosphere. No high-resolution three-dimensional model of the Earth’s crust and uppermost mantle is available yet for the northern Dinarides and adjacent southwestern Pannonian basin, and this is the motivation for undertaking this research. We used the Local Earthquake Tomography (LET) method to contribute to a better understanding of the crustal structure and its relationship to the upper mantle.

During the last 15 years, several geophysical experiments significantly contributed to the present understanding of the relationship between the Pannonian and Dinaridic crusts. In the present-day tectonic map (Fig. 1) after Schmid et al. (2008), the Dinaridic crust is built up by thrust slices derived from the Adriatic microplate, while the crust of the Pannonian basin comprises the Sava suture zone and tectonic units originally derived from the European plate (Tisza Mega-Unit).

Fig. 1
figure 1

Tectonic map of the wider area, according to Schmid et al. (2008). The superimposed boundaries of the Pannonian crust, Transition Zone, Dinaridic crust, and Moho fragmentation are based on gravity modelling (Šumanovac 2010)

The geophysical efforts to reveal the structure of the Earth´s crust in the Alpine–Carpathian–Pannonian–Dinaridic region were intensified by several international projects such as CELEBRATION 2000 and ALP 2002 (Guterch et al. 2003; Brückl et al. 2003). More recent investigations by teleseismic tomography also greatly contributed to the clarification of the lithospheric structure in the Dinaridic region (Šumanovac et al. 2017), especially in the northern Dinarides where a shallow fast velocity anomaly was discovered by Šumanovac and Dudjak (2016). A dense distribution of temporary and permanent seismic stations in the northern Dinarides allowed us to establish a more accurate tomographic model based on the local earthquake tomography in comparison to past seismic tomography investigations. The previous investigations were based on teleseismic tomography (e.g., Piromallo and Morelli 2003; Dando et al. 2011; Šumanovac et al. 2017), a method that is not convenient for resolving crustal structures due to the subvertical ray paths.

The primary objective of this investigation is the establishment of a 3D P-wave velocity model by means of local earthquake tomography. The 3D velocity model, used in conjunction with information on the crustal structure from previous studies (Brückl et al. 2007; Šumanovac et al. 2009, 2016; Šumanovac 2010), allows for better constraining the relationships between the crust and uppermost mantle in the survey area. Local events during two time periods were selected, and two databases were created to achieve a satisfactory spatial sampling. We used (1) data from temporary seismic stations within the ALPASS-DIPS project (Šumanovac et al. 2016) and (2) data from permanent seismic stations placed in Croatia, Slovenia, Bosnia and Herzegovina and Hungary.

A geological–geophysical model for a section across the northern Dinarides named Alp07 profile was constructed based on various geophysical methods (Šumanovac et al. 2016). The use of several geophysical methods diminishes an interpretation ambiguity and makes this profile a good reference for comparison with our model and to validate the results. The comparison with seismic refraction modelling and modelling based on receiver function (RF) analysis along the Alp07 profile allows for substantiating our velocity model, which is relevant for constraining the crustal and uppermost mantle structure. First-order anomalies such as the Moho boundary and the high-density body in a zone of transition between the crust of the External Dinarides and the crust beneath the Pannonian basin (Šumanovac 2010) will be resolved and compared with previous studies.

2 Study area and geological setting

The study area is located on the territories of southern Slovenia, central–northern Croatia, northwestern Bosnia and Herzegovina and southwestern Hungary. The region encompasses the contact zone between the Dinarides and the Pannonian basin (Figs. 1, 2). The strike of the Dinarides is NW–SE. The Dinarides formed after the collision of the Adriatic microplate with Europe-derived units such as the Tisza Mega-Unit in the Late Cretaceous and Early Paleogene (Tari and Pamić 1998; Ustaszewski et al. 2010), namely, during the SW verging folding and thrusting in the Late Eocene to Early Oligocene affecting the Adriatic margin (Korbar 2009).

Fig. 2
figure 2

Topographic map of the investigated Pannonian–Dinaridic region and histogram of hypocentre depths (below). Yellow dashed line delimits the area of the solution model where data coverage is sufficient. Bd Bjelovar depression, Kd Karlovac depression, Pd Požega depression, Mg Moslavačka gora Mountain, Sg Samoborsko Gorje Mountains

In the map, the Dinaridic mountain belt, which can be divided into the External and Internal Dinarides, links up with the Alps in the northwest. The External Dinarides are mostly composed of very thick thrust sheets, in places containing Mesozoic carbonate deposits with a stratigraphic thickness of more than 8000 m (Tišljar et al. 2002). The major part of the carbonate succession of the External Dinarides, namely, the Dalmatian and High Karst Units, is composed of deposits of the Adriatic Carbonate Platform (Vlahović et al. 2005). Two main thrust faults are commonly traced in tectonic maps of the survey area: the Ćićarija fault at the front of the External Dinarides (Fig. 1) and the more internal thrust system of the Velebit Mountains. Many authors trace a thrust fault in the base of the southwestern slope of the Mt. Velebit, but Tomljenović et al. (2017) questioned this fault. They argued that there is no evidence of its existence neither on the surface nor in the seabed.

The more internal parts of the NW Dinarides consist, going from SW, of the Pre-Karst and Bosnian Flysch Unit, non- to low-grade metamorphic units derived from a more distal part of the Adriatic, and the Western Vardar Ophiolitic Unit and the Sava zone (Pamić et al. 1998; Schmid et al. 2008). The Western Vardar Ophiolitic Unit and the Sava zone as defined by Schmid et al. (2008) are also known as the Dinaridic Ophiolite Zone and the Sava–Vardar zone, respectively, as defined by Pamić (2002). The Sava–Vardar zone represents a belt of ophiolitic, magmatic and metamorphic rocks and marks the Cretaceous–Palaeogene suture between the Internal Dinarides and the Tisza block (Schmid et al. 2008; Ustaszewski et al. 2009). All ophiolitic occurrences west of this suture belong to an ophiolite sheet of the Dinarides obducted in Jurassic times (Western Vardar Ophiolitic Unit, also known as Dinaridic Ophiolite Zone).

The eastern edge of the study area includes the western part of the Tisza block, mostly located in the subsurface of the Pannonian basin. Most authors agree that Middle Miocene crustal thinning and rifting in the Pannonian basin occurred as a result of back-arc extension in the context of the Carpathian subduction (e.g., Royden et al. 1983). An additional concept regarding the deep crustal configuration in the Pannonian basin was proposed by Dando et al. (2011). The authors suggested that higher-velocity mantle material became detached (or delaminated) from the crust beneath the Pannonian basin. Mantle downwelling and the detachment of higher-velocity material produced by previous subduction was considered to have triggered extension in the Pannonian basin.

The main depressions within the investigated part of the basin, separated by relatively shallow basement highs, are four WNW–ESE trending elongated sub-basins (Fig. 2): the Drava depression, the Bjelovar depression that continues into the Požega depression, the Sava depression, and the Karlovac depression (Pavelić 2001). These sub-basins are filled by Miocene sediments, which are thickest in the Drava and Sava depression areas (Saftić et al. 2003).

There are also several inselbergs. The Slavonian Mountains and Moslavačka gora Mountain, between the Sava and the Drava depressions, are composed of Palaeozoic–Mesozoic rocks. The outcrops of ophiolites are found on the northern edge of the Sava depression in the Medvednica and Samoborsko gorje Mountains (Pamić 1993; Pamić and Tomljenović 1998). The Southern Marginal Fault of the Pannonian basin (SMF) is thought to represent the SW tectonic margin of the Pannonian basin.

The Mid-Hungarian Shear Zone represents an approximately 30 km-wide zone within the Pannonian basin, which contains slices of South Alpine and Dinaridic origin (Kovács and Haas 2010). This transpressional zone is taken as the boundary between the Tisza and Alcapa Mega-Units, along which these two terranes rotated in opposite senses in Late Oligocene to Early Miocene times and became juxtaposed with each other (Csontos and Nagymarosy 1998). The Alcapa Mega-Unit includes the Eastern Alps, the West Carpathians and the Transdanubian ranges north of Lake Balaton (Csontos and Vörös 2004).

3 Previous geophysical investigations

Deep seismic sounding experiments mapped Moho depths of 40–48 km beneath the Dinarides and 22 km beneath the Pannonian basin (Dragašević and Andrić 1968). Many publications on the Earth´s crust and lithosphere structure originate from international projects that covered the Alpine–Pannonian–Dinaridic region (Brückl et al. 2003, 2007). Along the Alp07 profile (Fig. 2), which stretches in the SW–NE direction, i.e., perpendicular to the Dinarides, Šumanovac et al. (2009) interpreted three types of crust, based on 2D seismic and gravity modelling: (1) a two-layered Dinaridic crust, (2) a one-layered Pannonian crust and (3) a wide Transition Zone in between. Šumanovac (2010) used the same seismic model to calibrate densities and improve the resolution of the 2D gravity modelling in the entire Dinarides and surrounding regions. A gravity model based on the calibrated density set enabled the construction of a lithosphere model in the contact zone of the Dinarides with the Pannonian basin (Šumanovac 2015).

Active-source seismic surveys were followed by the passive seismic experiment ALPASS (Mitterbauer et al. 2011), and the ALPASS-DIPS was a part of this project (Šumanovac et al. 2016). In the framework of this experiment, Šumanovac et al. (2016) identified major crustal discontinuities by a receiver function analysis and constructed a geological model along the Alp07 profile based on the joint interpretation of seismic modelling, gravity modelling and receiver function analysis. Receiver function analysis was also applied to determine the lithospheric structure on seismic stations in Istria (Orešković et al. 2011) and in the Pannonian basin (Hetényi and Bus 2007).

The investigated area has been studied by numerous experiments investigating deep lithospheric structures. Positive velocity anomalies under the mountain chains, representing dipping slabs, were inspected by seismic tomography models on regional and global scales (Koulakov et al. 2009; Piromallo and Morelli 2003). In the first experiments using a greater density of seismic station distribution, the Dinaridic area is only partially covered (Mitterbauer et al. 2011; Dando et al. 2011). A fast, shallow anomaly (to depths of approximately 200 km) was imaged only in the central and southern Dinarides by Piromallo and Morelli (2003) and Koulakov et al. (2009), whereas Šumanovac et al. (2017) reported the presence of the fast anomaly within the same area, which extends to depths of up to 400 km. A fast anomaly beneath the northern Dinarides, interpreted as a lithospheric slab, was discovered for the first time due to the dense distribution of the seismic stations in that area (Šumanovac and Dudjak 2016). In recent work covering the Dinaridic region, Belinic et al. (2018) used the S-wave receiver function method to investigate the lithospheric thickness distribution. Authors reported a lithospheric thickness between 100 and 120 km below the northwestern Dinarides, thinning towards the Adriatic Sea and Pannonian basin, and anomalously thin lithosphere in the central Dinaridic region.

4 Data and method

Two datasets collected by seismic stations placed in the survey area were used. The first set comprises seismic waveform data from 15 temporary short-period seismic stations which were deployed in the framework of the ALPASS-DIPS project (Šumanovac et al. 2016). The stations were placed along the ALP07 profile, two in central Istria and one in northern Croatia (Fig. 2). These temporary short-period stations continuously recorded during the period from November 2005 to May 2007 with a sampling frequency of 50 Hz. The local events in this period were selected from the International Seismological Centre (ISC) Bulletin (2018). The second seismological database comprises data recorded by the Croatian State Seismological Network stations, including its Zagreb-Net subnetwork (operating instruments owned by the City of Zagreb), the Seismic Network stations of the Republic of Slovenia, the Hungarian National Seismological Network stations and the Seismic Network of the Republic of Srpska. The second database was obtained from the Croatian Seismological Survey and the EIDA–European Integrated Data Archive infrastructure within ORFEUS (EIDA 2017). The local events for the second database were selected during the period from January 2007 to June 2016.

The hypocentre parameters were not included in the inversion procedure for the solution model. The earthquake hypocentres and origin times were taken from the Croatian Earthquake Catalogue (CEC), and the parameters were calculated by the latest version of the HYPOSEARCH program based on a grid-search algorithm (Herak 1989). The program uses P- and S-wave arrival times and the crustal model from Bureau Central International de Séismologie (1972), which was modified by adding a thin layer on top. The CEC acquired its latest format with its revised release for the period 1908–1992 (Herak et al. 1996). The catalogue is being supplemented on a yearly basis, and its latest version is stored in the archives of the Department of Geophysics of the Faculty of Science, University of Zagreb. All data were processed manually. To determine the earthquake parameters, data from the Croatian Seismological Network were supplemented by readings reported in monthly bulletins of seismological stations in the neighbouring countries. We concluded that the parameters from the CEC are more accurate than those obtained by the relocation of the earthquake hypocentres with our dataset, especially for the shallowest earthquakes. The author of the HYPOSEARCH program demonstrated that the errors of the origin time are better than ± 0.02 s, with standard deviations < 0.02 s (Herak 1989). Although thousands of earthquakes were recorded within the model boundaries, we experienced problems such as noisy seismograms and a lack of data for some stations during several periods. This drastically reduced the number of observations used in the analysis. The analysis included 2583 first-arrival times from 119 earthquakes with magnitudes mostly ranging from 3 to 4. The azimuthal gaps varied from 6° to 329°, and the average gap was 166°. Phases were determined based on Bureau Central International de Séismologie (1972) with the assessed reliability added. The picking uncertainty in both datasets did not exceed 0.35 s, and the average was 0.08 s. The deepest earthquakes used in the analysis occurred at a depth of 30 km (Fig. 2).

Local Earthquake Tomography is able to derive a 3D velocity model of the shallow lithosphere with a horizontal resolution up to a few tens of kilometres. The model depth is not limited by the range of focal depth of the earthquakes because the method uses not only first arrivals but also refractions and reflections. The Moho can be reached even in the case of shallow events. The Local Earthquake Tomography method denotes the process of iterative inversion for the 3D velocity structure and hypocentre parameters using travel times from local earthquakes. First-arrival travel times (Pg and Pn phases) were used since we could determine them with greater precision than the other phases.

The forward problem of travel time prediction is performed with an algorithm implementing a robust grid-based eikonal solver, known as the fast marching method or FMM (Sethian and Popovici 1999). The inversion process is solved with an iterative non-linear subspace inversion method (Kennett et al. 1988). The inverse problem in seismic tomography can be presented in terms of minimizing an objective function consisting of a data term and one or more regularization terms. Solving the inverse problem requires the adjustment of the model parameters to satisfy the data, and the process is subject to regularization constraints such as damping and smoothing parameters. The subspace algorithm projects the entire problem into a smaller n-dimensional subspace. Local linearity is assumed in an inversion step, but the use of an iterative approach means that the nonlinear relationship between the travel time and velocity is reconciled.

The method allows the simultaneous inversion of multiple parameter classes such as velocity parameters and interface depth parameters and considers different dependencies of the objective function on various parameter types (Rawlinson and Sambridge 2003). To choose the best trade-off between minimizing the data misfit and generating the smoothest model, inversion processes with varying damping and smoothing parameters were performed. For the solution model, we used damping and smoothing parameters of 2.5 for the velocity inversion and 5.0 for the interface inversion.

The velocity and Moho structure are represented with cubic B-spline functions that are controlled by the grid of velocity nodes and the grid of interface nodes, respectively. The velocity grid consists of 25 radial grid points and 35 nodes in the north–south and east–west directions. This means that the grid allows for a spatial resolution of about four kilometres in the vertical and approximately 14 kilometres in the horizontal direction. However, resolution tests showed that with our dataset, we could not achieve the best resolution allowed by the model parametrisation. The number of interface nodes in both the north–south and east–west directions was 40. The discrete sampling of the velocity grid to perform the forward step is represented by a propagation grid, and there are 60 propagation grid nodes spanning a depth range of the model, which results in a grid spacing of approximately 1.7 km. In the north–south and east–west directions, the propagation grid spacing is approximately 8 km.

The tomographic inverse problem is generally underdetermined due to a lack of data and errors in the data. When the ratio of unknowns to data exceeds 12:1, the problem is simultaneously over- and underdetermined in different parts of the model (Sambridge 1990). The variation of the ray path densities is large (Figs. 3, 4), and areas within the model that are underdetermined are not considered in the solution model. The excessive 3D variation of the data constraints is a result of the limited distribution of seismic stations in the area and is typical for almost all datasets of this type. When using a subspace inversion method, the solution model remains generally unperturbed from its original state in regions with poor data coverage. The subspace method is a gradient-based inversion procedure that will not adjust parameters that have zero Fréchet derivatives (Rawlinson et al. 2008). Therefore, the solution only has meaning in areas of good path coverage, which are marked in the solution figures.

Fig. 3
figure 3

Path coverage across the model produced by the inversion of the recorded data (a few events used are outside the boundaries of the map, namely in the territories of Austria, Hungary, Serbia, Montenegro and Albania)

Fig. 4
figure 4

Path coverage across the solution model plotted in a vertical cross-section along the ALP07 profile

The inversion of the P-wave travel times was performed in six iterations, which resulted in a 95% reduction in the data variance, from 8.707 s2 (Root mean square error, RMS = 2.948 s) to 0.372 s2 (RMS = 0.552 s). The normalized χ2 value for the solution model is > 1, meaning that it does not fully satisfy the data, but this is to be expected when dealing with real applications (Rawlinson et al. 2010). The ray coverage does not support resolving small-scale anomalies. As pointed out by Diehl et al. (2009), this is indicated by the fact that the final RMS is still significantly larger than the average picking error. Regardless, the final data fit suggests that the recovered lateral heterogeneity is required by the data. With this dataset, we sought a smoother model, which means that we have less to interpret but the data has increased confidence.

Since the global velocity model ak135 (Kennett et al. 1995) is inappropriate as an initial reference velocity model in a local study, we used a multi-layered model, which is more appropriate for the area of research. However, it is important to highlight that by using ak135 and a constant gradient model as initial velocity models, the main features and their relative relationships were very similar in inverse models and in our final model. This indicates that the velocity perturbations in the solution model do not depend significantly on the initial velocity model. The consistency between anomalies common to all inversions suggests that the relatively high and low velocities actually represent the data.

5 Tomographic results

5.1 Resolution test

Although there have recently been discussions about the potential drawbacks of the synthetic reconstruction test (Rawlinson et al. 2014; Rawlinson and Spakman 2016), the spatial resolution of the tomographic solutions is still commonly evaluated with synthetic models. The synthetic test allows reaching conclusions about the solution robustness, which depends on the path coverage and data quality. The most widely used synthetic test is the so-called checkerboard test, in which the synthetic model is represented by an alternating high and low seismic velocity pattern. A uniform distribution and sufficient number of seismic propagation paths within the whole investigated model is hardly possible in real applications. Thus, in synthetic models, the forward modelling is solved using the same path coverage and phase types as in the observed dataset. The inversion method to the synthetic dataset is then applied to recover a synthetic model. The resolution is estimated from the capability of recovering positive and negative velocity anomalies from the synthetic model.

We constructed a synthetic model of the positive and negative velocity anomalies with a rectangular shape and peak amplitudes of dV = ± 0.8 km/s. The dimensions of the anomalies were 80 × 30 × 8 km. The space between every pair of anomalies, that is, a block with dV = 0, allows obtaining a clearer picture of how anomalous structures may be smeared out along ray paths. Travel times were calculated for the same source–receiver ray paths as in the observed dataset, and the inversion is controlled by model parameters very similar to those in the solution model. Random noise was not added to the synthetic travel times.

The synthetic and recovered models are shown by several horizontal (Fig. 5) and vertical (Fig. 6) cross-sections. The Alp07 profile and the profile stretching at 45.3°N are presented within the vertical sections. The quality of the recovered checkerboard pattern is generally good down to a depth of approximately 50 km in the northern Dinaridic area and in part of the Pannonian basin. The resolution is progressively lost below a depth of ~ 50 km. Most earthquakes in the survey area are located in the upper crust.

Fig. 5
figure 5

Three-dimensional resolution test using the checkerboard test method (triangles represent seismic stations). Left: depth slices of the synthetic model. Right: inversion slices for the same depth. A checkerboard pattern with positive and negative velocity anomalies is superimposed on the initial velocity model, which is used for the final model. The locations of profiles that are presented in Fig. 6 are shown on the 18 km depth slice of the synthetic model. Areas of poor resolution are blanked

Fig. 6
figure 6

Checkerboard recovery test results for two vertical cross-sections (for positions on Fig. 5). The black line within the model represents Moho. Areas of poor resolution are blanked

The full amplitudes of velocity anomalies are usually not recovered, but the general pattern of anomalies is recovered in the area 44.5–46.5°N and 14–18°E. This was expected since it is the part of the model with the best ray coverage. We conclude that the applied geometry can reveal potential structures of several tens of kilometres in latitude and longitude, with only minor smearing. The parts of the inverse model showing major smearing are blanked out in the figures.

The smearing in the slice at an 18 km depth is a result of the position between positive and negative anomalies (see vertical cross-section of Fig. 6). The structures above and below are vertically smeared. Vertical and horizontal smearing is present even with uniform distributions of rays and stations, as shown by synthetic tests (Rawlinson et al. 2014; Rawlinson and Spakman 2016). With synthetic reconstruction tests, Rawlinson et al. (2014) showed that despite the significant increase in the number of rays, the recovered model is not improved (see their Fig. 3). First arrivals are attracted to higher-velocity regions, and rays from close sources bunch together and do not contribute to better constraining a structure. Bezada et al. (2016) performed an experiment applying the real distribution of events and rays and an ideal (regular) distribution, which is impossible to obtain in practice. They showed that, for isotropic inversion, the solutions are very similar for the two distributions.

The smearing of structures can be more accurately estimated with discrete spikes in the sensitivity test. Figure 7 shows a cross-section through synthetic and recovered models comprising two anomalies with peak amplitudes of dV = ± 0.6 km/s. Gaussian noise with a standard deviation of 0.1 s is added to the synthetic travel times. The dimensions of the spikes are close to those of the main anomalies interpreted in the final model. In general, the input low- and high-velocity anomalies are recovered, but there are varying levels of distortion. The recovered structures seem to have greater dimensions than those in the input model. There is a tendency to smear structures vertically, which should be considered when interpreting a solution model. The presence of low background noise can be ascribed to the Gaussian noise added.

Fig. 7
figure 7

Synthetic test—input and recovered models of discrete spikes. Dashed contours outline the input anomalies

According to the results of the synthetic tests and ray coverage, a relatively good resolution is expected for the central part of the Alp07 and 45.3°N profiles down to 50 km.

A gentle curvature of the boundary representing the Mohorovičić discontinuity was added in the synthetic model (Fig. 6), but the used inversion parameters could not recover the curvature. However, this did not affect the velocity recovery, which is good since we infer the location of the Moho discontinuity from the high-velocity gradient in the lower crust and uppermost mantle. Overall, the tests show that the used dataset is acceptable for detecting a realistic range of crustal velocity anomalies of several tens of kilometres in horizontal directions. We can use the results to detect the crustal structure in the northern Dinarides area and a part of the Pannonian basin referred to as the Transition Zone (Šumanovac et al. 2009) into the Pannonian basin. Tomographic images are usually most easily interpretable in vertical cross-sections. The interpretable area of our solution model is delimited by the boundaries shown in Fig. 5 (dashed line on the horizontal depth slices).

5.2 Tomographic images

The solved P-wave 3D velocity model is presented in several horizontal (Fig. 8) and two vertical cross-sections (Fig. 9). The P-wave velocity distribution shows that the crust under the northern Dinarides is characterized by relatively stronger lateral and vertical velocity changes when compared to the crust in the Pannonian basin area. There is a shallow near-surface velocity anomaly under the northern Dinarides. A velocity greater than 6 km/s is noticed within the first 10 km of depth (HV in Fig. 9—on the part of the profile marked as belonging to the Dinaridic crust and HV anomalies at 8 km in Figs. 8, 11). Within the first few kilometres of depth, this high-velocity anomaly is continuous in the northern Dinaridic area. However, this anomaly seems to be separated by a low-velocity zone stretching perpendicular to the Dinarides from the Istria area across the narrow part of the northern Dinarides in the area of stations labelled 1–7, which is visible in the 8 km depth slice (Fig. 8). The velocity at a depth range of 15–25 km in the north Dinaridic region (Fig. 9) is relatively low, but there is a sharp transition towards higher velocities at lower crustal levels, which suggests that the crust could be two-layered. In the lower crust, the most noticeable feature is a large low-velocity region below the Dinarides, due to the changes in the Moho depth (LV in Fig. 8—depth slices 30 and 38 km and “thickening” in Fig. 9).

Fig. 8
figure 8

Depth slices across the P-wave inversion model. HDB high-density body from Šumanovac (2010), marked by a solid black line; HV high-velocity anomaly, LV low-velocity anomaly. Triangles denote seismic stations. Positions of two vertical cross sections are shown in the 20 km depth slice. All anomalies marked here are superimposed on the tectonic map in Fig. 11

Fig. 9
figure 9

Alp07 section and west–east section along 45.3 N across the P-wave inversion model. The location of the profiles is shown in Fig. 8 (depth slice at 20 km). The solution is blanked out where coverage is poor. HV high-velocity anomalies, KD Karlovac depression, SD Sava depression. Triangles on Alp 07 profile denote seismic stations. Red arrows indicate areas of relative crustal thickening and thinning. Black dashed lines represent an intra-crustal discontinuity and the Moho discontinuity

The Pannonian basin includes several depressions. Within the first few kilometres of depth, there are prominent low-velocity anomalies underneath the Sava depression (5.0–5.5 km/s), and the Karlovac depression (≈ 5.5 km/s). Our P-wave velocity model identifies the locations of these basin depressions and their relative depths, but the depressions reach greater depths than was expected. We infer that the velocities in the upper crust of the Pannonian segment are very influenced by depressions whose depths exceed 5 km in some parts (Saftić et al. 2003). Most of the hypocentres are in the upper crust, so a significant portion of the shallow seismic ray paths traverse the Karlovac, Sava, and Drava depressions. The low seismic velocities in the sedimentary rocks are probably smeared to a greater depth, as indicated by the resolution tests. We tried to include depression corrections by correcting the observed arrival times for the seismic stations placed in a depression within the Pannonian basin. Considering the depth of a depression (Saftić et al. 2003) and the lower velocity, the observed time was corrected for the rays arriving at these stations. For simplicity in the calculation, we assumed that all rays that arrive at a station have the same path length through the depression, although this is not strictly true because the depressions are asymmetric and elongated. This means that some rays have a longer paths across a depression than others and that the corrections cannot be the same for all rays at one station. The depth of the low-velocity anomalies in the Pannonian segment did slightly decrease after including the corrections, but the main features and relations in the model remained the same. We assume that this is a good indicator of inversion stability.

The smearing of low velocities partially “hides” high-velocity areas beneath the Sava depression. A high-velocity anomaly under the Sava depression is visible in the depth range 5–15 km (HV on Fig. 9—on the part of the profile marked as belonging to the Transition Zone). In the Transition Zone between the northern Dinarides and Pannonian basin, a high-velocity anomaly was determined by 2D seismic modelling, and a high-density body (HDB—marked on the 8-km depth slice in Fig. 8) was determined by 2D gravity modelling (Šumanovac et al. 2009; Šumanovac 2010). These features were interpreted as being related to a volume of oceanic crust beneath the Sava depression. In our model, the HV anomaly seems to be shifted to the southwest with respect to the HDB, and it partly overlaps with an anomaly from a 2D model (dashed line on the same depth slice). The high-velocity anomaly marked in the Transition zone of the 45.3° N vertical cross-section (Figs. 8, 9—15 km depth-slice at 16° longitude) looks like it was split into two parts because it is interrupted by a low-velocity area, probably caused by the effect of the Sava depression on the seismic velocities. The Sava depression is deepest just in the section at ≈ 16.5° longitude where this anomaly is discontinued in the vertical cross-section. In a tectonic framework, this high-velocity anomaly is within the ophiolitic zone (Fig. 11).

The most obvious feature in the model concerns the variations in the Moho depth. The criterion for the identification of the LET Moho discontinuity is a strong velocity gradient across the discontinuity. Diehl et al. (2009) showed that Moho can be resolved even with a rather coarse model parametrization (with a vertical extension of the model cells of 15 km). Depending on the vertical cell dimension, the velocity gradient across the Moho can be weak, which results in an atypical velocity value for the identification of the Moho in an LET model. Considering that the strongest velocity-depth gradient from a typical crustal velocity to a typical mantle velocity in our model is within four km, we can say that the Moho interface uncertainty is also within four km. There is a good correlation between the topography of the Moho defined along the Alp07 profile (Šumanovac et al. 2016) and the tomographic iso-velocity surface of Vp ≈ 7.0 km/s. The depth range of typical crustal velocities under the Dinarides indicates a greater Moho depth under the Dinarides, and a high velocity (> 7.0 km/s) at small depths in the Pannonian segment point to crustal thinning. These features are very well presented in the depth slices of 30 and 38 km. Additionally, the low-velocity area has the NW–SE oriented Dinaridic trend (Fig. 8). The greatest Moho depth is under the Dinarides, reaching approximately 45 km, and the shallowest (25–30 km) is in the Pannonian basin area (Fig. 9).

6 Discussion

We compared the results obtained from our LET 3D-velocity model with previous geophysical experiments in the investigated area. The investigated area is partly covered by the ALP01, ALP02 (see Fig. 2 from Brückl et al. 2007) and ALP07 profiles within the ALP 2002 project based on active-source seismic refraction experiments (Brückl et al. 2007; Šumanovac et al. 2009). We put an emphasis on the correlation with the 2D model along the Alp07 profile based on seismic refraction data, gravity modelling and receiver functions (Šumanovac et al. 2009, 2016).

A comparison of the main features shows the agreement of our results with those obtained by other methods, even though a direct comparison between different classes of seismic images is difficult and often needs calibration. Active-source seismic imaging better reveals discontinuities and faults, while tomography produces smooth variations in the seismic properties. Although converting the velocity field information into interface information is much easier for the Moho interface than for crustal discontinuities, the LET method can still obtain reliable results if quality data are used.

The Istria area has poor ray coverage from the NW, W, and SW directions, and it is at the very edge of the interpretable area in the LET model, as seen from the resolution test. It is also an area with poorer ray coverage in the ALP01 model (Brückl et al. 2007), and there are apparent differences between these two models. In the LET model, the upper crustal velocity is approximately 5.7 km/s, while in the ALP01 profile, the Adriatic foreland is characterized by velocities higher than 6 km/s (Brückl et al. 2007). Furthermore, the pronounced shallowing of Moho depth (from 40 to 28 km) at the southern end of the Alp01 profile by Brückl et al. (2007) was not obtained in the LET model. Velocities in the lower crust are similar in these two models (≈ 6.4 km/s). A better correlation can be noticed in the models along the Alp02 profile (Brückl et al. 2007), located in the Pannonian basin. In both models, the uppermost crust is characterized by velocities of approximately 5.5–5.9 km/s, and in the depth range of 10–30 km, the velocities are between 5.8 to 6.5 km/s. The Moho depth in the Alp02 model varies from 27 to 29 km in the Pannonian part, which is shallower by up to 4 km than the Moho depth according to our model.

The most detailed comparison of our LET data with those from controlled-source experiments can be made along profile ALP07 (Šumanovac et al. 2009). The main characteristics of the two models are shown in Fig. 10. The shape of the Moho discontinuity in the LET model is similar to that in the 2D model based on seismic and gravity modelling and receiver function analysis (Šumanovac et al. 2016). However, differences in the depth reach up to 5 km. In the LET model, the lateral and vertical velocity variations are more pronounced in the Dinaridic crust than in the Pannonian part, and it is inferred that the Dinaridic crust is two-layered. The intra-crustal discontinuity postulated on the basis of the LET 3D-model is a few km deeper but in general parallel to the one inferred from the 2D model by Šumanovac et al. (2009).

Fig. 10
figure 10

Alp07 cross-section through the P-wave inversion model with superimposed features defined in the 2D model. The topography is shown on the top, and the triangles denote seismic stations. HV high velocity, LV low velocity. The solid blue and red lines are crustal and Moho discontinuities inferred from the velocity model established by controlled-source experiments along the Alp07 profile of the ALP 2002 experiment (Šumanovac et al. 2009). The grey solid line represents the Moho inferred from receiver function modelling (Šumanovac et al. 2016). The dashed blue and red lines are discontinuities interpreted on the basis of our LET model. The dashed black line denotes the estimated mean depth for the Moho obtained by combining all three models. Note the vertical exaggeration. Numbers in green colour denote Vp seismic velocities from the 2D seismic refraction model, while HDB and crustal faults are represented by solid green and red dotted lines, retrospectively (Šumanovac et al. 2009)

It is important to emphasize that the lateral changes of low and high velocities in both, the LET and the Alp07 model (Šumanovac et al. 2009), show a similar pattern. In the Pannonian basin, the low velocities in the uppermost few kilometres are partly affected by the velocity distortions caused by the depressions of the bedrock below the Miocene basins. We suspect that the relatively deep depressions cause vertical smearing of the low-velocity anomalies in the LET model. The smearing probably dampens the high-velocity anomaly beneath the Sava depression, where high-velocity and high-density anomalies were obtained in previous investigations (Šumanovac et al. 2009; Šumanovac 2010). The high-density body shown by Šumanovac (2010) seems to be displaced to the northeast relative to the HV anomaly in the LET model (Fig. 11). The high-density body appears in a block on the gravity profiles Alp07, GP-1 and GP-2 (Šumanovac 2010) and is thickest on the ALP07 profile (up to 15 km). The thickness of the HDB decreases suddenly to approximately 8 km on the other two profiles, GP-1 and GP-2, SE from Alp07 (see Figs. 4, 5 from Šumanovac 2010). The thinning of the SE part of the HDB body and vertical smearing of the low velocities beneath the Sava depression cause the LET velocity model to be unable to reveal that part of the HV body. Elongated and asymmetric depressions probably distort the ray travel times. This can cause an apparent displacement of the HV body with respect to the HDB body.

Fig. 11
figure 11

Overview of the anomalies in the survey area superimposed on the tectonic map. Boundaries of the major tectonic units are after Schmid et al. (2008) and Šumanovac (2010). MF Moho fragmentation

We decided to use a relatively smooth structure for the Mohorovičić discontinuity, even though our tomography model indicates that the Moho flanks below the Dinarides could be steep and with a sudden increase in depth on both sides of the Dinarides (Fig. 10). The results of the resolution test indicate that we cannot discuss in detail the slope of the Moho flanks, but a steep north-eastern flank, i.e., asymmetrical flanks for the deepest Moho below the Dinarides, point to the over-thrusting of the Pannonian segment. The under-thrusting of the Adriatic microplate below the Pannonian segment is commonly accepted in the literature. A gravity study assumed Moho fragmentation based only on structure geometry below the Dinarides (Šumanovac 2010). We have generally limited the interpretation to a depth of 50 km, and the highest ray coverage in the LET model is obtained in the uppermost mantle of the central part of the Alp07 profile (Fig. 4). Therefore, we can attempt to explain the HV and LV features beneath the Dinarides. In the depth range of 40–55 km, the low-velocity anomaly is surrounded by high-velocity anomalies (Fig. 10). These anomalies in the 52 km depth-slice (Fig. 8) are superimposed on the tectonic map (Fig. 11). The low-velocity anomaly in the upper mantle corresponds well with the Moho fragmentation, as determined by the geometrical relationships and the asymmetry of the flanks at the Moho below the Dinarides (Šumanovac 2010). The low velocity could be caused by the fragmentation in the uppermost mantle, pointing to the contact of the Pannonian and Adriatic mantles due to the sinking of the Adriatic slab beneath the Dinarides, which was recently discovered by teleseismic tomography (Šumanovac and Dudjak 2016; Šumanovac et al. 2017). According to the velocity distribution (Figs. 10, 11), the subduction zone can be located on the NE flank of the Dinarides.

Overall, the LET velocity model with its main velocity anomalies shows a good correlation with previous models (Šumanovac et al. 2009, 2016; Šumanovac 2010) in the contact zone between the northern Dinarides and the southwestern margin of the Pannonian basin. We think that one can talk with great certainty about crustal thickening below the Dinarides and thinning under the Pannonian basin. The intra-crustal discontinuity is slightly different in different models, but it is certain that a two-layered structure is only found in the Dinaridic crust. We expect that the application of this method for establishing a 3D velocity model when applied to the entire Dinaridic region and the Pannonian marginal zone will give a more precise velocity distribution, and calibration with previous models will make geological interpretation easier.

7 Conclusions

A new 3D velocity model for the crust and uppermost mantle of the northern Dinarides and southwestern Pannonian basin was obtained by the tomographic inversion of source–receiver travel times. The results of the LET tomographic experiment presented in this paper provide evidence that significant changes in the crustal structure occur in the contact area between the northern Dinarides and southwestern Pannonian basin. We only interpreted relatively large structures in the survey area, according to the estimated resolution, on the basis of the applied resolution tests. However, this 3D model is correlated and calibrated with previous 1D and 2D geophysical data, thereby improving the effective resolution. The LET velocity model will significantly diminish the ambiguities in the geological interpretation within the survey area.

The velocity maps show good agreement with the most prominent geological features. The data allow for the detection of low-velocity and high-velocity bodies within the crust, pointing to the main geological structures. The high mountain chain in the northern Dinarides is associated with relatively high velocities at shallow depths (< 10 km), but in the depth range of 10–20 km, lower crustal velocities are inferred. There is a high-velocity anomaly obtained in the depth range of 5–15 km in the area between the Karlovac and Sava depressions (belonging to the Internal Dinarides). This anomaly is placed somewhat deeper and shifted to the southwest in comparison to the high-density body defined by Šumanovac (2010) and could be interpreted as a block consisting of oceanic rocks.

High vertical velocity gradients occur in the crust beneath the Dinarides at a depth of approximately 20 km, which suggests that the Dinaridic crust could be two-layered. The intra-crustal discontinuity seems to be somewhat deeper but is in general parallel to the one detected in the 2D seismic model by Šumanovac et al. (2009). In the Pannonian basin area, the high-velocity gradient in the uppermost crust can be correlated with deep sedimentary basin fills in the Sava and Drava depressions.

A deep low-velocity zone beneath the Dinarides at depths of 30–40 km is interpreted as crustal thickening. The model is characterized by relatively high velocities at these depths in the Pannonian basin area, pointing to crustal thinning. The characteristic NW–SE trending of the crustal thickening is visible on the horizontal cross sections, and the Moho shape can be determined in vertical cross sections on the basis of the highest vertical velocity gradient. The estimation of the Moho shape and depth is based on the iso-velocity surface Vp ≈ 7.0 km/s. It generally matches well with the existing 2D geological–geophysical model in the survey area. The low-velocity anomaly beneath the Dinarides in the uppermost mantle could be interpreted as being a result of Moho fragmentation, localized at the contact between the Pannonian and Adriatic mantles.