Margarete A. Jadamec

Reconciling surface plate motions with rapid three-dimensional mantle flow around a slab edge

The direction of tectonic plate motion at the Earth’s surface and the flow field of the mantle inferred from seismic anisotropy are well correlated globally, suggesting large-scale coupling between the mantle and the surface plates. The fit is typically poor at subduction zones, however, where regional observations of seismic anisotropy suggest that the direction of mantle flow is not parallel to and may be several times faster than plate motions. Here we present three-dimensional numerical models of buoyancy-driven deformation with realistic slab geometry for the Alaska subduction–transform system and use them to determine the origin of this regional decoupling of flow. We find that near a subduction zone edge, mantle flow velocities can have magnitudes of more than ten times the surface plate motions, whereas surface plate velocities are consistent with plate motions and the complex mantle flow field is consistent with observations from seismic anisotropy. The seismic anisotropy observations constrain the shape of the eastern slab edge and require non-Newtonian mantle rheology. The incorporation of the non-Newtonian viscosity results in mantle viscosities of 1017 to 1018 Pa s in regions of high strain rate (10-12 s-1), and this low viscosity enables the mantle flow field to decouple partially from the motion of the surface plates. These results imply local rapid transport of geochemical signatures through subduction zones and that the internal deformation of slabs decreases the slab-pull force available to drive subducting plates.

Enabling scientific workflows in virtual reality

To advance research and improve the scientific return on data collection and interpretation efforts in the geosciences, we have developed methods of interactive visualization, with a special focus on immersive virtual reality (VR) environments. Earth sciences employ a strongly visual approach to the measurement and analysis of geologic data due to the spatial and temporal scales over which such data ranges. As observations and simulations increase in size and complexity, the Earth sciences are challenged to manage and interpret increasing amounts of data. Reaping the full intellectual benefits of immersive VR requires us to tailor exploratory approaches to scientific problems. These applications build on the visualization methods strengths, using both 3D perception and interaction with data and models, to take advantage of the skills and training of the geological scientists exploring their data in the VR environment. This interactive approach has enabled us to develop a suite of tools that are adaptable to a range of problems in the geosciences and beyond.

Overriding plate controls on subduction evolution

Geologic and geophysical observations indicate that the thickness, density, and strength of the lithosphere vary on the Earth. However, the role of the overriding plate lithosphere properties on the evolution and morphology of subduction is not well understood. This paper presents time-dependent numerical models of subduction that vary the overriding plate thickness, strength, and density and allow for a plate interface that evolves with time via an anisotropic brittle failure rheology. We examine the effect of these parameters on subduction evolution, in particular, the emergence of (a) asymmetric versus symmetric subduction, (b) trench retreat versus advance, (c) subduction zone geometry, (d) slab stagnation versus penetration into the lower mantle, and (e) flat slab subduction. Almost all of the models presented result in sustained asymmetric subduction from initiation. Trench advance occurs in models with a thick and or strong overriding plate. Slab dip, measured at a depth below the plate boundary interface, has a negative correlation with an increase in overriding plate thickness. Overriding plate thickness exerts a first-order control over slab penetration into the lower mantle, with penetration most commonly occurring in models with a thick overriding plate. Periods of flat slab subduction occur with thick, strong overriding plates producing strong plate boundary interface coupling. The results provide insight into how the overriding plate plays a role in establishing advancing and retreating subduction as well as providing an explanation for the variation of slab geometry across subduction zones on Earth, where similar patterns of evolution are observed.

Styles of rifting and fault spacing in numerical models of crustal extension

Extension of the Earths crust can result in differing styles of rifting, such as horst-and-graben, half-graben, metamorphic core complexes and areas of distributed crustal thinning. Faulting patterns can range from either being distributed to highly localized. Observations indicate that the factors controlling the extensional deformation, symmetry, and fault spacing include rheological aspects such as the yielding mechanism and strain softening, and physical aspects such as initial heterogeneities and the strength of the lower crust compared to the upper crust. Time-dependent numerical models of extension are presented, which investigate the influence of the yielding mechanism, lower crust strength, strain weakening, and initial heterogeneity in the crust have on (a) the style of rifting, (b) fault spacing, and (c) integrated strength in the upper crust. Models with an anisotropic yielding mechanism result in more realistic lithospheric strength profiles, slip plane angle distributions, and fault interaction than models with an isotropic yielding mechanism. Heterogeneity type and yielding mechanisms have the largest effect on the resulting symmetry of deformation, whereas the amount of strain weakening has the greatest influence on asymmetry. The likelihood of the metamorphic core complex mode occurring is primarily controlled by lower crust strength. Crustal thinning is encouraged by both low amounts of strain weakening and a strong lower crust. Lower crust viscosity exerts the primary control on whether the resulting deformation is distributed or localized. The degree of strain weakening has the largest influence on the average strength of the upper crust and the slope of the strength profile in the upper crust.

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.

The zone of influence of the subducting slab in the asthenospheric mantle

Due to the multi-disciplinary nature of combined geodynamics and shear wave splitting studies, there is still much to be understood in terms of isolating the contributions from mantle dynamics to the shear wave splitting signal, even in a two-dimensional (2D) mantle flow framework. This paper investigates the viscous flow, development of lattice preferred orientation (LPO), and predicted shear wave splitting for a suite of buoyancy driven subduction models to shed light on the nature of the slab driven asthenospheric flow field. The models show that the spatial dimensions of the slab driven zone of influence in the mantle, LPO fabric, and resulting synthetic splitting are sensitive to the slab dip and slab strength, with the maximum LPO fabric development occurring for models with a moderate initial slab dip and weaker slab. The models show that olivine LPO fabric in the asthenosphere generally increases in alignment strength with increased proximity to the slab, but can be transient and spatially variable on small length scales. The results suggest LPO formed during initial subduction may persist into the steady-state subduction regime. The vertical flow fields in the asthenosphere can produce shear-wave splitting variations with back-azimuth that deviate from the predictions of uniform trench-normal anisotropy, a result that bears on the interpretation of complexity in shear-wave splitting observed in real subduction zones. Furthermore, the models demonstrate that the corner flow paradigm should not be equated with a 2D subduction framework.

Distance field computation for geological slab surface data sets

The three-dimensional shapes of tectonic plates that sink into the Earth’s mantle (slabs) are the starting point for a range of geoscience studies, from determining the forces driving the motion of tectonic plates, to potential seismic and tsunami hazards, to the sources of magmas beneath active volcanos. For many of these applications finite element methods are used to model the deformation or fluid flow, and therefore the input model parameters, such as feature geometries, temperature or viscosity, must be defined with respect to a smooth, continuous distance field around the slab. In this paper we present a framework for processing sparse and noisy seismic data (earthquake locations), defining the shape of the slab and computing a continuous distance function on a mesh with variable node spacing. Due to the inhomogeneous volumetric distribution of earthquakes within the slab and significant inaccuracies in the locations of earthquakes occurring hundreds of kilometers below the Earth’s surface, the seismicity data set is extremely noisy and incomplete. Therefore, the preprocessing is the major part of the framework consisting of several steps including a point based smoothing procedure, a powerful method to use other observational constraints on slab location (e.g., seismic tomography or geologic history) to extend of the slab shape beyond earthquake data set and continuous resampling using moving least squares method. For the preprocessed point data we introduce approaches for finding the three-dimensional boundary of the slab and a subdivision of the slab into quadric implicit polynomials. The resulting distance field is then compiled from distances to the piecewise continuous approximation of the slab and distances to slab boundary.