Full 3-D pseudo-transient finite difference modelling of stress distribution around continental plateaus

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. Due to mass conservation, thinning of the crust implies horizontal spreading of the plateau towards the lower altitude surroundings. This spreading is observed in GPS velocity records. For example, around the Tibetan Plateau, horizontal surface velocities are in the order of 2 cm/yr. Crustal flow also generates differential stresses in and around the plateau. Estimating these stresses and their spatial variation in 3-D is a computational challenge and necessitates new advanced computing methods.

 

Here, we present a new 3-D numerical algorithm to solve the Stokes equations under gravity. The algorithm is based on an Eulerian pseudo-transient finite difference method. To test the algorithm, we consider a simplified plateau geometry and density structure. Two different simulations are performed, one pseudo-2D simulation and one full 3-D simulation considering the corner region of a plateau. The pseudo-transient method allows an explicit solution of the Stokes equations. When the pseudo-transient time derivatives approach zero, a steady-state solution for the velocity field is obtained. In the first simulations, the rheology is linear viscous. The crust-mantle interface and the interface between crust and overlying sticky air are tracked with a Cahn-Hilliard diffuse interface model. To test our results, we compare the results of the pseudo-transient model with 3-D results obtained with the implicit numerical algorithm LAMEM (bitbucket.org/bkaus/lamem) and with 2-D results obtained with a Lagrangian finite element model. We have developed two versions of the algorithm, one in Matlab (mathworks.com) and one in Julia language (julialang.org). The aim of the Julia version is to eventually utilize a parallel GPU computing environment. Furthermore, we also present first results for cylindrical coordinates for the same plateau geometries. The aim of the model with cylindrical coordinates is to quantify the impact of the Earth’s curvature on the stress state.