Numerical simulation of the Late Jurassic closure of the Vardar Tethys

We report updated results of our ongoing research on constraining geodynamic conditions associated with the final closure of the Vardar branch of the Tethys Ocean by means of application of numerical simulations (previous interim results reported in EGU2020-5919).

The aim of our numerical study is to test the hypothesis that a single eastward subduction in the Jurassic is a valid explanation for the occurrence of three major, presently observed geological entities that are left behind after the closure of the Vardar Tethys. These include: ophiolite-like igneous rocks of the Sava-Vardar zone and presumably subduction related Timok Magmatic Complex, both Late Cretaceous in age as well as Jurassic ophiolites obducted onto the Adriatic margin. In our simulations we initiate an intraoceanic subduction in the Early/Middle Jurassic, which eventually transitions into an oceanic closure and subsequent continental collision processes.

In the scope of our study numerical simulations are performed by solving a set of partial differential equations: the continuity equation, the Navier-Stokes equations and the temperature equation. To this end we used I2VIS thermo-mechanical code which utilizes marker in cell approach with finite difference discretization of equations on a staggered grid [Gerya et al., 2000; Gerya&Yuen, 2003].

The 2D model consists of two continental plates separated by two oceanic slabs connected at a mid-oceanic ridge. Intraoceanic subduction is initiated along the ridge by assigning a weak zone beneath the ridge. Time-dependent boundary conditions for velocity are imposed on the simulation in order to model a transient spreading period. The change of sign in plate velocities is found to be useful for both obtaining obduction / ophiolite emplacement [Duretz et al., 2016] and causing back-arc extension. Changes in velocities are linear in time. Simulations follow a three-phase evolution of velocity boundary conditions consisting of two convergent phases separated by a single divergent phase where spreading regime is dominant. Effect of duration and magnitude of the second phase on model evolution is also explored.

Our so far obtained simulations were able to reproduce the westward obduction and certain extension processes along the active (European) margin, which match the existing geological relationships. However, the simulations involve an unreasonably short geodynamic event (cca 15-20 My) and we are working on solving this problem with new simulations.