Three-dimensional simulations of geometrically complex subduction with large viscosity variations

The incorporation of geologic realism into numerical models of subduction is becoming increasingly necessary as observational and experimental constraints indicate plate boundaries are inherently three-dimensional (3D) in nature and contain large viscosity variations. However, large viscosity variations occurring over short distances pose a challenge for computational codes, and models with complex 3D geometries require substantially greater numbers of elements, increasing the computational demands. We modified a community mantle convection code, CitcomCU, to model realistic subduction zones that use an arbitrarily shaped 3D plate boundary interface and incorporate the effects of a strain-rate dependent viscosity based on an experimentally derived flow law for olivine aggregates. Tests of this implementation on 3D models with a simple subduction zone geometry indicate that limiting the overall viscosity range in the model, as well as limiting the viscosity jump across an element, improves model runtime and convergence behavior, consistent with what has been shown previously. In addition, the choice of method and averaging scheme used to transfer the viscosity structure to the different levels in the multigrid solver can significantly improve model performance. These optimizations can improve model runtime by over 20%. 3D models of a subduction zone with a complex plate boundary geometry were then constructed, containing over 100 million finite element nodes with a local resolution of up to 2.35 km, and run on XSEDE. These complex 3D models, representative of the Alaska subduction zone-transform plate boundary, contain viscosity variations of up to seven orders of magnitude. The optimizations in solver parameters determined from the simple 3D models of subduction applied to the much larger and more complex models of an actual subduction zone improved model convergence behavior and reduced runtimes by on the order of 25%. One scientific result from 3D models of Alaska is that a laterally variable mantle viscosity emerges in the mantle as a consequence of variations in the flow field, with localized velocities of greater than 80 cm/yr occurring close to the subduction zone where the negative buoyancy of the slab drives the flow. These results are a significant departure from the paradigm of two-dimensional (2D) models of subduction where the slab velocity is often fixed to surface plate motion. While the solver parameter optimization can improve model performance, the results also demonstrate the need for new solvers to keep pace with the demands for increasingly complex numerical simulations in mantle convection.