Interpretation of the tectonic evolution of the western part of the Sava Depression: structural analysis of seismic attributes and subsurface structural modeling

This paper presents the results of an extensive structural investigation of the western part of the Sava Depression (SW part of the Pannonian Basin System) and provides insights in the tectonic evolution of the investigated area. Structural analyses were based on the 3D seismic volume of the study area, analysis of seismic attributes, and construction of 3D structural model and structural maps. Our results pinpoint to three tectonic phases in the structural development of the Sava Depression. The first tectonic phase is characterized by extensional tectonic features developed from the Early to Middle Miocene. The second tectonic phase follows thermal subsidence and general deepening of the study area and inherited tectonic features during the Late Miocene, while the final tectonic phase is characterized by structural reactivation and tectonic inversion of inherited and newly formed tectonic features from the Pliocene to the Quaternary. ARTICLE HISTORY Received 7 January 2019 Revised 24 August 2019 Accepted 29 August 2019


Introduction
The study area, western part of the Sava Depression (WSD) is located in the central part of NW Croatia, in the vicinity of the town Zagreb ( Figure 1). It covers an area between Medvednica Mt. in the NW, and Moslavačka Gora Mt. in the SE, while the SW and W boundaries are formed by the NE banks of the Sava River ( Figure 2). Previous (throughout the twentieth century) and recent geological investigations in the NW Sava Depression could be clustered into surface and subsurface geological studies. Surface investigations are generally presented by the sheets of the Basic Geological Map at 1:100,000 scale (Šikić, Basch, & Šimunić, 1977(Šikić, Basch, & Šimunić, , 1979Basch, 1980;Pikija, 1987aPikija, , 1987b. Subsurface research incorporated numerous studies by the Croatian Petroleum Company INA d.d. that comprised data collected from many exploration wells and acquired geophysical datasets.
The aim of this study is to improve current knowledge about the local features of the tectonic evolution of the SW Pannonian basin, based on the results of structural analysis and subsurface modeling of the WSD. The results are presented as a Neogene and Quaternary tectonic reconstruction of the tectonic evolution of the study area. Reconstruction encompasses structural analysis and interpretation of geophysical and exploration well datasets, provided by the Croatian Petroleum Company INA d.d. Furthermore, seismic attribute analysis was performed, with subsurface mapping of three stratigraphic horizons and corresponding structure contour maps.
Formation of the Croatian part of the Pannonian Basin System (PBS; Figure 1) is initiated during Early Miocene, originating in the area of Hrvatsko zagorje and the western part of the Drava Depression, with the ENE-WSW directed extension conditioned by a regional stress regime with an ENE striking P-axis (Tomljenović & Csontos, 2001). This stress field allowed the formation of a system of NNW striking normal faults that enabled the initial ENE directed extension and deepening of Hrvatsko zagorje and the western part of the Drava Depression that also continued during the Ottnangian and Karpatian towards the Sava, the Slavonija-Srijem Depressions and remaining parts of the Drava Depression (Pavelić, 2001;Lučić et al., 2001;Tomljenović & Csontos, 2001;Velić et al., 2002). Sedimentation of Miocene deposits commenced within half-grabens, formed in normal fault hanging walls, and it was characterized by heterogeneous freshwater sedimentation, that generally comprised coarse grained clastics and subordinately alteration of sandstones and siltstones (Pavelić & Kovačić, 1999;Saftić et al., 2003;Velić et al., 2002). Marine flooding of the Croatian part of PBS during Badenian (Ćorić et al., 2009), caused a transition from fluvial and lacustrine environments towards predominantly marine environments (Lučić et al., 2001;Pavelić & Kovačić, 2018;Pavelić, Miknić, & Sarkotić Šlat, 1998). Extensional phase and deposition of the syn-rift sediments can be correlated to the syn-rift tectonic phase (Pavelić, 2001) or to the first depositional megacycle (Velić et al., 2002) in basin evolution of Croatian part of PBS, associated to 300-2000 m thick rock succession in the Sava Depression (Saftić et al., 2003).
Termination of rifting processes appeared to be diachronous process in the entire PBS, however, it ended in the Middle Miocene (Balázs, Matenco, Magyar, Horváth, & Cloetingh, 2016). In the Middle Miocene regional stress regimes within the PBS significantly changed, being firstly manifested with regional thermal subsidence and deepening of the PBS basement and secondly characterized by local compressional events and partial basin inversion (i.e. Pogácsás, 1984;Saftić et al., 2003;Tari, 1994).
In the Croatian part of the SW Pannonian Basin (Figure 1), WNW and E directed Late Miocene thermal subsidence and deepening commenced by extensional reactivation of major inherited NNW striking normal faults (Tomljenović & Csontos, 2001). These faults, that represent the marginal 'opening' faults of existing  Schmid et al., 2008). Study area is shown with patterned polygon within the Croatian part of the Pannonian basin system. depressions enabled deposition of Late Miocene strata that range from 2000 to >4000 m in thickness at their deepest parts (Horváth & Royden, 1981;Pavelić, 2001;Saftić et al., 2003;Velić et al., 2002). The sediments were deposited within the Lake Pannon and encompassed clastic successions with a source area in   Rögl, 1996;Saftić et al., 2003).
Pliocene and Quaternary strata characterize sedimentation within the final tectonic phase in the evolution of the Croatian part of the SW Pannonian Basin. Characteristic formation of new NW-SE striking reverse faults, and reverse reactivation of former normal faults, accompanied tectonic inversion of the Croatian part of the SW Pannonian Basin. In addition, the whole PBS is conditioned by continuous convergence and CCW rotation of the Adria Microplate with respect to the European plate (D'Agostino et al., 2008;Jarosinski, Beekman, Matenco, & Cloetingh, 2011: Tomljenović & Csontos, 2001. Therefore, during the Pliocene and Quaternary, the horizontal stress field pattern (S Hmax ) in the PBS represents transferred far-field intra-plate stresses from collisional zones between the Adria indenter and stable Europe. In the Croatian part of SW Pannonian Basin principal horizontal stress axis aligns either N-S or NNW-SSE in the central part of NW Croatia, or NE-SW towards the E (Grenerczy, Fejes, & Kenyeres, 2002;Grenerczy, Kenyeres, & Fejes, 2000;Grenerczy, Sella, Stein, & Kenyeres, 2005). Consequently, because of contemporaneous tectonic uplift and subsidence of the deepest parts of the depressions, deposition occurred within remnants of Lake Pannon and marshes that altered into fluvial depositional systems during the Quaternary (Pavelić & Kovačić, 2018;Saftić et al., 2003 and references therein). With thicknesses between 500 and 1500 m drilled in the central part of the Sava Depression, the Pliocene and Quaternary strata can be correlated with the third depositional megacycle (Saftić et al., 2003;Velić et al., 2002).

Analysis of seismic attributes
In this study, initial determination and spatial interpretation of fault planes and fault zones within the study area was achieved by delineation of several seismic attributes. This analysis encompassed delineation of four seismic attributes: Chaos, Variance, Cosine of instantaneous phase and Sweetness.
The Chaos attribute (Randen et al., 2000) is conceived on analysis of dominating orientations (λ), where smooth reflector will have one dominating direction (λ max >> λ mid ≈ λ min ), bent reflector will have two strong directions (λ max ≈ λ mid >> λ min ) and 'damage' zone will have gradients pointing in all directions (λ max ≈ λ mid ≈ λ min ) (Randen, Pedersen, & Sønneland, 2001). Using this measure, regions with low consistency typically correspond to areas with chaotic textures. Usually, they present faulted zones, but also low consistency sedimentary textures, like carbonate reefs and channels (Marfo, Omosanya, Johansen, & Abrahamson, 2017;Randen et al., 2001;Sahran, 2017). On seismic data this attribute facilitated separation in horizontal sections of chaotic zones that hypothetically represent major structural discontinuities in the study area.
The seismic attribute of Variance (s 2 t ) is described (Azevedo, 2009) as edge imaging and detection method, which presents measure of how spread is the dataset around the mean value. The variance attribute uses an algorithm that computes the local variance of input seismic signal through a multitrace window with user-defined size: where x ij is the sample value at horizontal position, i and j is vertical time sample, w j-t is the vertical smoothing term over a window of length L. In seismic data studies, this seismic attribute is used for depicting depositional and structural features (Khair et al., 2012;Marfo et al., 2017). In our study, we used this attribute to delineate structural discontinuities. The Cosine of Instantaneous Phase seismic attribute is cosine of instantaneous phase angle (φ) and its measured in degrees and it is independent on the amplitude: where f(t) is real part of analytical seismic trace and f*(t) represents the imaginary component of the seismic trace, calculated through Hilbert Transform (Taner, Koehler, & Sheriff, 1979). Concept of use of this attribute is that amplitude peaks and troughs on seismic will remain at the same position but strong and weak events will exhibit equal strength, giving attribute capability to enhance reflector continuity and visual appearance of edges (Azevedo, 2009 (Hart, 2008) to detect sand channels in 3D seismic data. Attribute is derived by dividing reflection strength by the square root of instantaneous frequency (Radović & Oliveros, 1998). In our study Sweetness seismic attribute was conducted primarily in interpretation of vertical seismic reflection sections.
3.1.1. Subsurface modelling of the study area The 3D geometric properties of faults in the study area, their kinematics and timing of activity, during the Neogene and Quaternary, the values of relative displacements along stratigraphic horizons were all determined by construction of a 3D structural subsurface model. For this purpose, the 3D seismic volume of the western part of the Sava Depression with eight exploration well datasets were organized within the Petrel Seismic to Simulation Project database. A constructed 3D structural subsurface model integrated subsurface mapping of three stratigraphic horizons; which correlate to regional electro-log markers (Malvić & Saftić, 2008): Base Neogene -Tg and PTc markers, Base Pannonian -Rs7 marker and Base Plio-Quaternary deposits -α' marker and interpretation of 28 fault planes along inline and crosslines of the 3D seismic volume. Because of complex geological setting and unreliable time-depth data, presented structural maps (Maps 1-3 (Main map)) are constructed in the time domain (Two-way time, TWT, in milliseconds).

Early to middle Miocene extension within the WSD
Results of structural investigation and analysis of the seismic attributes and the constructed subsurface model indicate three distinct tectonic phases in the tectonic evolution of the western part of the Sava Depression during the Neogene and Quaternary. With the onset during the Early to Middle Miocene, the tectonic evolution of the WSD was characterized by extension and the formation of initial depocenters. Results show that the extensional phase and newly formed depocenters occurred along numerous mapped normal faults. Interpreted faults dominantly strike towards the NW (R-13, R-15, R-16, R-18, R-19, R-20, R-21, R-22, R-23), while subordinate faults strike towards the S (R-10, R-11, R-26) and NE (R-8) (Figures 5 and 6, Maps 1 and 2).
The intensity of fault's tectonic activity with extension and deepening is associated with significant vertical displacement values along the interpreted faults (Table 1). With values that occasionally offset pre-Neogene basement units around 150 ms (Table  1), normal faults are usually characterized by an enéchelon pattern accommodating predominant NE-SW directed extension (see faults R-13, R-18, R-19, R-20, R-21 and R-22 in SW part of study area; Figures  5 and 6, Map 1). Furthermore, within this study, several conjugate NW striking fault pairs (e.g. R-9 and R-15, R-11 and R-23; see Figures 4 and 5, Map 1) have been also distinguished, suggesting the possibility of local scale extension that could be associated with local trans-tension accommodated along the same faults. In the WSD, syn-tectonic deposits with thicknesses of c. 500 ms (Maps 1 and 2) are very often found as the first sequences in these normal faultbounded 'half-graben like' structures ( Figures 5 and  6). The same deposits represent eroded material of the uplifted footwall basement blocks that were locally transported and deposited in the hanging-wall of basement blocks. Beside the inherited palaeorelief (Map 1), the differential NE-SW extension with heterogenous uplift rates, subsidence and tilting of pre-Neogene basement blocks along NW striking faults produced a substantial deepening gradient that increases from the NW towards the SE. The range of TWT depth values on the structure contour map of the Base  Neogene are between 1000 and 2400 ms (Map 1) with two major depositional areas within the western part of the Sava Depression (Map 1).

Late Miocene thermal subsidence and deepening of the WSD
The second stage in the tectonic evolution of the western part of the Sava Depression was identified as Late Miocene thermal subsidence and deepening. Subsidence is evidenced by an excessive accumulation of a thick post-rift sedimentary succession of between 800 and 1300 ms of deposits (see Figures 5 and 6). With dominantly inherited structures, the structure contour map of the Base Pannonian surface (Map 2) of the WSD shows similar structural patterns with the NW striking normal faults depicted at the Base Neogene  Figure 6 and R-24 in Figure 5). These characteristics indicate that beside the dominant subsidence and deepening processes in the study area during the Upper Miocene, there are indicators of the onset of tectonic inversion because of change in the regional field of stress regimes (Saftić et al., 2003;Tari, 1994;Tomljenović & Csontos, 2001). The constructed Base Pannonian structure map indicates dominant sedimentary infill  Prelogović et al., 1998). B. Schematic NE trending geological cross-section within the western part of the Sava Depression. Note three major depositional megacycles on the geological cross-section (after Saftić et al., 2003).
from the NW, with subordinate contributions from the NE (Map 2).

Pliocene and quaternary tectonic inversion within the WSD (5,6 Ma-recent)
The youngest structural surface of the Basal Pliocene and Quaternary deposits were interpreted only in the S, SW, central and SE parts of the WSD. This was conditioned by the low quality of the available 3D seismic data in the NW part of the study area, where Plio-Quaternary strata in places outcrop at the surface. Hence, to overcome these difficulties, the structural surface of the Base Pliocene and Quaternary strata in the north-western part of was supplemented with the available structural data collected by Basch (1980) at the outcrops of the mapped geological units within the WSD (Figure 1). The resultant surface of the Base Plio-Quaternary strata has a different structural appearance from the previous two maps (Map 3). Results show that Pliocene and Quaternary structural reactivation and tectonic inversion within WSD is evidenced in the northern, central and eastern part of the study area. A predominant compressional stress regime is manifested as a series of reverse faulted structures associated with reactivated inherited NW striking normal faults (R-12 and R-23; Figures 5 and 6, Map 3) as well as to several newly formed NW striking reverse faults (R-5, R-6, R-15, R-24, R-37; Figures 5 and 6, Map 3). Besides reverse faulting, structural reactivation and tectonic inversion in the study area also occurred with folded Upper Miocene and Pliocene-Quaternary sedimentary successions. The best examples are uplifted structures within the areas of the NNW striking fault zones of the subvertical R-11 and R-23 faults, and R-15 and R-9 faults (Figures 5 and 6). In the Western and SW part, the intensity of the Pliocene and Quaternary tectonic inversion is less pronounced, seismic reflectors are sub-horizontal or gently dipping toward the SW, while identified NW striking faults (R-1, R-2, R-3, R-4, R-5, R-6; see Figures 5 and 6, Map 3) experience minimal vertical displacement (≤10 ms; see Table 1). Pliocene and Quaternary compression also resulted in the formation of an anticline structure, with NW-SE stretching positioned in the hanging wall of fault R-24. TWT values on the structure map are in the range of 0-925 ms indicating that the deepest parts of the WSD are in the southern and south-eastern parts. This is in general correlation with the previous locations of the local depocenters identified on the constructed structural contour maps.

Discussion and conclusions
The paper presents results of the structural investigation and analysis of the processed seismic attributes and 3D subsurface modeling of the western part of the Sava Depression. Complexity of structural setting of the study area prompted us to use seismic data processing trough seismic attributes that were tested in other studies in delineation of fault structures and other structural discontinuities, i.e. seismic reflectors (see references in Methods section). The principal application of Sweetness attribute was in identification and interpretation of stratigraphic boundaries. The Base of Neogene horizon interpretation was challenging as it presents palaeorelief surface with unclear contacts of basement rocks and overlaying basinal strata, however, Sweetness attribute enabled us to track change in coherency of seismic data. Chaos attribute facilitated good separation, in horizontal sections, of chaotic zones that on vertical seismic sections in Variance and Cosine of Instantaneous Phase attributes were associated to the structural discontinuities. Cosine of Instantaneous Phase also facilitated the lateral continuity of seismic reflectors and seismic facies in vertical sections. Concerning primarily structural investigation of seismic data, these four attributes derived by seismic data processing were sufficient to complete structural interpretation. The best results using seismic attributes were achieved in the deeper parts of the seismic reflection sections, where the amplitudes of seismic traces were more prominent. This is because of the target area i.e, hydrocarbon accumulations that are positioned in the deeper part of the WSD. Locations of interpreted stratigraphic horizons are set in respect to drilled contacts within exploration wells. In that sense, we followed seismic reflectors which correspond to picked well log horizon. Most of observed normal faults were formed and active during principal sin-rift extensional tectonic stage, however, particular normal faults continued extensive tectonic activity in the Middle Miocene. which beside Sava Depression correlate to other areas within Croatian part of Pannonian Basin System. Map of the Base of Neogene (Map 1) with number of normal faults shows palaeorelief in the base of syn-rift sediments or first depositional megacycle. Map of the Base of Pannonian (Map 2) convey inherited structures that delineate extensional tectonics during the rift phase, however, smaller number of normal faults and their vertical offsets indicate termination of syn-rift phase, which may be correlated to change from first to second depositional megacycle. Map of the Base of Plio-Quaternary (Map 3) shows completely different structural style, being characterized by compressional (reverse) tectonics, i.e. structural activity, which correlate to latest inversion and onset of third depositional megacycle.
Consistent with the tectonic evolution of the Croatian part of the Pannonian Basin, the first tectonic stage in western part of Sava depression is characterized by Early to Middle Miocene extension of the area with formation of the principal depocenters. Local extension was achieved along the en-échelon system of NW striking normal faults with general NE-SW directed extension ( Figure 5). During the syn-tectonic stage the most intensive deepening was identified in the SW, W and central parts with vertical displacement values up to 500 ms (Map 1). Similar depositional structures were documented in several portions of the Sava Depression and other Croatian parts of the PBS (see Cvetković et al., 2019 with references for details). In Sava Depression, extensional structures that are related to the opening of the Croatian part of PBS are described in studies of Tomljenović and Csontos (2001) and Cvetković (2013), whereas in Drava Depression similar extensional features are described in studies of Lučić et al. (2001) and Matoš (2014). At the regional scale, in the Hungarian part of the PBS Early to Middle Miocene extensional structures are described in detailed studies of Tari and Horvath (2006), Horváth et al. (2015), Balázs et al. (2016).
The second tectonic stage in the tectonic evolution began in the Early Pannonian and is characterized by further deepening, with deposition of more than 1300 ms of Pannonian and Pontian deposits, tectonic activity was allocated to the inherited NW striking normal faults. Overall, tectonic activity had decreased in intensity, however, it was mostly concentrated in the western and central parts of the study area (Map 2). This tectonic phase is often referred as Late Miocene extension (i.e. thermal subsidence phase) and it is recognized in other Croatian parts of PBS, e.g. in the Karlovac and Zagorje Sub-Depressions (see Tomljenović & Csontos, 2001), in the Drava and Mura Depressions (Lučić et al., 2001 andMatoš, 2014 with references) and Hungarian parts of the PBS (Balázs et al., 2016;Horváth et al., 2015;Tari & Horvath, 2006) During the Pliocene and Quaternary, the tectonic evolution of the western part of the Sava Depression was characterized by structural reactivation and tectonic inversion. A change in regional stress field, with an overall predominance of compressional and transpressional stress regimes yielded a structural reactivation of the former NW striking normal faults into reverse faults as well as the formation of a new set of NW striking reverse faults (Maps 1-3; Figure 5). The largest vertical offset of 200 ms is observed on reactivated normal fault ( Figure 5). The tectonic inversion is the most prominent within the NE and E parts of the study area, where beside reverse faulting and tectonic uplift along conjugate, almost subvertical pairs of NW striking reverse faults, a possible horizontal movement could be inferred (Maps 1-3). This could also indicate that beside dominant Pliocene and Quaternary compression and reverse faulting in the western part of the Sava Depression, a local transpression probably yielded a heterogeneous structural reactivation that change fault deformation characteristics along its strike. The Plio-Quaternary transpressionalcompressional tectonic phase is recognized at the numerous locations in PBS, however, structurally reactivated and tectonically inverted structures are especially pronounced along the margins of the PBS (e.g. Bada, Horváth, Gerner, & Fejes, 1999;Balázs et al., 2016;Horváth et al., 2015;Matenco & Radivojević, 2012;Prelogović et al., 1998;Tari & Horvath, 2006;Tari & Pamić, 1998).
Based on our results we conclude that study area presents a holotype of tectonic evolution for western portion of the Sava Depression with identified tectonic stages that are correlative to the tectonic evolution of the PBS. Though the study areais relatively restricted local scale feature, interpreted subsurface structures reflect regional scale tectonic events. In this study area structural reactivation and tectonic inversion shows different level of intensity that is expected taking into account position of surrounding inselbergs, fault framework, depression depth distribution and local stress field within the SW part of PBS (e.g. Matoš, Zajc, Kordić, Tomljenović, & Gosar, 2017;Prelogović et al., 1998;Rukavina, Matoš, Tomljenović, & Saftić, 2016;Tomljenović & Csontos, 2001 with references).

Software
Petrel is a software platform that gives access to various subsurface data through different types of visualization and analytical possibilities. Through this, the user can make technical and scientific results and interpretations which can provide new understanding about the subsurface architecture. For our study we only used several processes available within the Petrel E&P Software Platform. Mapping of this horizons together with faults was done using the Seismic Interpretation tool in the Geophysics process. Interpreted horizons and faults represented the input for structural modeling of our study area. Modelling was done with tools in the Corner point gridding process, with which we could model the geometry of fault structures in relation to interpreted horizons.

Disclosure statement
No potential conflict of interest was reported by the authors.