2D Numerical Simulation of Intraoceanic Subduction during the Upper Jurassic Closure of the Vardar Tethys

In this study we report interim results of our ongoing research that involves the application of numerical modeling for constraining the geodynamic conditions associated with the closure of the Vardar branch of the Tethys Ocean. The study is aimed at better understanding the ultimate fate of the Balkan ophiolites, namely at addressing the question whether these ophiolites represent relicts of an ocean that completely closed during Upper Jurassic/lowermost Cretaceous time (Vardar Tethys) or they also contain remnants of the ocean floor of a Late Cretaceous oceanic realm (Sava – Vardar) [Schmid et al., 2008].

In our numerical models we try to simulate a single intraoceanic subduction that commences in the Lower/Mid Jurassic and ends in the Lower Cretaceous, transitioning into oceanic closure processes and subsequent collision between Adria and Eurasia plates. These convergent-collision events should have led to the formation of ophiolite-like igneous rocks of the so-called Sava - Vardar zone.

A series of numerical simulations were performed with varying parameters. In the scope of our numerical simulations, the set of equations is solved: the continuity equation, the Navier-Stokes equations and the temperature equation. Marker in cell method was incorporated in solving this system with finite difference discretization of the equations on a staggered grid. To utilize this numerical method a thermo-mechanical code I2VIS [Gerya et al., 2000; Gerya & Yuen, 2003] was used for obtaining the final results. 

Our actual 2D thermo-mechanical models cover the crust and the upper portion of the mantle with varying starting widths of the Vardar Ocean in the Lower Jurassic. The ocean is modeled with two segments: the western subducting slab and the eastern overriding slab. Slabs with different ages and thicknesses were used and the convergence rate is varied. The intraoceanic subduction is assumed to have been initiated along the mid oceanic ridge. Two continents (i.e. Adria and Eurasia) with different thicknesses of the continental lithosphere and crust are also modeled adjacent to a single oceanic realm between them.

The parameter study is in function of defining conditions under which the hypothesized scenario occurs. So far, we have succeeded in reproducing westward obduction onto the Adriatic margin, which is in accordance with the geological observations, i.e., with the top-west emplaced West Vardar ophiolites [see Schmid et al., 2008 for references]. However, our model is yet to produce sufficient amounts of back-arc extension along the Eurasian active margin and that is crucial for explaining the formation of the igneous provinces occurring within the Late Cretaceous Sava – Vardar zone and the Timok Magmatic Complex.