INTRODUCTION

One of the most important geophysics problems is to obtain a self-consistent geodynamic model of the modern Earth. It is necessary to explain the origin of global tectonic structures such as mid-ocean ridges, continental rifts, volcanic provinces, collision belts, subduction zones, lithospheric plates, and other structures. The theoretical basis of global geodynamics is to study the mantle convection of the real Earth. This study implies, at a minimum, the solution of the Stokes equation, which describes a three-dimensional flow of a viscous fluid dependent on the viscosity–temperature relationship. We do not consider the geodynamic evolution of the Earth, which implies the solution of, along with the Stokes motion equation, the heat transfer equation, as well as the equations for the transfer of chemical components. Our study is restricted to the modern “instantaneous” pattern of flows in the mantle based on a three-dimensional pattern of the Earth’s seismic tomography recalculated into the temperature and density distribution in the mantle.

The mathematical modeling of three-dimensional mantle convection in the spherical Earth for an arbitrary temperature field and variable viscosity is possible only by numerical methods [1]. In recent years, many numerical experiments have been carried out to study convection in the mantle in two- and three-dimensional versions in terms of various parameters of the Earth and in relation to different geodynamic settings (for example, [2–8]). The instantaneous 3D geodynamic models of the Earth were also calculated based on the seismic tomography, for example [9, 10]. Based on the SMEAN 2 global seismic tomography model [11], we construct a geodynamic spherical model of the modern Earth and present its consequences for the Antarctic region concerning potentially dangerous processes of accelerated glacier sliding into the ocean and glacier destruction processes observed in recent decades.

NUMERICAL MATHEMATICAL MODEL

The instantaneous spherical structure of global viscous mantle flows can be calculated using the Stokes equation based on the seismic tomography data. The SMEAN 2 seismic tomography model [11] we used is one of the best and contains variations in the shear seismic velocities throughout the entire mantle. In order to use the data of this model in numerical calculations, it is necessary to convert shear wave velocity variations in the mantle into temperature anomalies. Variations in seismic velocities in the mantle Δ\({{{v}}_{{\text{s}}}}\) are recalculated into the material density according to the relation Δρ = 0.3 × Δ\({{{v}}_{{\text{s}}}}\), based on the geophysical considerations [12]. In turn, the field of temperature variations in the mantle is determined from density variations using the thermal expansion formula: ΔT = –(1/α) × (Δρ/ρ), taking into account the relationship between the thermal expansion coefficient α and the depth. In this case, the thermal expansion coefficient α varies with depth according to the dependence α = (3‒4.44 × (1 – r)) × 10–5, where r is the dimensionless radius of the Earth. Thus, the thermal expansion coefficient α changes from 3 × 10–5 on the Earth’s surface to 1 × 10–5 at the mantle bottom at the boundary with the core. When the adiabatic and potential temperature are added to the temperature variations thus calculated, the total temperature in the Earth’s mantle is obtained.

We model the Earth’s mantle in the Boussinesq approximation in the three-dimensional spherical geometry. The mantle is heated from the core and from the inside due to the decay of radioactive elements (internal heating). In the calculations, we use the model with a simplified viscosity–temperature relationship according to the Arrhenius law [13]:

$${{\eta }_{T}} = \exp (2E{\text{/}}(T + {{T}_{{{\text{bot}}}}}) - 2E{\text{/}}({{T}_{{{\text{ref}}}}} + {{T}_{{{\text{bot}}}}})),$$

where E is a dimensionless parameter that determines the viscosity drop in the model, T is a dimensionless superadiabatic temperature, and Tref = 0.5, Tbot = 1 (temperature at the mantle bottom). In this work, E = ln104.5 = 10.36, which roughly corresponds to the activation energy of wet olivine. In our model, the jump in viscosity at the boundary between the upper and lower mantle is taken as 30. The presented viscosity law yields variations in viscosity, both in depth and along the lateral, by several orders of magnitude. Near the surface, a weak oceanic lithosphere is formed due to the low temperature, while the continents present highly viscous regions down to a depth of 200 km (due to the lower temperature). The Rayleigh number is Ra = 1 × 108. In the numerical simulation, we used a modified version of the CitcomS software [14]. The code, which use data from the SMEAN 2 seismic tomography model, was added to the software. In this model, the data were decomposed into spherical harmonics so as not to depend on the computational grid. The calculations were made with a grid of 170 × 170 × 59 nodes in angles and depth, respectively, with a uniform depth step of 50 km. In the CitcomS software we used, the grid by the angles is uneven and consisted of 12 spherical segments. For such a grid, there are no singularities in the polar regions. The grid is detailed in [14].

At the input, 58 files are recorded: they contain seismic velocity variations in the Earth distributed by the spherical harmonics and those converted into temperature variations at each grid point. Thus, the initial temperature field is formed. Further, the adiabatic and potential temperature are added to it. Then, the momentum transfer equation (Stokes equation) is solved for the flow velocities in the natural variables such as velocity and pressure by the finite element method using the Uzawa algorithm [15]. This approach makes it possible to obtain a solution even when the material viscosity changes by many orders of magnitude.

RESULTS AND DISCUSSION

For the purpose of interpreting the results and representing them in graphical form, we chose a characteristic NS cross section between the poles, passing through East Africa and the central part of the Pacific Ocean (40°, 220° E) (Fig. 1). Figure 1 shows the temperature variations in the mantle in the spherical NS section of the Earth through the poles at 40° and 220° E. It can be seen that global hot ascending flows rise from the Earth’s core under the Pacific Ocean and Africa. A weaker and partially upper mantle flow rises to the surface in the Arctic region. The descending flow between Africa and Eurasia (Fig. 1) corresponds to the subduction zone in the Mediterranean Sea. South of Alaska, the descending flow corresponds to the Aleutian subduction zone. There is a descending mantle flow under the East Antarctic region.

Fig. 1.
figure 1

Distribution of temperature and velocity anomalies in the Earth’s mantle in the cross section along 40° and 220° E.

Full size image

Figure 2 shows temperature variations in the Antarctic region at a depth of 100 km. Black contours show the dynamic topography, while black dots indicate volcanoes. The East Antarctic region is characterized by high negative temperature anomalies of up to 400° and a dynamic topography of up to –900 m. These data confirm a descending mantle flow under the East Antarctic region.

Fig. 2.
figure 2

Temperature variations in the Antarctic region at a depth of 100 km. Black outlines show a dynamic topography (m), while black dots indicate volcanoes.

Full size image

The Antarctic Peninsula, the West Antarctic Rift System, the Ross Ice Shelf, and the Ross Sea are characterized by positive temperature anomalies of up to 350° and a positive dynamic topography of up to 600 m. These features correspond with the data on the increased surface heat flow and recent volcanism in the West Antarctic Rift System and part of the Antarctic Peninsula. Meanwhile, in another part of the West Antarctic region, the Filchner–Ronne Glacier, the temperature anomalies and dynamic topography are negative. This result is consistent with the observed evidence of no rifting at present, a measured small heat flow under this area, and the absence of volcanoes.

CONCLUSIONS

The results obtained by numerical modeling based on the SMEAN 2 seismic model of variations in shear seismic velocities yield the temperature and viscosity distribution and the instantaneous structure of global mantle flows calculated using the Stokes equation. The calculation results are consistent with the theory of plate tectonics. Descending mantle flows and negative temperature anomalies occur under the continents, except for the regions of East Africa, Southeast and East Asia, and the West Antarctic. Under the East Africa region, there is a positive temperature anomaly and an ascending upper mantle flow responsible for the rift system on the surface of the African continent. A similar high-temperature anomaly and rising currents were also found in the West Antarctic region. Global mantle ascending flows are under the Pacific Ocean and South Africa. The Earth’s mantle consists of two regions with high lateral differences in temperature and viscosity: (1) the difference between the hot asthenosphere and the cold continental lithosphere of cratons at a depth of about 200 km and (2) the difference between hot ascending plumes and cold zones of descending slabs in the D'' layer at the boundary with the core (Fig. 1).

The Antarctic region is clearly divided into two parts by the temperature anomalies and the dynamic topography. The West Antarctic region, except for the Filchner–Ronne Glacier and the northern part of the Transantarctic Mountains bordering on the Ross Sea, is characterized by positive temperature anomalies in the sublithospheric mantle and a positive dynamic topography. Meanwhile, the temperature is lower in the sublithospheric mantle and the dynamic topography is negative under the East Antarctic region and the Filchner–Ronne Glacier (Fig. 2).

The numerical results obtained above are in good agreement with the observed heat flow on the surface of the Antarctic area [16] and a high range of volcanoes in the West Antarctic Rift System and part of the Antarctic Peninsula [17]. Hence, the calculated model of modern mantle flows provides a completely logical explanation for the West Antarctic Rift System and many volcanoes in the West Antarctic region. The increased heat flow and volcanism under the West Antarctic Rift System contribute to melting of the West Antarctic Ice Sheet base thus facilitating the ice sliding from the West Antarctic interior into the sea along the bedrock. The rapid glacier sliding and destruction against the background of an increased heat flow is due to the trigger effect of the activated subglacial volcanoes caused by deformation waves arriving at the West Antarctic region; these waves are generated by strong earthquakes in the subduction zones surrounding the Antarctic area [18–20]. This process can lead to the large-scale rapid sliding of sheet glaciers into the sea (for example, the Doomsday Glacier) and a global rise in the World Ocean level by a few tens of centimeters and even a few meters.

The presented global model of mantle flows requires further refinement. In particular, our model does not take into account variations in the chemical composition in the mantle, continental crust, and D'' layer. More detailed calculations of the flow structure in the Earth’s upper mantle, the introduction of plate rheology on the Earth’s surface, etc., are needed.