GPU-based finite-difference solution of 3-D stress distribution around continental plateaus in spherical coordinates

On Earth, different geodynamic features form in response to a tectonic event. Continental plateaus, such as the Tibetan Plateau, are formed in a collisional environment and they are characterized by an unusually large crustal thickness, which generates lateral variations of gravitational potential energy per unit area (GPE). These GPE variations cause the thickened crust to flow apart and thin by gravitational collapse. Although plateau and lowland are in isostatic equilibrium, the lateral GPE variations must be balanced by horizontal differential stresses, which prevent the plateau from flowing-apart instantaneously. However, the magnitude and distribution of differential stress around plateau corners for three-dimensional (3-D) spherical geometries relevant on Earth remain disputed. Due to the ellipticity of the Earth, the lithosphere is mechanically analogous to a shell, characterized by a double curvature. Shells exhibit fundamentally different mechanical characteristics compared to plates, having no curvature in their undeformed state. Understanding the magnitude and the spatial distribution of strain, strain-rate and stress inside a deforming lithospheric shell is thus crucial but technically challenging. Resolving the stress distribution in a 3-D geometrically and mechanically heterogeneous lithosphere requires high-resolution calcuations and high-performance computing.

Here, we present numerical simulations solving the Stokes equations under gravity. We employ the accelerated pseudo-transient finite-difference (PTFD) method, which enables efficient simulations of high-resolution 3-D mechanical processes relying on a fast iterative and implicit solution strategy of the governing equations. The main challenges are to guarantee convergence, minimize the iteration count and ensure optimal execution time per iteration. We implemented the PTFD algorithm using the Julia language. The Julia packages ParallelStencil.jl and ImplicitGlobalGrid.jl enable optimal parallel execution on mulitple CPUs and GPUs and ideal scalability up to thousands of GPUs.

The aim of this study is to quantify the impact of different lithosphere curvatures on the resulting stress field. To achieve this, we use a simplified plateau geometry and density structure implemented in a spherical coordinate system. The curvature is modified by varying the radius of the coordinates system, without altering the initial geometry. We particularly focus on stress magnitudes and distributions in the corner regions of the plateau.