Deep Structure and Dynamics of the Central Balkan Peninsula from Seismic Data

Abstract—Analysis of P- and S-receiver functions for 19 seismic stations on the Balkan Peninsula has been performed. Half of the stations are in Bulgaria. The crustal thickness varies from 28–30 to 50 km. The ratio of longitudinal and shear wave velocities in the upper crust reaches 2.0 in some places. In the southwest of the study area, the 410-km seismic boundary is uplifted by 10 km relative to nominal depth. The elevation may be caused by hydration and/or cooling of the mantle transition zone under the influence of the Hellenic subduction zone. A low S-wave velocity layer related to the 410-km boundary may be located atop this boundary. In the northwestern part of the study area this layer is present in spite of the absence of the 410-km boundary. A similar paradox has been previously noted in central Anatolia. Indications of a low-velocity layer are present at a depth exceeding 410 km. The simultaneous inversion of the receiver functions of the two types (P and S) and the Rayleigh wave phase velocities reveals a large (7–9%) decrease in the S-wave velocity in the upper mantle of southern Bulgaria and northern Greece. The thickness of the low-velocity layer (asthenosphere) is about 50 km. The lithosphere-asthenosphere boundary (LAB) is at depths of 40 to 60 km. In terms of tectonics, this zone is characterized as the South Balkan extension system. To the north of 43° N, the S-wave velocity in the upper mantle is usually at least 4.4 km/s and the LAB is not detected or is detected at a depth of over 80 km. The SKS analysis of azimuthal anisotropy reveals lateral zoning in the upper mantle that is correlated to velocity zoning. Probably, the mechanically weak low-velocity mantle of the South Balkan system is easily deformed, and the azimuth of the fast direction of anisotropy (20°) indicates the direction of extension. At the northern stations, the fast direction (about –30°) may be a reflection of an older process.


INTRODUCTION
Since the Late Jurassic, the tectonic processes in the Balkan Peninsula ( Fig. 1) have been determined by the collision of the Arab, Asian and European plates. Between the main plates, there are fragments of continents formed by the rifting of Pangea in the early Mesozoic. The Alpine orogeny resulted from the collision of one of these fragments (Apulian) with the southeastern margin of Europe (Dinter, 1994). In the Balkan Peninsula, the ophiolitic Vardar zone is considered the locus of the elements of the Alpine Fault. To the west and southwest of the Vardar zone, there are deformed elements of the Apulian fragment-the Dinaric and Hellenic Alps. They run parallel with the eastern coast of the Adriatic Sea and constitute the marginal part of the Hellenic subduction system, which was active throughout most of the Neogene. Between the Carpatho-Balkan and Dinaro-Hellenic branches of the Alpine orogeny, there are tectonic zones of the Rhodopes, Sredna Gora, Stara Planina, Pre-Balkan Thrust, and the Moesian Platform.
In the north, the Rhodope Massif borders the Sredna Gora area, which is part of a large Late Cretaceous magmatic belt formed over the north-dipping Neo-Tethys Plate. During the Late Cretaceous, the area was undergoing an extension phase. Later on, this zone was involved in the Alpine compression phase along with the Balkan folded zone to the north. In the south of the Balkan Peninsula, there is the South Balkan extension system (Burchfiel et al., 2008), which is located north of the North Anatolian Fault Zone and south of the poorly defined tectonic boundary in northern Bulgaria. The crust in this area underwent thickening, heating, and a decrease in viscosity. Extension prevailed in the area from the Late Paleogene through the Holocene. Today, this area is separated from the more active extension zone in the IZVESTIYA Aegean region south of the North Anatolian Fault Zone. The deep-seated structure of the crust and mantle of the Balkan Peninsula is generally known from geophysical research data. The deep seismic sounding has shown that the crustal thickness ranges between 30 km for the Moesian Platform and 50 km for the Rhodope Massif (Volvovsky et al., 1985). Estimates of velocity inhomogeneities for P-waves were obtained to a depth of about 250 km (Botev and Spassov, 1990). The analysis of dispersion curves of surface waves provided estimates of the crustal thickness and S-wave velocity in the upper mantle to a depth of about 300 km (Raykova and Panza, 2015). The joint analysis of seismic and gravity data was used to evaluate the topography of the 410-km boundary (Yegorova et al., 1998), which represents the top of the mantle transition zone. This work estimates that the 410-km density boundary beneath the Sredna Gora zone is 50 km higher than the normal depth (410 km) in adjacent areas. The structure of the crust and the transition zone of the mantle were considered in the articles (Georgieva, 2015;Georgieva and Nikolova, 2013).
Our research is aimed at a more detailed analysis of inhomogeneities of the crust, upper mantle, and mantle transition zone of the Balkan Peninsula based on seismic data. The main tool of our analysis is the receiver function technique, which allows us to explore seismic boundaries and analyze the deepseated structure at a higher resolution than was previously possible. The most common version of this technique uses P-wave receiver functions (PRF). The main difficulty in using PRF is that transmitted Ps-waves from boundaries in the upper mantle are recorded almost simultaneously with strong multiples in the crust, which act as noise. This problem is essentially solved by using S-wave receiver functions (SRF) (Farra and Vinnik, 2000), but the SRF technique has its own challenges. Among the transmitted Sp-waves in the SRF technique, the wave from the Moho boundary has the largest amplitude. As a result of long-period filtering, this waveform acquires a side lobe that resembles in time and amplitude the Sp-wave from the lithosphere-asthenosphere boundary (LAB) and causes interpretation errors. Both problems (multiple reflections and side lobes) are automatically solved by a simultaneous inversion of PRF and SRF. Therefore, the method version we apply uses Pand Swave receiver functions together. In addition, we use Rayleigh waves and analyze the azimuthal anisotropy of the upper mantle based on SKSand SKKS-wave recordings.

RECEIVER FUNCTIONS: TECHNIQUES AND SEISMIC DATA
In our work we used the recordings of 19 broadband seismic stations (Fig. 1). The stations belong to different networks and have been in operation for several years. About half of the stations are in Bulgaria. The epicenters of most teleseismic earthquakes used here are located east of the Balkan Peninsula (Fig. 2). For each seismic station, we built between 45 and 90 receiver functions of longitudinal waves (PRF) and between 25 and 50 receiver functions of shear waves (SRF) at epicentral distances of 30° to 90° for PRF and 65° to 95° for SRF.
PRF is built in the LQ coordinate system. The L-axis is parallel to the main direction of displacements in the P-wave in the wave propagation plane. The Q-axis is perpendicular to the L-axis in the same plane. To suppress noise, the recordings are filtered with a low-pass filter with an angular period of about 5 s. The waveforms of different earthquakes are standardized by deconvolution in the time domain (Berkhout, 1977), and useful signals are selected by applying migration to the standardized components. The move-out time corrections for the migration are calculated to detect converted waves from the boundaries in the depth range from 0 to 800 km.
When SRF is built, the Q-axis is parallel to the main direction of motion in the S-wave in the wave propagation plane. The L-axis is perpendicular to the Q-axis in the same plane. The low-pass filter has an angular period of 8 s. Similarly to PRF, the waveforms are standardized by deconvolution in the time domain and detected by migration. The move-out time corrections for the migration are calculated as the product of the deviation of the Sp-wave slowness from the Swave slowness (differential slowness) and the deviation of the epicentral distance from the reference distance (differential distance).
The results of building PRF are shown in Fig. 3. The seismic stations are divided into 4 groups: 1-GRG, KKB, KNT, SKO, VAY; 2-ALN, EDRB, RZN, PLD; 3-DJES, MPE + PLVB, TRAN, VTS; 4-TIRR, PRD, PSN, PVL. Using groups instead of single stations improves the signal-to-noise ratio. The migration in Fig. 3 allows us to observe the distinct arrivals of the P410s and P660s waves. The standard error of the P410s and P660s time estimates in the groups evaluated by the Bootstrap technique (Efron and Tibshirani, 1986) is about 0.2 s. The standard arrival times of these waves for the IASP91 model (Kennett and Engdahl, 1991) at a slowness of 6.4 s/deg or an epicentral distance of 67° for the surface source are 44.0 and 67.9 s, respectively. The differential time (the interval between the arrivals of P660s and P410s waves) is 23.9 s. The wave arrival times in the IASP91 model serve as a good approximation for data from most continental stations .
The group 1 times are 43.0 and 68.0 s for P410s and P660s, respectively. The differential time is 25.0 s, which is 1.1 s longer than the standard time. These data mean that the 410-km boundary is elevated by about 10 km, while the 660-km boundary is almost at standard depth. As a result, the differential time is increased by about 1 s, and the P410s time is reduced by the same value. The differential time for the second group is 23.9 s, which equals exactly the standard value. The P410s and P660s seismic phases start at times of 45.0 and 68.9 s, i.e., with a delay of 1 s relative to standard values. The delay is due to the low wave velocity at depths of under 410 km.
The P410s seismic phase in the second group data is preceded by a phase with approximately the same amplitude and opposite polarity. This phase is formed at the upper boundary of the low-velocity (partially melted) layer located above the 410-km seismic boundary. In the third group, the P410s wave is not detected due to signal weakness, but the wave from the top of the low-velocity layer above the 410-km boundary at a time of 39.0 s is clearly identified. The depth of this boundary is about 360 km. A very similar wave field (a wave with negative polarity and no P410s) is observed at the stations in central Anatolia (Vinnik et al., 2014). The formation of a low-velocity layer above the 410-km boundary is explained by dehydration during the upward flow of the mantle substance and the olivine-wadsleyite phase transition during the crossing of the 410-km boundary (Bercovici and Karato, 2003). The unusual wave field in the third group indicates a complex process not only above the 410-km boundary, but also below it. Finally, group 4 has the P410s and P660s phases with a differ- The results of SRF migration for the four groups are shown in Fig. 4. The established sign rule requires negative polarity of the Sp wave formed at the boundary between media with a low S-wave velocity above the boundary and a high velocity below the boundary. The wave with negative polarity at a time of -3 s is formed at the Moho boundary. The phase preceding this wave, with positive polarity and lower amplitude, represents the side lobe of the Moho wave. At a time of about -53 s, the receiver functions detect a wave with negative polarity and a maximum amplitude of about 0.03 with a slowness of 0.4-0.6 s/deg. This is an Sp-wave formed at the 410-km boundary. The pulse with negative polarity is preceded in groups 2, 3 and 4 by a wave with a comparable amplitude and positive polarity. We interpret this peculiarity as an Sp-wave formed at the bottom of the low-velocity layer in the depth range of 450 to 510 km. The large amplitude, comparable to the amplitude of the S410p, does not allow us to interpret this wave as a side lobe of the S410p seismic phase. The indications of this layer were previously found in several other regions (e.g., (Vinnik et al., 2012)). In our case, the anomaly of the mantle transition zone is under the Black Sea. The theoretical S410p phase time for the IASP91 model in the four groups is -52.3, -52.9, -52.3, and -52.3 s, respectively. The actual wave arrives 0.8, 0.5, 0.8, and 1.5 s earlier. The early arrival can be caused by an increased ratio of longitudinal to shear wave velocities or by a depressed by several kilometers 410-km discontinuity (Vinnik et al., 2010).

SIMULTANEOUS INVERSION OF PRF, SRF,
AND P-AND S-WAVE TRAVELTIME ANOMALIES We use the previously described technique for a simultaneous inversion (e.g., (Vinnik et al., 2007)). We assume that in the vicinity of the seismic station, the medium under study can be represented by a laterally homogeneous, isotropic stack of layers. The lateral inhomogeneity can be described as a mosaic of homogeneous blocks. The search for optimal models is performed by means of an iterative procedure similar to the simulated annealing procedure (Mosegaard and Vestergaard, 1991). Trial synthetic PRF and SRF for building the target function are computed by the Thomson-Haskell method (Haskell, 1962) with flattening (Biswas, 1972). The stack typically consists of 9 layers, each described by three free parameters: longitudinal wave velocity (Vp), shear wave velocity (Vs) and thickness. The density is expressed through Vp according to Birch's law. We consider several randomly chosen initial models, for each of them an iterative sequence of 10 5 trial models is generated. We assume that the mantle velocity inhomogeneities are at depths not exceeding 300 km. Therefore, the Vp and Vs values at a depth of 300 km and deeper are fixed at the parameters of the standard IASP91 model (Kennett and Engdahl, 1991). The last 5% of the models are used to estimate the posterior distribution of the parameters. We break the parameter space of the model into cells and represent the results of inversion by the number of hits in each cell.
The receiver functions are sensitive to changes in velocity with depth and, to a lesser extent, to absolute velocity values. In order to increase the sensitivity to absolute values, it is reasonable to apply a simultaneous inversion of the receiver functions and independently measured deviations of Pand S-wave traveltimes of teleseismic earthquakes from the standard values. The simplest way to estimate anomalies is based on the knowledge that in many regions of the world the main seismic boundaries in the transition zone (410-km and 660-km) are at almost the same depth. Significant lateral inhomogeneity of these boundaries is known mainly in subduction zones and hotspots. The greatest lateral inhomogeneity of the upper mantle is usually found at depths of under 300 km. In this case, dTps-deviation of the arrival time of the Ps-wave formed at the boundary in the transition zone from the standard value can be expressed as dTps = dTs -dTp, where dTs and dTp are the teleseismic anomalies of Sand P-waves. This equation can be rewritten as where K-the ratio of S and P wave traveltime anomalies. Estimates of the K value have been repeatedly discussed in the literature (e.g., ). Most published K values range from 3.0 to 4.0. A significant, about 1 s, deviation of the differential time from the standard value (23.9) indicates lateral inhomogeneity of the transition zone. In this case, estimates of teleseismic anomalies of Pand S-wave traveltimes can be obtained otherwise, from observations of surface Rayleigh waves.
We obtained estimates of traveltime anomalies for our network using data on the phase velocity of the fundamental mode of the Rayleigh wave in the range of periods from 35 to 150 s (Ritzwoller et al., 2002). A similar analysis of observations in the area of the central Anatolian plateau was performed earlier (Vinnik et al., 2014). The phase velocity values were obtained in the nodes of the coordinate grid with an interval of 2°. The shear wave velocity values as a function of depth were obtained for the same nodes by the method (Hermann and Ammon, 2002). The obtained velocity sections have good sensitivity to phase velocity at depths of up to 200 km, i.e., in the depth range of the lithosphere and asthenosphere. The teleseismic anomalies dTs (Fig. 5) are determined by subtracting the time for the IASP91 model. The anomaly dTp is determined from the equation dTs/dTp = K.
An example sequential inversion of the receiver functions of the two types and the dTs and dTp travel- time anomalies is shown in Fig. 6 for the ALN station. Figure 6a shows the result of inversion without dTs and dTp traveltime anomalies. In this case, the section for S-waves shows the LAB boundary at a depth of 60 km with an S-wave velocity jump of about 0.4 km/s (9%). The lower boundary of the low-velocity layer (asthenosphere) is at a depth of 100 km. The Moho boundary is at a depth of 30 km. The ratio of Vp and Vs velocities in the upper 10-km layer of the crust is close to 1.9. The width of the green corridor for Vp and Vs is more than 1 km/s. When dTs and dTp are used (Figs. 6b, 6c), the width of the green corridor for Vs is reduced by several times, while for Vp it remains practically unchanged. The main peculiarities of the Vs section, which include the depth to the main boundaries and the Vp/Vs ratio in the upper crust, are preserved. The brown corridor becomes more distinctive. Changing K = dTs/dTp from 3.0 (Fig. 6b) to 4.0 (Fig. 6c) has hardly any effect on the inversion result. Further on we consider the Vs results at K = 3 in Fig. 7 to be the most reliable. The upper mantle region covered by the seismic data is displaced about 100 km northeast of the station. The SKO station is an exception, because in this case the seismic sources are predominantly located west of the station. The first group of obtained models includes Vs sections for the GRG, KKB, KNT, SKO, VAY stations (Fig. 7a) in the southwestern part of the study area. The Moho boundary for all the sections of this group is at a depth of 28-30 km.
The LAB is at a depth of 40, 55, 60, 55, and 40 km, respectively. The S-wave velocity in the asthenosphere is 4.0-4.2 km/s, which is 7-9% lower than the IASP91 standard value. The second group consists of sections for the ALN, EDRB, RZN, PLD stations

AZIMUTHAL ANISOTROPY
We supplemented the study of mantle inhomogeneity in the Balkan Peninsula with an analysis of azimuthal anisotropy, partly using the same seismic network. Because of the poor quality of the recordings, several stations have been replaced by new ones. Table 1 shows the parameters of the recordings used and the results obtained, which are also shown in Fig. 8. An example of the recording analysis is shown in Fig. 9. In the case of azimuthal anisotropy, the shear wave splits into two orthogonally polarized quasi-shear waves. The polarization of quasi-shear waves is determined by the elastic moduli of the medium in which the shear wave propagates. The main mineral that influences the shear wave splitting in the upper mantle is olivine. SKS and SKKS wave recordings are used to determine  (Ritzwoller et al., 2002     α and δt, the polarization direction of the fast split wave and the delay time of the slow split wave (Vinnik et al., 1984). The obtained estimates of azimuthal anisotropy characterize the area within a radius of several tens of kilometers from the station.
Our estimates are obtained by the method (Vinnik et al., 1989) assuming a transverse-isotropic model with a horizontal axis of symmetry. The anisotropy parameters are determined by minimizing the target function representing the mean square of the difference between the observed SKS or SKKS T-component and the theoretical T-component synthesized from the observed radial (R) component. The optimal α and δt parameters are determined by exhaustively searching through all possible values in the nodes of the dense grid and for all SKS and SKKS recordings with acceptably low noise levels. The standard error of estimating fast direction from a single good quality recording is about 10°, but we usually use several recordings of the same station and get a smaller error.
The analysis results (Fig. 8)    in this region should be easily deformed and the azimuth of the fast wave polarization can indicate the direction of extension. A fast direction of about -30°i s associated with the high-velocity zone.

DISCUSSION AND CONCLUSIONS
Our analysis reveals a ~10 km upward flow of the 410 km seismic boundary in the Rhodope region. The elevation could have occurred as a result of temperature decrease or hydration (Karato, 2011) of the mantle transition zone. The temperature decrease required to elevate the boundary by 10 km is about 100°С. The Hellenic subduction zone extending along the Adriatic coast is oriented favorably for moving the cold material of the lithosphere 200-300 km to the east, where the temperature decrease in the transition zone is seen in the seismic data.
The second noteworthy observation is related to the low-velocity layer immediately above the 410-km boundary. The origin of this layer is associated with the elevation of the mantle substance, during which the olivine-wadsleyite phase transition takes place, accompanied by the release of water and partial melting (Bercovici and Karato, 2003). The normal structure of the transition is demonstrated by the data of the group of ALN, EDRB, RZN, PLD stations located southeast of the study area (Fig. 7b). In this case, the P410s seismic phase with positive polarity is preceded by the same phase, but with negative polarity. This is the wave from the upper boundary of the partially molten layer. There is an unusual situation with the phase transition in the area of the DJES, MPE+PLVB, TRAN, VTS stations. In this case, we clearly observe the phase with negative polarity, but there is no normal P410s wave with positive polarity, although both phases are related and should exist simultaneously. A similar anomaly was found at the stations of central Anatolia (Vinnik et al., 2014).
Another noteworthy observation is made in the analysis of the S410p wave recordings. In this case, we observe the S410p precursor wave with opposite polarity (Fig. 4). To explain this precursor, it should be assumed that there is a low-velocity layer at depths from 450 to 510 km, Keshav et al. (2011) reported a sharp temperature drop in the solidus of the carbonate mantle at pressures typical of this layer. The 520-km boundary, the origin of which is debatable, may correspond to the base of this layer (Vinnik et al., 2012).
The simultaneous inversion of the receiver functions of the two types (P and S) and the Rayleigh wave phase velocities reveals a sharp (7-9%) decrease in the S-wave velocity in the upper mantle of southern Bulgaria and northern Greece. Based on the tectonic features, the low-velocity region is characterized as an extension system (Burchfiel et al., 2008), and the decrease in the shear wave velocity may be caused by partial melting during decompression. The thickness of the partially molten layer (asthenosphere) is about 50 km. The lithosphere-asthenosphere boundary (LAB) is found at depths of 40 to 60 km. In other words, partial melting of the mantle can occur near the Mohorovičić discontinuity. To the north of 43° N, the S-wave velocity in the upper mantle is usually at least 4.4 km/s and the LAB is not detected or is detected much deeper than in the south. Specifically, the S wave velocity decrease characteristic of the asthenosphere is not found at the DJES, TIRR, and VTS stations or is distinguished at a depth of 85 km (TRAN station) and at a depth of 200 km (VRI station).
Our SKS analysis of azimuthal anisotropy reveals zoning in the upper mantle that is similar to velocity zoning. In the southern group of stations, the fast wave polarization azimuth α is about 20°. At most of the other stations, the azimuth α is about -30°. The majority of the southern group of stations is located in the region that is characterized as an extension system based on the tectonic features. The upper mantle of this region is in a partially melting state. The fast wave polarization azimuth in this region probably characterizes the main direction of extension. The second group of stations is located in the Moesian platform influence zone with a thick lithosphere. A similar orientation of the anisotropy was previously found in the neighboring regions of the East European Platform (Dricker et al., 1999). To understand the dynamics of   these regions better, we should study anisotropy as a function of depth.