Towards numerical modelling of natural subduction systems with an application to Eastern Caribbean subduction

My thesis aims to resolve two main outstanding problems that occur when trying to numerically simulate the 3D complexity of lithosphere subduction evolution through geological time: high computational cost and the efficient constructing of complicated 3D initial conditions mimicking a subduction setting in the geological past or at the present-day. The issue of high computational cost is mostly due to the required use of nonlinear rheologies which are needed to realistically include temperature-, pressure- and stress dependent material properties of rock deformation. The resulting nonlinear Stokes equations require an iterative solving strategy, which can require a lot of iterations with the simple to implement Picard iteration. The Picard iteration is robust, but slow to converge. The convergence rate can be greatly improved by using a Newton iteration, which may fail to converge if the iteration solution far from the real solution, but will converge very fast when it is near the real solution. We discovered that the linear system produced by the Newton iteration may in some cases not be solvable. We provide a general remedy for this issue by forcing the system to be solvable at a minimum amount of convergence loss. We also show that this solution works well for large 3D tectonic settings. The issue concerns setting up initial conditions for the geometric complexity of realistic 3D subduction settings involving creating 3D fields of initial temperature and of materials (e.g. crust, mantle, tectonic plates, faults and weakness zones, 3D slab geometry). An initial condition may concern the earliest stage of a subduction or other geodynamic system as well as an advanced stage, depending on our knowledge of the tectonic evolution. In this thesis, I developed a generic way for setting up initial numerical models leading to a stand-alone code library called the Geodynamic World Builder (GWB), which implements this. We show that the code library is already successfully coupled to three different geodynamic codes. The new generic approach together with the code library renders relatively easy initial-model construction as well as modification of initial settings. In addition, it creates a platform to make initial setups reproducible for use in different geodynamic codes. In the last main chapter, I use these new developments to investigate a complicated and realistic 3D subduction evolution that resembles the eastern Caribbean subduction setting since 10 Ma. The GWB allowed for flexible construction of the initial settings of tectonic plates, weakness zones, and 3D slab geometry at 10 Ma, while the Newton solver allowed for obtaining modelling results relatively fast. The modelling results show that slab dragging, i.e. lateral transport of the slab through the mantle by the trench-parallel motion of the North American plate, has a strong effect on the slab geometry and the slab stress field. Slab stress is displayed with the direction of maximum shear stress in the slab, which may align with potential fault slip directions of earthquakes. This shows that when slab dragging occurs major stress orientations are generally not trench-perpendicular, as is usually perceived for slab stress.