High-resolution lithosphere viscosity and dynamics revealed by magnetotelluric imaging

See allHide authors and affiliations

Science  30 Sep 2016:
Vol. 353, Issue 6307, pp. 1515-1519
DOI: 10.1126/science.aaf6542


An accurate viscosity structure is critical to truthfully modeling lithosphere dynamics. Here, we report an attempt to infer the effective lithospheric viscosity from a high-resolution magnetotelluric (MT) survey across the western United States. The high sensitivity of MT fields to the presence of electrically conductive fluids makes it a promising proxy for determining mechanical strength variations throughout the lithosphere. We demonstrate how a viscosity structure, approximated from electrical resistivity, results in a geodynamic model that successfully predicts short-wavelength surface topography, lithospheric deformation, and mantle upwelling beneath recent volcanism. We further show that this viscosity is physically consistent with and better constrained than that derived from laboratory-based rheology. We conclude that MT imaging provides a practical observational constraint for quantifying the dynamic evolution of the continental lithosphere.

During the last decade, progress in both seismic imaging (1, 2) and geodynamic modeling (3, 4) produced many useful insights into the large-scale dynamics of the continental lithosphere. However, the driving forces for fine-scale (<200-km) tectonic processes, such as those within the western United States, remain heavily debated (510). Large-scale deformation generally results from mantle convection-induced sublithospheric stress, such as dynamic topography and mantle traction (3, 4, 7, 1012). Consequently, the existing debate on the fine-scale tectonic deformation is largely due to the uncertain lithospheric buoyancy and viscosity structures. Fortunately, recent geophysical measurements constrain the density distribution of the crust and mantle lithosphere very well (13, 14). However, the detailed viscosity structure of the lithosphere still remains poorly understood.

Taking the western United States as an example, its lithospheric density structure could be derived from multiple geophysical observations, including seismic velocity and heat flow measurements (13), as well as the geometry of lithospheric discontinuities (14). In contrast, the spatial pattern of lithospheric viscosity is more difficult to infer observationally because of the intimate involvement of time. Previous attempts included matching geodetic deformations with gravitational potential energy (GPE), by assuming either isostatic equilibrium (9) or considering large-scale mantle convection (10), and fitting the postseismic relaxation of earthquakes (15). However, considerable discrepancies still exist among these inferred lithosphere viscosity structures.

Magnetotelluric (MT) imaging represents a promising and independent approach to inferring a detailed lithospheric viscosity structure. This approach builds on direct geophysical and geological observations, requiring no a priori assumptions about lithospheric dynamics. Changes in the physical state that alter electrical conductivity generally affect viscosity as well, which allows us to use MT sounding to investigate the spatial distribution of viscosity. In reality, the viscosity depends on temperature, strain rate, grain size, and composition (1621). Among these, compositional effects, specifically trace amounts of water and melt, are the most difficult to determine, but play a crucial role in the reduction of strength (17). Fluids, including melts, accumulate in regions of strong deformation (17). The presence of melts greatly reduces rock strength while enhancing electrical conductivity (18). Hydration of nominally anhydrous minerals also increases electrical conductivity [(22) and references therein]. Likewise, laboratory experiments on major minerals (e.g., olivine, pyroxene, and garnet) indicate weakening in the presence of water (1921).

We established a theoretical framework between electrical conductivity and viscosity based on their similar controlling factors. Effective viscosity is the ratio of stress to strain rate, Embedded Image, where σ is stress and Embedded Image is the strain rate. A general form of this relation isEmbedded Image (1)where d is grain size; COH, water fugacity; E* and V*, activation energy and activation volume, respectively; P and T, pressure and temperature, respectively; R, the ideal gas constant; and A, m, n, and r, all laboratory-derived parameters (23). The electrical resistivity model takes a similar formEmbedded Image (2)where the parameters are also like those in Eq. 1 (24).

Using Eq. 1, one could first estimate stress (σ), independent of viscosity, from local buoyancy variations (38) and/or far field mechanical conditions (9, 10, 15). One could then numerically calculate the accurate stress, including that induced by lithosphere deformation for a given viscosity structure. In contrast to the viscosity formulation, the omission of grain size in the formulation for resistivity (Eq. 2) could likely be naturally restored in the structural information of MT sounding (Fig. 1), as demonstrated below. Consequently, we suggest, using the similarity between Eqs. 1 and 2, that electrical resistivity is positively correlated with effective viscosity (ρ ∝ η). Previous conceptual models of rifting that highlight variations in strain appear to mimic patterns seen in the midcrustal electrical conductivity under similar circumstances (25). We propose that these high-conductivity patterns (Fig. 1, A and B) also represent low-viscosity regions owing to consequences of past deformation of the continental lithosphere.

Fig. 1 Geophysical characteristics of the Basin and Range and Colorado Plateau.

(A) Shaded topography, surface heat flow interpolated from the global heat flow database (13), crustal thickness estimated by seismic receiver functions (14), MT station locations, and late Cenozoic volcanism. (B) Inverted resistivity model from (25), along the 2D profile shown in (A). The vertical gray bars at the surface represent MT stations. Note the correlation of faults and weak zones with low resistivity and that of cratonic lithosphere with high resistivity.

We think that the above theory and observations support using MT images as a proxy for mechanical strength throughout the lithosphere. However, we recognize several difficulties limiting the direct application of the above relations in accurately estimating viscosity. First, large uncertainties in laboratory models result in large uncertainties in viscosity. Second, developing a rigorous model is restricted by additional uncertain factors, such as composition and grain size. Third, MT models do not yield precise estimates of the absolute conductivity, as they are sensitive to conductance (the product of conductivity and thickness) and are insensitive to the resistivity of resistors. Furthermore, the spatial sensitivity of magnetotellurics to electrical conductivity variations with depth and the smoothing effect due to inversion cannot capture the true scale of viscosity variations. With these uncertainties in mind, we processed the MT image in a way similar to that in which seismic tomography is converted into density anomalies (2, 3). We converted the resistivity structure into a viscosity profile using an empirical relation,Embedded Image (3)where η is effective viscosity; η0, reference viscosity (1020 Pa·∙s); ρ, apparent resistivity; ρ0, reference resistivity (102 Ω m); and C0 and C1, numerical coefficients. Because of the unknown values of C0 and C1, both the regional average and maximum contrast of the viscosity structure, controlled by these two respective coefficients, are uncertain. Subsequently, we updated these coefficients and, thus, viscosity properties using observational constraints, including surface topography and intraplate deformation. Overall, we looked for an effective lithospheric viscosity structure resembling the spatial pattern of the MT image that satisfies the available constraints. For simplicity and to show the effectiveness of the technique, we limited ourselves to a two-variable problem where C0 and C1 are constant over space.

We used a recent high-resolution MT inversion for the tectonically active western United States (25) and estimated the lithospheric buoyancy by considering thermal and compositional effects. The thermal structure was estimated from surface heat flow by solving for a two-dimensional (2D) steady-state thermal diffusion problem with the observed heat flow as the upper boundary condition (Fig. 2 and fig. S1). We assumed that the 1350°C isotherm approximately follows the 100 Ω m contour of the apparent resistivity (Fig. 1B) (26). The compositional buoyancy of the continental crust was approximated using the recent estimate on Moho depth (14), highlighting thicker crust beneath the Colorado Plateau (CP) than the Eastern Great Basin (EGB) and the Transition Zone (TZ) (Fig. 2). A uniform crustal density of 2850 kg·∙m−3 was adopted for the study region (27).

Fig. 2 Models with a 1D viscosity profile predicting intraplate deformation, surface topography, and mantle flow.

(A and B) Intraplate deformation and surface topography. The black curves in (A) and (B) represent the observed Basin and Range extension rate and surface topography, respectively. The green and red colors represent results from a weak and strong lithosphere, respectively. (B) The corresponding contributions to topography from crustal buoyancy (orange), thermal buoyancy (blue), and their combined effects (red) are shown. (C) Mantle flow. The background viscosity is for the strong lithosphere model, but the mantle flow from both the weak (green) and strong (red) lithosphere models are shown. Gray lines represent the geotherms, and magenta line marks the Moho.

We verified the resulting lithospheric density structure with both analytical and numerical solutions (Fig. 2B). By assuming a simple viscosity profile with only depth-dependence, similar to traditional geodynamic models (24), we found that the topography contribution owing to crustal buoyancy (orange line in Fig. 2B) is a mirror image of the Moho (magenta line in Fig. 2C), consistent with Airy compensation. The resulting topographic relief owing to thermal buoyancy (blue line in Fig. 2B) also resembled that from a recent global geodynamic model (4), a resemblance that validated our calculation. More important, the predicted total topography also matched a smoothed version of the observed topography. We note that the regional model used here could not properly reproduce the baseline topography of the region because of the undefined sea level.

This model (Fig. 2), although it satisfied the smooth profile of surface topography, failed to reproduce the westward block motion (3 mm/year) of the EGB relative to the CP (Fig. 2A), as revealed geodetically (28) (fig. S2). If we assume an average lithospheric viscosity of 1022 Pa·∙s over a lithospheric thickness of 50 km, as suggested by earlier studies (9), the model predicts negligible lithospheric deformation (Fig. 2A). A weaker lithosphere (with a viscosity of 1021 Pa·∙s) leads to larger intraplate deformation, but the resulting velocity profile differs from that observed, with most deformation occurring inside the CP. Physically, the deformation mechanism here is similar to earlier GPE models, so that the mismatch suggests that either some far-field tectonic stress is needed (9, 10) or the viscosity structure is inappropriate or possibly both. Furthermore, these models predict limited magnitudes of mantle upwelling flow beneath the region, inconsistent with the recent volcanic activities within the TZ.

Adopting a spatially variable viscosity structure, converted from the MT image using Eq. 3 (Fig. 3), greatly improved the prediction of surface topography and lithosphere deformation. To show the validity of this approach on inferring viscosity by matching tectonic observables, we performed a parameter test by systematically changing the values of C0 and C1, which control the lithosphere’s average strength and spatial viscosity variation, respectively. Among the five models showing different amounts of spatial viscosity contrast (Fig. 3A and fig. S3), the smallest amount (10×) of viscosity variation generated surface predictions similar to those in the 1D viscosity model (Figs. 2 and 3), with strong lithospheric deformation inside the CP and the TZ. The deformation localized more toward the EGB and TZ with increasing viscosity variations. In addition, a perfect step-profile of intraplate velocity emerged as the viscosity contrast increases to 107, similar to the observed lithospheric deformation (28). A larger viscosity contrast (109) produced progressively decreasing surface velocity from the EGB to the east, which was inconsistent with observations.

Fig. 3 Sensitivity of tectonic constraints on the MT-converted lithosphere viscosity.

(A and B) Dependence of (A) intraplate deformation rate and (B) surface topography on the amount of spatial variation of lithosphere viscosity. (C) The background color shows the viscosity with a 10-fold spatial variation, and the three mantle flow fields correspond to cases with 10 (blue), 105 (orange), and 109 (green) of viscosity variations, respectively. (D and E) Dependence of tectonic constraints on the average (or net) viscosity of the lithosphere. (F) Background color represents viscosity of the strongest (C0 = 8.0) lithosphere case, and the three flow fields are from cases with C0 of 8.0 (blue), 2.0 (orange), and 0.5 (green), respectively. The red bars on the surface mark the MT stations.

Besides intraplate deformation, the predicted surface topography also showed clear improvements relative to the 1D viscosity model. The MT-based viscosity structures predicted large variations of surface topography on length scales as small as ~20 km (fig. S3). We further corrected the flexural effect on these raw topography signals (26) while considering the spatially varying lithospheric elastic thickness (29) (fig. S4). In general, predicted topography smooths with the degree of viscosity homogeneity (Fig. 3B). A 10-fold lateral variation in viscosity results in topography that (blue line in Fig. 3B) largely resembled the 1D viscosity model (pink line in Fig. 3B), but a clear difference was that the former predicted more uplift (subsidence) in the TZ (CP) because of a larger contribution from subsurface convection (Fig. 3C). This topographic difference became more pronounced with increasing viscosity contrast, as shown by the increasing amplitude of topographic roughness. We suggest that the short-wavelength topography reflects shear displacement of crustal blocks along the resolved faults, driven by lithospheric strain variation resulting from flow at greater depths. This relation can be seen from the striking correlation of localized elevation peaks (Fig. 3, B and E) with the strong crustal blocks right below (Fig. 3, C and F), as well as the independence of topography on local lithospheric density structure (fig. S5). The model that best matched lithosphere deformation (Fig. 3A, with viscosity variations of 105 to 107) also well predicted the fine-scale topography, with remarkable similarity in both wavelength and amplitude (Fig. 3B).

Similar tests on the average (or net) lithospheric viscosity by varying the parameter C0 demonstrated an almost linear correlation with the magnitude of intraplate deformation (Fig. 3D) and topography variation (Fig. 3E and fig. S6) (26), which implied a strong constraint. Collectively, we found that intraplate deformation and surface topography represented sufficient and complementary constraints on the average strength and spatial variation of lithosphere viscosity structure (Fig. 3). The model that best matched all observations had a spatial viscosity variation between 3 × 1023 Pa·∙s and 3 × 1017 Pa·∙s (Fig. 4A).

Fig. 4 Effective viscosity structures derived from the MT image and laboratory-based rheologies.

(A) Best-fitting model with an MT-converted viscosity structure, including six orders of magnitude of viscosity variations. (B) Viscosity structure and model predictions using a power-law rheology (n = 3 in Eq. 1 and other parameters listed in supplementary materials) and pseudoplasticity (yield stresses of 40,150 MPa for above and below 40 km in depth, respectively). The general similarity between the two viscosity structures and model predictions confirms the physical validity of the MT-converted viscosity.

Another robust prediction from our models was a region with vigorous mantle upwelling at the lower-crust and uppermost mantle depths beneath the TZ, where the peak ascending velocity exceeded 5 cm/year (Figs. 3 and 4). The locally reduced viscosity and the large lateral gradient of buoyancy in this region facilitated the mantle upwelling responsible for late Cenozoic volcanism, high surface heat flow (Fig. 1), and rough TZ topography (Figs. 3 and 4). The apparent eastward encroachment of flow beneath the western edge of the CP was consistent with the volcanic history (5) and the proposed destabilization of the plateau’s cratonic root (8).

In our best-fitting model (Fig. 4A), the dramatic lithospheric viscosity reduction to the west of the CP led to a sharp change of lithospheric deformation rate, where almost nonexistent internal deformation inside the CP increased rapidly to ~3 mm/year westward within a short distance of <100 km, which was close to that observed (28). Toward the western end of the MT survey line, the predicted deformation rate dropped slightly to ~2.5 mm/year. This decrease was likely related to an apparent lithospheric strengthening because of the limited resolution of the MT image toward the edge of the inversion domain. The overpredicted topography within the western CP may reflect surface erosion attributable to the local river system (Fig. 4A). Alternatively, this may be due to a lithosphere density anomaly not considered in our simple buoyancy structure or local inaccuracy of the MT-image in representing lithosphere viscosity.

We further investigated the physical consistency of the MT-converted viscosity with that forwardly derived from laboratory-based rheologies (23). By solving the mantle flow using the power-law rheology with other factors in Eq. 1 assumed uniform (26) (fig. S7), we obtained a viscosity structure that resembled the MT image at large scales, but the model failed to predict both fine-scale topography and lithosphere deformation (fig. S7). Inclusion of pseudoplasticity in the model improved the fit in both viscosity and lithosphere deformation (Fig. 4B). However, deformation within the CP was overpredicted, likely related to the artificial return flow close to the eastern boundary of the model; this model also lacked localized crustal weak zones and, thus, the short-wavelength topography within the EGB and TZ. These mismatches should reflect the missing effects of compositional variation and grain size in the viscosity calculation.

The general consistency between the forward [using rheological laws (Fig. 4B and fig. S7)] and their inverse [using MT (Fig. 4A)] viscosity calculations showed that the latter properly captured the effects of strain-rate dependence and plastic deformation. The fact that the MT-converted viscosity best matched observations implied that the MT image also captured the effects of other viscosity-influencing factors like composition and grain size, which are difficult to infer using a forward approach. We suggest that future research including MT as a quantitative constraint will better unearth the physics of lithospheric rheology. Consideration of mineral physics and geochemical and geological data in high-fidelity geodynamic models will help to push this research frontier forward. Application of the MT-based modeling approach to other tectonic regions (30) should help to better understand the detailed dynamic evolution of continents including topography, deformation, and volcanic history.

Supplementary Materials

Materials and Methods

Supplementary Text

Figs. S1 to S7

References (3134)

References and Notes

  1. Methods and data processing techniques are available in the supplementary materials on Science Online.
  2. Acknowledgments: The authors thank J. Hu for helping with the non-Newtonian rheology. L.L. is supported by NSF grants EAR-1554554 and ACI-1516586. Collection of MT data and support for D.H. are provided by NSF grants to P. E. Wannamaker EAR81-16602, 84-17765, and 02-30027. L.L. conceived the project and performed the geodynamic calculations. D.H. provided expertise on MT-sounding and lithosphere buoyancy. Both authors contributed to manuscript preparation. MT data are available in reference (25). Additional information on the geodynamic modeling is available in the supplementary materials.
View Abstract

Stay Connected to Science

Navigate This Article