PROBLEMS OF NUMERICAL MODELING OF LARGE-SCALE MANTLE CONVECTION IN THE SUBDUCTION ZONE

1. INTRODUCTION

In geodynamics, subduction is "the process by which an oceanic lithospheric plate/slab sinks into the mantle and moves under a continent or an island arc" [Stern, 2002, р. 5] where the main driving mechanisms of the Earth dynamics are thermal energy within the Earth [Khain, 2010] and density contrast between cold subducting plate and surrounding mantle. The reliability of this statement is confirmed by the fact that the mantle substance behaves as a solid at rapidly changing loads but can flow as a viscous fluid at long-term loads [Monin, 1977]. According to the theory of lithospheric plate tectonics, slab is a seismotomographically identified fragment of an 80–100 km thick oceanic lithospheric plate sinking (subducting) into the mantle [Koulakov et al., 2011; Sizova et al., 2014].

The importance of studying interaction between mantle convection and subduction lies in the understanding of the pattern of the distribution of substances in the geosphere and the structure and properties of the peripheral envelopes of the Earth. Subduction zones are one of the most tectonically active parts of the Earth related to large-scale mass and energy transport [Stern, 2002; Frost, 2006]. Of particular interest is the study of the processes in rheologically contrasting media: convective in relatively low-viscous asthenosphere and conductive in rigid lithosphere.

Though the numerical models of upper mantle convection have a wide coverage [Cristensen, 1984; Schubert, Anderson, 1985; Zhong, Gurnis, 1992; Dobretsov, Kirdyashkin, 1997; Lithgow-Bertelloni, Richards, 1997; Rychkova, Tychkov, 1997; Van Keken et al., 2011; Gerya, Yuen, 2003], there are still many open problems concerning changes in mantle convection structure induced by slab subduction; discontinuity and fragmentation of subduction; accumulation of slab material at the upper-lower mantle boundary; adaptation of modern computing methods for representation of upper mantle convection related to subduction zone. It is accepted that the lithosphere is entirely underlain by the asthenosphere whose viscosity is one or two orders of magnitude less than that of the lithosphere [Dobretsov et al., 2001]. But it is here that there are the conditions for melting of lithospheric material and generation of deep-focus earthquakes [Trubitsyn et al., 1993].

The aim of this work is the analysis of modern methods for solving some of the above-mentioned problems and the construction of 2D subduction dynamics models that take into account phase transitions in the subducting slab which cause a significant slab density increase. Consideration is being given to the schemes for solving nonlinear model equations and to the problems of software implementation of numerical modeling apparatus.

2. MANTLE CONVECTION AND SUBDUCTION MODELS

The study of subduction assumes that the character of interaction between mantle convection and slab is as follows. The lithospheric model description provided in [Zanemonets et al., 1974; Myasnikov, Markaryan, 1977; Bercovici et al., 1993] incudes highly viscous conductive quasi-fluid layer with zero velocity of fluid substance motion. The other part of the description consists of incompressible fluid whose rheological properties are characterized by the form for recording viscosity values and accepted initial and boundary conditions. Mantle flows usually assume constant slip at the computational area boundaries. However, conditions can arise for cyclic activation of magmatism at the asthenosphere-lithosphere boundary.

Mantle convection assumes to have chemical-dense nature [Keondzhyan, 1980; Monin, Sorokhtin, 1981]. Then, due to high viscosity, its modeling usually involves the Stokes equation in the Boussinesq approximation [Slezkin, 1955; Landau, Lifshits, 2001]. Accepting this approximation provides an opportunity to introduce the temperature factor into the equation for vorticity and, therefore, for velocity.

The tasks are solved for different forms of representation of viscosity. In terms of vorticity–stream function, heat convection was studied in viscosity continuum [McKenzie et al., 1974] and in the continuum where viscosity increases exponentially with depth [Gurnis, Davies, 1986].

A dynamic model of slab where it is represented by a negatively buoyant thin elastic plate is considered in [Turcotte, Schubert, 1985]. It is assumed that the plate bending is related to external vertical loads and ratios of its bending moments. Numerical model of subduction consists of the dimensionless equation for plate bending under the effect of applied moments and vertical loads. It is also assumed that a subducting plate can carry water deep into the upper mantle [Shubert et al., 2001]. This approach was adapted for the Amur microplate subduction through the upper mantle in [Gavrilov, Kharitonov, 2022], where phase transitions correspond to the key points in [Trubitsyn V.P., Trubitsyn A.P., 2014]. The eligibility of such approach is only limited by the first upper mantle layers where Kirchhoff’s hypothesis of normals remaining normal to the deformed middle surface of the plate and an unchanging plate thickness is applicable [Timoshenko, Woinowsky-Krieger, 1966]. However, the results of seismic sounding show that there is a decrease in the plate thickness [Fukao et al., 2009]. With plate subduction, there is a corresponding increase in pressure and temperature which leads to partial loss of light plate substance fractions, losses of water and gases from hydrous minerals and mineral lattice rearrangements, subsequent phase transitions, and slab fragmentation during deep subduction.

Alternative 2D models of upper mantle convection in the subduction zone are based on the statements of continuum mechanics and hydrodynamics where the equations of motion are expressed in the Stoke approximation [Gurnis, Davies, 1986; Doglioni et al., 2006; Van Keken et al., 2011; Duretz et al., 2012]. In these models, the process of subduction is mainly caused by excessive slab density. The mantle density follows the Boussinesq approximation. Some models differ in the form features of recording mantle and slab viscosity values, dependence of slab dynamics on temperature and tangential stress values [Van Keken et al., 2011], and in ways of impact on the plate and compositional variations [Kincaid, Sacks, 1997].

In [Gerya, Yuen, 2003], there is presented a 2D model of thermomechanical mantle convection in the Boussinesq approximation where the upper mantle viscosity is determined by temperature and pressure, and the model equations follow the Stokes approximation. Radioactivity and adiabatic and dissipative heating in the mantle act as independent spatiotemporal heat sources which determine the kinetics of slab subduction into the mantle. The computations are also based on using the marker and cell method. The solution is interpolated back to the Eulerian grid configuration at each time step. Such an approach requires considerable computational resources. Suffice to say that modeling of the Precambrian continental collision [Zakharov et al., 2015] was performed using the Lomonosov Moscow State University supercomputer complex.

In [Dobretsov, Kirdyashkin, 1997; Dobretsov et al., 2001, 2009; Koulakov et al., 2011], a study was made of thermophysical characteristics of only crustal layer which is a boundary layer between the mantle and the slab. However, these models do not consider the influence of mantle flows on slab dynamics.

In [Trubitsyn, 2008], consideration is being given to the case when the equation of temperature of the endothermic phase transition segment at the upper-lower mantle boundary is complemented by an analytic representation of a phase boundary shift/phase distribution function It is shown that an endothermic phase transition at a depth of 660 km does not lead to convection stratification. There has been observed slowdown and deformation of abnormally cold descending fluids. A slab model is a field of low temperature environment, and slab dynamics is determined by the same equations of mantle convection. The lack of a slab model as a representation of a real object of the upper mantle – slab system may be due to the drawback of this approach.

In [Bobrov, Baranov, 2014], the viscosity ratio considers deformation plasticity which results in abrupt mode of subduction and abrupt changes in stress fields depending on the slab separation stage. In the areas which do not contain subducting slabs the stresses are very low. The overall picture for temperature dynamics shows the slab segments as dense isothermal patterns. The lack of a slab model as a representation of a real object of the upper mantle – slab system may be due to the drawback of this approach.

In [Lobkovskii, Ramazanov, 2021], Galerkin method application and Fourier series expansion provide a basis for an analytical representation of a stationary solution to model equations of convection in the mantle whose substance is determined by a cold fluid. There are presented the results of computational experiments for structural variations in convective cells depending on the Rayleigh number. The lack of a slab model as a representation of a real object of the upper mantle – slab system therein may also be due to the drawback of this approach.

The above-mentioned publications do not consider the problems of transformation of hydrous minerals of the slab that occurs due to squeezing water out of the upper mantle rocks. Besides, in many cases the slab is not regarded as an independent geodynamic object.

3. SLAB FORMATION AND SUBDUCTION DYNAMICS

A model of mantle dynamics in the subduction zone consists of interrelated systems of equations, the first of which determines mantle dynamics, and the second – a colder and a denser oceanic plate subduction against the background of the converted mantle. The mantle is modeled by a two-dimensional incompressible Newton fluid in the Cartesian coordinate system. Due to high viscosity of the mantle, the quasi-stationary two-dimensional convective fluids are determined by the well-known Stoke equations, whose dimensionless form in terms of vorticity – stream function is derived as [Torrance, Turcotte, 1971; Turcotte, Schubert, 2002]:

(1)

(2)

(3)

(4)

where dimensionless variables stand for: x, y – axes of the Cartesian coordinate system (x=y=0 left corner), y-axis is directed downward;

 – 2D Laplace operator;

μ(T,y), ξ(x,y,t), ψ(x,y,t), ρ(x,y,t) – dynamic viscosity of the mantle, vorticity, stream function of mantle flow and mantle density; U(x,y,t)=ψ_y, V(x,y,t)=–ψ_x, T(x,y,t) – lateral and vertical velocities of mantle flow, temperature distribution in the mantle;

 – the Rayleigh number for upper mantle;

 – gravitational acceleration, typical densities and velocities of mantle flows, coefficient of thermal expansion and upper mantle depth; dT=T_max–T_min – drop in temperature typical of the upper mantle; χ_ρ – diffusion coefficient as a function of density for rocks; lower indices of the variables show at the corresponding partial derivatives of the same name.

The velocity-vorticity relation (1) allows avoiding the obvious presence of pressure.

Nondimensionalization of (1)–(3) was performed in a standard way (the primes denote dimensionless variables):

where the variables with upper bars are typical values; λ and χ are thermal conductivity and thermal diffusivity of the mantle, and С_Р is isobaric specific heat capacity; the subsequent dimensionless equations omit the primes.

The computational domain is a two-dimensional rectangular region of the Cartesian coordinate system. The boundary conditions imply the absence of heatflow and substance on the lateral faces of the computational domain and the presence of heatflow and substance on the upper g_T,0, g_ρ,0 and lower g_T,H, g_ρ,H faces of the domain:

(5)

where x_max is a dimensionless lateral extent.

Stream function assumes no-slip condition on the lateral extents: ψ(0,y,t)=ψ(x_max,y,t)=0, ψ_n(0,y,t)=ψ_n(x_max,y,t)=0, where n is a normal [Lobkovsky et al., 2013, 2021]. The portion of the upper mantle where an oceanic plate moves towards the subduction zone, as well as the entire bottom of the mantle, assume slip condition:

when y=0, interval X is a horizontal portion of the upper mantle where the plate moves towards its subduction zone, and when y=1, interval X is lateral component x; at solid boundaries (including the boundaries with the slab), stream function is equal to zero; the vorticity value is determined by Thom’s formula [Roache, 1980]

(6)

Analysis of seismotomographic profiles [Stern, 2002; Fukao et al., 2009] shows that the extension of the slab is much greater than the thickness which allows providing its model representation via a flexible thin bar and implies the independence of horizontal subduction velocity on vertical y coordinate and that of vertical subduction velocity on x coordinate. It is also assumed that the slab is in hydrostatic state, and the slab density can be calculated from the density of the mixture of heavy ρ_H=const and light ρ_L=const components (solid solution). The slab subduction occurs at ρ^*–ρ˃0. Then the subduction model can be determined by the following dimensionless equations

(7)

(8)

(9)

(10)

(11)

The B(C,ρ^*) formulation set forth herein assumes the following. Seawater seeps into oceanic lithosphere through fractures and pores, and reacts with minerals in the crust and mantle to form hydrous minerals (such as serpentine) that store water in their crystal structures [Turcotte, Schubert, 1985]. That is why during the slab subduction, there are various transformations of such minerals (in particular, the formation of more compacted units of minerals due water squeezing) which can be interpreted as a partial phase transition of the slab substance. In this situation, large amounts of gas- and water-saturated formations (fluids) accumulate at the mantle – slab boundary [Dobretsov et al., 2017], which subsequently can cause seismicity and partial melting of the subducting slab substance [Gerya et al., 2004]. Then, during the slab subduction, there is an increase in pressure and temperature of a medium which causes a removal/ascent of a light component, resulting in an increase in heavy component concentration and a subsequent weighting of the slab substance. The mechanism of this process is as follow. At early stages of subduction of the slab, when its substance is oversaturated with seawater (and, consequently, light components) seeped therein, an increase in pressure and temperature of a medium promotes an intensive ascent and a subsequent formation of excess volumes of high-temperature and gas-saturated components. Before the slab density becomes similar to the density of the medium, the slab subduction slows down, thus forming stagnant areas of accumulation of the slab substance. The fact is that the effect of mantle flow and an irregular distribution of the slab density lead to a difference in stresses and subsequent deformations thereon and, finally, to slab fragmentation. There is so called slab stagnation whose boundaries are clearly defined along the hypocenters of deep-focus earthquakes. Their maximum recorded depth, particularly under South America and Central America, does not exceed 700 km [Stern, 2002; Fukao et al., 2009]. Then, with an increase in heavy component concentration in the slab, there is a break of the stagnation zone, with a subsequent subduction of the fragments down to the Earth’s core [Agrusta et al., 2017].

A simple record of such mechanism herein is determined by relation:

(12)

The upper-mantle boundary condition for subduction model (6)–(10) is determined by an initial segment of subduction where typical subduction velocity is assigned. At the upper mantle bottom there is an accumulation of the slab substance.

Seismic tomography of the upper mantle shows the mantle – slab boundary irregularities which are mainly due to the influence of mantle flow on the slab. The thickness of the slab therewith is much less than the extension, which allows providing its model representation via a flexible and thin bar, affected by the mantle flow. In this situation, the need for getting partial derivatives to solve (7)–(10) in a curvilinear coordinate system with dynamic and irregular boundaries allows the use of adequate meshfree methods. Here it is appropriate to use the Lagrangian meshfree particle method (SPH) [Monaghan, 2005]. This method is based on representation of the object of study (slab herein) as a finite set of elements/particles [Harlow, 1967], carrying density, spatial location and velocity properties, and as a set of additional attributes.

The essence of the method lies in the fact that a certain physical characteristic of a single particle in point is determined by the weighted sum of this value for other particles

(13)

which lie in neighborhood are mass and density of the particle with number а, and h (so called smoothing length); is an analytic kernel through which the particle properties are smoothed and their mutual influences are taken into account. The smoothing function determines the amount of impact of each particle on the calculated value. As the particle, for which the attribute value is calculated, is closer to a nearby particle, the effect of the latter on the result becomes more pronounced, and the smoothing function increases more rapidly is this position [Aliev, 2008; Medin, Parshikov, 2010]. In this situation, the calculation of partial derivatives

comes down to their calculation for analytic function:

where b, a are ordinal numbers of particles, b, a=1÷N(t); N(t) – is the current amount of particles in the method; x_ba=x_b–x_a, y_ba=y_b–y_a and  – is kernel function which provides the basis for analytical calculations of the derivatives set in square brackets. The Laplacian of  is provided after [Brookshaw, 1985]. For the kernel function use has been made of the 3^rd degree spline curve

at 0<R<h,

at R>h;

[Monaghan, 2005].

With the SPH method the heat transfer formula for the b^th particle is given by

(14)

where c_P,b, T_b, T_a, k_b, k_a, x_a, y_a are heat capacity, temperature, thermal conductivity coefficient and coordinates for corresponding particles. For t+1 moment, they are determined by:

(15)

where are lateral and vertical particle velocities determined by the SPH method from (7) and (8) with the given initial and boundary conditions.

4. COMPUTATIONAL EXPERIMENTS

A regular computational grid:

is generated as follows: x_i=x_i–1+dx, y_i=y_i–1+dy where dx and dy are step sizes along the axes; N_x and N_y are the amounts of steps along the coordinate axes; a dimensionless temporal step dt is assumed const. The pressure is assumed to be hydrostatic, and the relationship between dimensionless viscosity for the upper mantle and slab is derived as:

where a_T=–0.2 and b_Y=1.33 [Trubitsyn et al., 1993; Bobrov, Baranov, 2014].

Initial distributions for temperature and density follow equations:

where T^(L)(x,y) and ρ^(L)(x,y) were derived based on the sources referred to herein [Dziewonski, Anderson, 1981; Rychkova, Tychkov, 1997; Sorokhtin, Ushakov, 2002; Honda, Yuen, 1994].

Analysis of Fig. 1 shows the smooth flow where the stream function values in the cell center are seven times higher than those at the domain boundaries. There is a slight variation in heat and mass transfer components due to the balance of their flows between the geospheres. Under the statement of the problem, no-slip condition is fulfilled and the stream function is equal to zero on the domain boundary. Since the mantle – slab boundary is impermeable, the stream function also turns to zero thereon which causes smooth flowing around the slab profile.

The distribution of temperature T, density ρ and stream function ψ after 1500 temporal iterations (more than 48 million years) is shown in Fig. 2. For convenience and clarity, the dynamics of the spatiotemporal distributions of temperature and density is presented as deviations from corresponding values for t=500.

The comparison between temperature and density for different temporal iterations shows some cooling and softening of the bottom layers in the upper mantle which is assumably due to non-zero flows of heat and substance at the upper mantle boundaries (condition (5)). The analysis of Fig. 2, b, shows somewhat non-linear density distribution: between the layers with negative density deviation for dimensionless t=500, there is a layer with positive density deviation. This is caused by convective mixing in the mantle. Also of note are temporal evolution of the stream function and, therefore, the impact on the slab which gives rise to the occurrence of different deformations in the slab body and subsequently leads to its partial fragmentation.

The distribution of temperature, density and streamlines after 1900 temporal iterations (more than 91million years) is shown in Fig. 3.

The distribution of temperature shows some cooling of the bottom layers in the upper mantle which is attributed to the heat efflux. Also of note is some softening of the bottom layers due to convective mixing. As a result of plate sinking into the mantle (subduction), the initial convective cell is split into two components. The slab therein acts as a real vertical "partition". The distribution pattern of isolines in Fig. 2, c, and Fig. 3, c, shows that in the immediate vicinity of the slab there are conditions for the compaction of the stream function values and the formation of small-scale (Karig) vortices [Gavrilov, Kharitonov, 2017]).

As the slab approaches the bottom of the lower mantle, it appears to be flattening out which is due to an increase in its temperature and decrease in its viscosity therein. Together with denser underlying layers, this leads to deformation of the frontal part of the slab under the effect of convective flows.

The manifestation of features of the formation of small-scale vortices, as well as that of flattening out and fragmentation of the slab, becomes illustrative in the scaled neighborhood of the slab in Fig. 4.

The effect of mantle flows on the slab and an irregular density distribution thereon cause its partial fragmentation. The slab fragments reach the top of the lower mantle and accumulate there (diffusion of particles at this boundary is assumable).

5. CONCLUSION

The paper deals with the key principles for modeling subduction in the upper mantle. Consideration is being given to a complex model consisting of a model of mantle convection (vorticity equations, stream functions, heat and mass transfer equations) and a spatiotemporal model of slab dynamics (velocity equations, heat and mass transfer equations, adapted SPH method for solving the problem). The paper proposes a model for quantitative accounting of phase transitions in the slab substance.

The effect of mantle flows on the subduction dynamics is taken into account by supplementing the slab model with the sets of identical particles from the slab-adjacent mantle cells.

The effect of mantle flows on the slab and an irregular density distribution thereon causes its partial fragmentation.

6. DISCLOSURE

The author has no conflicts of interest to declare. The author read and approved the final manuscript.