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

Abstract

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 strainrate 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 im- * Corresponding author.Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. proved 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...

Tectonic drivers of the Wrangell block: Insights on fore‐arc sliver processes from 3‐D geodynamic models of Alaska

“…The 3‐D variable resolution finite element mesh, temperature structure, and plate boundary interface (Figure and Movie S1 in the supporting information) were made with SlabGenerator, a C/C++ code written by Jadamec and Billen [], Jadamec and Billen [], and Jadamec et al []. SlabGenerator uses the sea floor age grid [ Müller et al , ] and a half‐space cooling model [ Turcotte and Schubert , ] to determine the subducting plate thermal structure.…”

“…A lower‐viscosity shear zone represents a weak, poorly coupled boundary or fault zone, whereas a higher‐viscosity shear zone represents a strong, highly coupled plate boundary or fault zone. CitcomCU was modified to read in and implement narrow 3‐D shear zones representative of the plate boundary and of intracontinental fault zones within the upper plate [ Jadamec et al , ]. Each shear zone is blended into the background viscosity, η eff , using a sigma function, such that the shear zone viscosity, η wk , is defined by where η 0 is the reference viscosity (Table ) and A wk is a scalar weak zone field delineating the shear zone [ Jadamec and Billen , ; Billen and Jadamec , ; Jadamec and Billen , ; Jadamec et al , ].…”

A Visual Survey of Global Slab Geometries With ShowEarthModel and Implications for a Three‐Dimensional Subduction Paradigm

“…As data assimilation and working with numerical modeling output from large‐scale parallel simulations become an integral part of modern scientific research, new tools are required to analyze and conceptualize the information (Anderson et al, ; Billen et al, ; Damon et al, ; Erlebacher et al, ; Jadamec et al, ; Kellogg et al, ; Kreylos et al, ; Suarez et al, ; Wessel et al, ; Zhou & Liu, ). With the new digital framework for research, scientific contributions are taking the form of multiple formats from contributed software (Hwang et al, ; Zhong et al, ) and workflows (Jadamec, ), to large‐scale data sets and integrated parsers (Anderson et al, ; Haxby et al, ; Hutko et al, ; Trabant et al, ), to fully visual formats to communicate the digital and 3‐D information, as in this paper.…”