ON THE SUBDUCTION OF THE APULIAN LITHOSPHERIC MICROPLATE UNDER THE EURO-ASIAN ONE AND THE MANTLE WEDGE THERMAL CONVECTION AS A POSSIBLE MECHANISM OF HYDROCARBONS UPWARD TRANSPORT IN THE PANNONIA AND THE VARDAR BASINS Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

GEO-INFORMATION SYSTEMS

UDC 530.311

Gavrilov S.V., Kharitonov A.L. On the subduction of the Apulian lithospheric microplate under the Euro-Asian one and the mantle wedge thermal convection as a possible mechanism of hydrocarbons upward transport in the Pannonia and the

Vardar basins

Gavrilov Sergei Vladilenovich

Doctor of physical and mathematical sciences, Main scientist of the laboratory 102, Schmidt Institute of Physics of the Earth

of the Russian Academy of Sciences Kharitonov Andrey Leonidovich, Candidate of physical and mathematical sciences, Leading scientist of the Main magnetic field laboratory, Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radio Waves Propagation of the Russian Academy of Sciences

Abstract. For the non-Newtonian mantle rheology case the 2D thermal viscous dissipation-driven thermal convection in the mantle wedge above the Apulian lithospheric microplate subducting under the Euro-Asian plate is modeled numerically. The effects of the 410 km and 660 km phase transitions are taken into account. Within the framework of the model constructed the horizontal extent of the 2D heat flux anomaly observed in the rear of the Dinarides mountain belt corresponds to subduction velocity of ~10 mm per year which is close to that observed with the help of geodetic means. In the case of non-Newtonian rheology the upwelling convective flow transporting heat to the Earth's surface locates at the distance from the trench corresponding to the actually observed 2D heat flux anomaly, the velocity in the convective vortices being from ~10 mm per year to ~10 m per year for the water content in the mantle wedge from 0.3X10-1 to 3X10-1 weight percent respectively. The convection cell dimension is of the order of the horizontal scale of the heat flux anomaly observed in the Pannonia basin and Vardar zone. Upwelling mantle wedge convective flow is indicated to be able to provide the mantle wedge hydrocarbons transport to the Earth's surface for the mantle wedge mantle content over 0.3X10-1 weight percent.

Keywords: mantle wedge 2D thermal convection, non-organic mantle hydrocarbons transport, subduction angle and velocity, mantle rheology, phase transitions in the mantle.

Рецензент: Сагитов Рамиль Фаргатович, кандидат технических наук, доцент, заместитель директора по научной работе в ООО «Научно-исследовательский и проектный институт экологических проблем», г.

Оренбург

Introduction

According to [3], the subduction of the Apulian lithospheric microplate under the Dinarides, Pannonia basin and Vardar zone is sufficiently flat and during the last ~ 45 Ma occurred at the angle of

~ 25°, which remained unchanged during this time interval. The genesis of the Dinarides mountain belt (with the transversal horizontal extent of ~ 300 km) apparently is of the thrust and fold nature, associated with the former collision and subduction of the oceanic branches of the Neo-Tethis and Alpine-Tethis as the result of thrusting of the African plate under the Eastern and Western Europe during the last 55 - 35 Ma [3]. In [1] numerous papers are referred, containing contradictory estimates of the relative motions of the Apulian lithospheric microplate and Euro-Asian plates, made on the basis of seismic, geophysical and geodetic data. In fig. 3 in [Op. cit.] the velocity of subduction of the Apulian lithospheric microplate under the Euro-Asian one is seen to amount to ~5 - 8 mm^a'1 according to referred works, while in [1] this velocity is estimated to be of the order of ~5 mm*a-1 according to calculations based on geodetic observations. In [12] the Pannonia basin and Vardar zone are noted to be the zones of the Middle Miocene extension occurred ~14 - 11.6 Ma ago, which led to the lithosphere thinning, these zones being the back-arc basin characterized by the back-arc spreading. At that time the single mountain belt parallel to Apulian shore was split into Carpathians and Dinarides and the shallow Pannonian sea was formed, which existed approximately to 600 thousand years ago. Presently the Pannonian oil- and gas-bearing basin is situated in this region. Here the conditions are clarified under which the centre of the back-arc spreading initiates as the result of convective instability, driven by the dissipative heat release in the mantle wedge above the subducting Apulian lithospheric microplate.

According to [4, 8, 9] two types of dissipation-driven small-scale thermal convection in the mantle wedge are possible, viz. the 3D finger-like convective jets, raising to volcanic chain, and 2D transversal Karig vortices, aligned perpendicularly to subduction. These two types of convection are shown to be spatially separated due to the pressure and temperature dependence of mantle effective viscosity, the Karig vortices, if any of them formed, being located behind the volcanic arc [4]. There are contradictory judgments on the velocity of subduction of the Apulian lithospheric microplate under the Euro-Asian one, although the order of magnitude of the present-day subduction velocity (~ 10 mm*«-1) can apparently be regarded as established sufficiently reliably. The mountainous massif Dinarides locates parallel to the north-eastern shore of the Adriatic sea, and probably is of the thrust-and-fold nature. The 2D maximum of the heat flux anomaly of ~100 mWxm-2 observed in the rear of the Dinarides massif in the Pannonia basin and the Vardar zone [10] can be assumed to owe its origin to the convective heat supply from the mantle wedge. Numerical modeling of the 2D convection, occurring in the mantle wedge in the form of the Karig vortices and presumably transporting heat upwards, allows to judge about the mean water content in the mantle wedge and to assume the mantle hydrocarbons to be transported to the Earth's surface by the upwelling convective flows. Numerical convection models accounting for the pressure-, temperature- and stress-dependence of viscosity fit best to observational data in the case of non-Newtonian rheology at the mantle water content of ~ (0.3 - 3)x10-1 weight % for the velocity of

subduction of ~ 10 mmxa"1 during the Middle Miocene. In [15] such rather a high water content (and even greater one, up to 3 weight %) can be observed in the mantle wedge in the mantle transition zone. The Middle Miocene subduction velocity of ~10 mmxa-1 during the formation of the Pannonia basin is of the order of the observed presently, or, can somewhat exceed it because of the gradual diminution of the velocity of convergence of African and Euro-Asian plates.

Model

Thermo-mechanical model of the mantle wedge between the base of the overlying Euro-Asian plate and the upper surface of the Apulian lithospheric microplate subducting under the Euro-Asian one with a velocity V at an angle P is obtained for the infinite Prandtl number fluid as a solution of non-dimensional 2D hydrodynamic equations in the Boussinesq approximation

(d2zz-d2xx)x ^x(d2zz-d2xx)x ^+4x52xz^x52xz w=Ra* Tx-Ra(410)xrx(410)-Ra(660)xrx(660),

(1)

dt T = AT- (w z x Tx) + (Wx x Tz) + (Di / Ra) x (x2lk / 2xrj) + Q ,

(2)

for the stream-function w and temperature T. Here r is dynamic viscosity, d and indices denote partial derivatives with respect to coordinates x (horizontal), z (vertical) and time t, A is the Laplace operator, r(410) and r(660) are volume ratios of the heavy phase at the 410 km and 660 km phase boundaries, the velocity components Vx and Vz are expressed through w as

Vx = W z , Vz = - W x

(3)

while non-dimensional Rayleigh number Ra, phase numbers Ra(410), Ra(660) and dissipative number Di are

Ra = [(a xpxgxd3xT1) / (rc x x)] = 5.55x108;

Ra(410)=[(5p(410)x g x d3) / (rc x x)]= 6.60x108; (4)

Ra(660) = [(5p(660)x gxd3) / (rc x x)] = 8.50x108; Di = [(a x gx d ) / Cp] = 0.165,

where a= 3x10-5 K-1 is the thermal expansion coefficient, p = 3.3 g x cm-3 is the density, g is gravity acceleration, Cp = 1.2*103 J x kg-1x K-1 is specific heat capacity at constant pressure, T1 = 1950oK is the temperature at the base of mantle transition zone (MTZ) at depth 660 km regarded the lower boundary of the model domain, Q = 6.25x10-4 mWx m-3 is the volumetric heat generation in the crust, m is the viscous stress tensor, d = 660 km is vertical dimension of the modeled domain, r]s = 1018 Paxs is the viscosity scaling factor, x = 1 mm2 x s-1 is thermal diffusivity, <p(410) = 0.07p and 8p (660) = 0.09p are

the density changes at the 410 km and 660 km phase boundaries respectively. In (1), (2) the scaling

_^

factors for time t, coordinates x and z, stresses m and the stream-function y are (d**x_1), d, nx' d and x respectively. Assuming rheology be linear for the diffusion creep deformation mechanism dominating in the mantle at depths over ~ 200 km [2], we accept the temperature- and lithostatic pressure p dependent viscosity as [15]

j = ( n / 2 * A) * (h / b*)m * { exp [ (E* + p* V*) / (R *T) ] },

(5)

where for "wet" olivine A = 5.3x1015 s-1, m = 2.5, the grain size h = 10-1 - 10 mm, b*= 5x10-8 cm is the Burgers vector [14], E = 240 kJ'mol-1 is activation energy, V=5x103 mm3'mol-1 is activation volume, H = 300 GPa is the shear modulus normalizing factor, R is the gas constant. At the constants chosen and the grain size h =1.6 mm, non-dimensional viscosity also denoted r is

r = 5 x 10-7x exp {[14.8 + 6.72 x (1 - z)] / T },

(6)

where T is non-dimensional temperature, non-dimensional z normalized by d is pointing upwards from the MTZ base and x is pointing against subduction along the MTZ base. The aspect ratio of the model domain is 1:2.25 thus the subduction angle being J3= 24° if subduction is assumed to take place along the model domain diagonal. Non-dimensional subduction velocity V = 10 mm*a_1 normalized by (d -1*x) equals V = 0.208*103, i.e. non-dimensional velocity components of subducting Adriatic micro-plate are Vx = - 0.190*103 and Vz = - 0.085*103.

To check as to how the estimate of velocity of subduction of the Adriatic micro-plate is sensitive to the accepted linear rheological law here we make extra computations for the non-Newtonian rheology, in which case the viscosity formulae (5)-(6) are rewritten as:

r = ( 1 / 2 x A x Cwr x tn-1) x (h / b*)m x { exp [ (E* + px V*) / (R xT) ] },

(7)

where according to [13] for "wet" olivine n = 3, r = 1.2, m = 0, t = (t2^)172 , E*= 480 kJ'mol"1, V* = 11x103 mm3xmol-1, A = 102 c_1x(MPa)"n, Cw > 10-3 for "wet" olivine is the weight water concentration (in %%). It should be noted the constants in (7) vary considerably in the papers referred to by [13] and heretofore we gave averaged values of constants. At Cw =10-3 on accounting for

Tik2 = (4 X ,2) x [ (^zz - wxx)2 / 2 + 2 x ¥xz2 ],

(8)

non-dimensional viscosity is

, = { 1.0 / [ (wzz - wxx)2 / 2 + 2 x ¥xz2 ]1/3} x exp{ [10.0 + 5.0 x (1 - z)] / T },

(9)

Following [13] we assume the phase functions T(0 as

r(l) = (1/2)x {1 - th[z - z(l)(T)] / w(l)}; z(l)(T) = zjl) - {[fl*(T - To(l))] / (p x g)},

(10)

where the signs are changed as z-axis is pointing upwards, z(l)(T) is the depth of the l-th phase transition (l = 410, 660), zo^ and To(l) are the averaged depth and temperature of the l-th phase transition, y(410) = 3 MPaxK-1 and y(660) = -3 MPaxK-1 are the slopes of the phase equilibrium curves, w(l is the characteristic thickness of the l-th phase transition, T0(410) =1800° K, To(660)=195 0° K are the mean phase transition temperatures. The heats of phase transitions are neglected in (2) as insignificant in the case of developed convection as in [13]. From (10) it follows

rx(l) = - (y(l) / 2xpxgxw(l)) x Tx x ch 2{[(z - zo(l)+ J(l) x (T-To(l)) ) / (pxg)] / w(l)},

wherefrom it is clear the phase transition with y(l) > 0 facilitates convection (at l = 410), while the phase transition with y(l < 0 hinders convection (at l = 660). In non-dimensional form Zo(410) = 0.38, Zo(660 = 0, w(l) = 0.05, y(410) = 2.5x109, y(660) = -2.5x109 , T/410) = 0.92, T/660) = 1 and in (1)

rx®=-(5p^ y(l)/2^p^Ra(l)^w(l))^Tx^ch-2{[z-zo(l)+y(l)(dp(l)/p-Ra(l))^(T-To(l))]/w(l)],

(12)

Equations (1)-(2) are solved for the isothermal horizontal and vertical boundaries regarded no-slip impenetrable ones except for the "windows" for in- and outgoing subducting plate, where the plate velocity is specified. Vertical boundary distant from subduction zone is assumed penetrable at right angle, the latter boundary condition appears not too imposing in the case of very flat subduction. Q in (2) is non-zero in the continental and oceanic crust 40 and 7 km thick. Initial vertical boundaries temperature is calculated for the half-space cooling model for 109 yr and 108 yr for Euro-Asian (continental) plate and Apulian lithospheric microplate (paleooceanic plate) respectively. It is convenient to express dimensionless T2tK in (2) through the stream-function w as in (8):

Tk2 = (4 x n2) x [ (¥zz - Wxx)2/ 2 + 2 x ¥xz2 ] Results and discussion

Assuming the heat flux maximum q is formed above the convective flow, upwelling towards the Pannonia basin and Vardar zone, and the convection cell dimension is equal to horizontal scale of the heat flux anomaly zone, the convection cell dimension can be estimated of ~300 km. To more accurately compute the consistent model of small-scale convection in the mantle wedge between the overriding Euro-Asian plate and subducting Apulian lithospheric microplate it is necessary from the computational point of view first to specify in (1)-(2) vanishing non-dimensional numbers Ra^0, Di = 0, i.e. to ignore convection and viscous dissipation. This approach is applied as convection with Ra and Di (4) passes through rather vigorous stages, and the time steps in integrating (1)-(2) become too small thus making it difficult to model the thermal structure of the plates. Solving (1)-(2) by the finite element method in space on the grid 104x104 and the 3-rd order Runge-Kutta method in time one obtains for Ra^-0, Di = 0 and V=10 mm a year non-dimensional quasi steady-state w and Tshown in fig. 1, where the streamlines are depicted with the step 5 and the isotherms with an interval 0.05.

2.2 2.0 1.8 1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0

Figure 1. Schematic cross section of the region of subduction of the Apulian lithospheric microplate under the Euro-Asian plate with no effects of viscous dissipation and convection taken into account. (1) - the quasi steady-state distribution of non-dimensional temperature and (2) - the quasi steady-state distribution of stream-function, the streamlines above the subducting slab correspond to the mantle

wedge flows ("A" and "B") induced by subduction.

Subducting plate was considered rigid, while the viscosity at the zone of plates friction (at temperatures below 1200o K) was reduced by 2 orders of magnitude as compared to (7). The latter viscosity reduction at the plates contact zone accounts for lubrication effected by deposits partially entrained by the subducting plate. Such a lubrication prevents the overriding Euro-Asian plate from gluing to the subducting one [9]. Fig. 1 shows the results of computation for the formulae (7) - (9) for non-Newtonian rheology case for the water content Cw = 3x10-1 weight %. The velocity V=10 mm per year is chosen as resulting in the best fit of the model convective zone size to horizontal extent of the observed heat flux anomaly in Pannonia basin and the Vardar zone. The Apulian lithospheric microplate subducting with a given velocity V is considered rigid and is shown in fig. 1(2) by the equidistant diagonal streamlines. The induced mantle wedge flow above the subducting plate is seen to occur in the form of two vortices A and B (located one above another), the latter 2 vortices being considerably compressed in the vertical direction and the upper one (with > 0) revolves clockwise while the lower one (with < 0) revolves counterclockwise. In fig. 1(2) the upper induced flow "A" is seen to be firmly pressed to the

subducting Apulian lithospheric microplate, the strain rate in the zone of contact of the opposite flows (i.e. flow "A" and subducting slab) being very high thus resulting in the viscosity (7) drop by several orders of magnitude. It should be noted that in the case of non-Newtonian rheology the greater the subduction angle the more extensive is this contact zone, in which the dissipative heat release mainly occurs. This may serve the reason why the Karig vortices (and the resulting back-arc spreading) are formed in the zones of comparatively steep subduction. The opposite flow "A" in fig. 1(2) apparently is induced by the flow "B", forced by subducting plate.

Assuming Ra = 5.55x108 and Di = 0.165, i.e. switching dissipation and convection on, taking into account the effects of phase transitions, from (1)-(2) the convection at Cw = 3x10-1 weight % is found to destroy the induced mantle flows in the mantle wedge during the time interval ~ 0.6*10"6 (in dimensional form ~ 0.1 Ma) and to assume the quasi steady-state form shown in fig. 2, in which the streamlines in the convection vortices are depicted with the interval 4x104

Distance, km

~ ' 1452. 1320. 1188. 1056. 924. 792. 660. 528. 396. 264. 132. 0.0

Figure 2. Quasi steady-state non-dimensional stream-function distribution in the zone of subduction of the Apulian lithospheric microplate under the Euro-Asian plate with the effects of dissipative heating and convection taken into account for non-Newtonian rheology for the water content Cw=3x10-1 weight %. Parallel equidistant streamlines represent subducting Apulian lithospheric microplate. Arrow "A" shows possible direction of the upward convective transport of dissipative heat and mantle

hydrocarbons to the Earth's surface.

These convective vortices are seen actually to correspond to a single convection cell aroused at subduction velocity V = 10 mm*«"1. The latter convection cell dimension is of the order of ~ 300 km, i.e. is very close to the observed horizontal extent of the heat flux anomaly observed in the Pannonia basin and Vardar zone. The velocity in convective vortices may exceed 10 m* a"1 for the water content of 3x10-1 weight %. The direction of a possible upward transport of mantle hydrocarbons and dissipative heat to the Earth's surface is shown by the arrow "A". It should be noted that in the case of Newtonian

rheology the convection in the mantle wedge cannot be aroused at the velocity of subduction of 10 mmx^r1 and subduction angle of 25°.

In fig. 3 the quasi steady-state stream-function ^ is shown for the mantle water content of 0.3x10" 1 weight %, in which case the dissipation-driven convection is aroused in essentially a single vortex with a characteristic velocity of ~10 mmxor1.

Distance, km

Zfl452 1320. 1188. 1056 924 792 660 528. 396 264 132 0.0

t 1 _ L _ _L _ 1 _ Л_I _ I _ I _ I _1 _I _ J « л

Figure 3. Quasi steady-state non-dimensional stream-function distribution in the zone of subduction of the Apulian lithospheric microplate under the Euro-Asian plate with the effects of dissipative heating and convection taken into account for non-Newtonian rheology for the water content Cw = 0.3x10-1 weight %. Parallel equidistant streamlines represent subducting Apulian lithospheric microplate. Arrow "A" shows possible direction of the upward convective transport of dissipative heat and mantle

hydrocarbons to the Earth's surface.

This convective flow may serve a means of the upward mantle hydrocarbons and dissipative heat transport, which direction is shown in fig. 3 by the arrow "A". Such a convective flow may cause the back-arc spreading in the Pannonia basin and the Vardar zone. At the mantle water content Cw < 0.3x10-1 weight percent the dissipation-driven convection cannot be aroused in the mantle wedge even in the case of the non-Newtonian rheology for subduction velocity of ~10 mmx^r1 and angle of 25°.

It is worth noting that in the case of Newtonian rheology the 2D transversal convective rolls in the mantle wedge, as in fig. 2, can be formed only at sufficiently small angles of subduction. Thus, at в = 30° the transversal convective rolls are not formed even at the velocity of subduction of 100 mmx a_1

[5, 9]. In the case of the non-Newtonian rheology the transversal rolls (2D Karig vortices) can be aroused at greater subduction angles and sufficiently small subduction velocities owing to viscous friction in the zone of contact of the opposite induced flow ("A" in fig. 1) and subducting slab. It should be noted that numerous thermo-mechanical mantle models in the zones of subduction (see, e.g. [8, 9] and the vast number of references there) showed convection in the form of transversal rolls never to occurre as the models with extremely small subduction angle and sufficiently great subduction velocity were not investigated.

Conclusions

The size of the cell of 2D mantle wedge dissipation-driven convection in the case of the realistic non-Newtonian rheology equals ~ 300 km at the subduction velocity 10 mm xor1 , in which case a single convection cell is aroused. This explains the formation and horizontal extent of the only 2D heat flux anomaly observed in the rear of the Dinarides. The water content sufficient for the 2D convection to be aroused is ~ 3x10-1 weight %, or, alternatively, it is ~3x10-2 weight %, but the 2D convection is aroused as a single Karig vortex. The velocity in convective vortices in the non-Newtonian rheology case is ~10 m per year at the water content Cw ~ 3*10-1 weight percent and ~10 mm per year at the water content Cw ~ 0.3*10-1 weight percent in the mantle wedge. The upwelling convective flow may be sufficient to provide upward transport of mantle wedge hydrocarbons to the Earth's surface.

References

1. Babbicci, D., C. Tamburelli, M. Viti, E. Mantovani, D. Albarello, F. D'Onza, N. Cenni, E. Mugnaioli. Relative motion of the Adriatic with respect to the confining plates: seismological and geodetic constraints // Geophys. J. Int. - 2004. - V. 159. - P. 765 - 775.

2. Billen, M., G. Hirth Newtonian versus non-Newtonian Upper Mantle Viscosity: Implications for Subduction Initiation // Geophys. Res. Lett. - 2005. - V. 32. - (L19304), Doi:10.1029/2005GL023458.

3. Carminati, E., M. Lustrino, C. Doglioni Geodynamic evolution of the central and western Mediterranean: Tectonics vs. igneous petrology constraints // Tectonophysics. - 2012. - V. 579. - P. 173 - 192.

4. Gavrilov, S.V. Investigation of the island arc formation mechanism and the back-arc lithosphere spreading // Geofizicheskie Issledovaniya. - 2014. - V. 15 (4). - P. 35-43.

5. Gavrilov, S.V., D. H. Abbott Thermo-mechanical model of heat- and mass-transfer in the vicinity of subduction zone // Physics of the Earth. - 1999. - V. 35. - P. 967-976.

6. Gavrilov, S.V., A.L. Kharitonov Velocity of the subduction of the Russian platform under the Siberian one in the Paleozoic by the distribution of mantle hydro-carbon upwelling zones in Western Siberia // Geofizicheskie Issledovaniya. - 2015. - V.16. - P. 36-40.

7. Gavrilov, S.V., A.L. Kharitonov Subduction velocity of the Russian plate under the Siberian one at Paleozoic: a constraint based on the mantle wedge convection model and the oil- and gas-bearing zones distribution in Western Siberia // Modern Science. - 2016. - V.16. - P. 155-160.

8. Gerya, T.V., J.A.D. Connolly, D.A. Yuen, W. Gorczyk, A.M. Capel Seismic implications of mantle wedge plumes // Phys. Earth Planet Interior. - 2006. - V.156. - P. 59-74.

9. Gerya, T.V. Future directions in subduction modeling // J. of Geodynamics. - 2011. - V.52. - P. 344-378.

10. Lenkey, L., P. Dovenyi, F. Horvath, S.A.P.L. Cloetingh Geothermics of the Pannonian basin and its bearing on the neotectonics // EGU Stephan Mueller Special Publication Series. - 2002. - N. 3. -P. 29 - 40.

11. Royden, L.H., F. Horvath, A. Nagymarosy, L. Stegena Evolution of the Pannonian basin system: 2. Subsidence and thermal history // Tectonics. - 1983. - N. 2. - P. 91 - 137.

12. Schubert, G., D. L. Turcotte, P. Olson Mantle Convection in the Earth and Planets // New York, Cambridge University Press, - 2001. - 940 p.

13. Trubutsyn, V.P., A.P. Trubitsyn Numerical model of formation of the set of lithospheric plates and their penetration through the 660 km boundary // Fizika Zemli. - 2014. - N.6. - P. 138- 147.

14. Zharkov, V.N. Geophysical Researches of the Planets and Satellites // Moscow: Joint Institute of Physics of the Earth RAS, - 2003. - 102 p.

15. Zharkov, V.N. Physics of the Earth's Interiors // Moscow: Nauka i obrazovanie, - 2012. - 384

p.