2 7
Chapter Deformation in transcurrent and extensional environments with widely spaced weak zones1
2.1 Abstract
Previous mechanical models of the western U.S. have concluded that plate boundary forces cannot generate far-field deformation. Such models have ignored preexisting large-scale lithospheric strength variations, an assumption that appears to be inconsistent with seismically determined variations in lithospheric structure. We have formulated a three-dimensional viscous flow model with imposed plate motions, but include lateral zones of low viscosity. These models show that strain rates are concentrated in weak zones with adjacent blocks experiencing little deformation. Deformation can extend far inboard of plate boundaries, contrary to the result of previous studies with rheologically homogeneous plates, and apparently compatible with the variation in seismic velocity and GPS determined deformations in the western U.S. These results suggest that plate boundary forces cannot be neglected in the deformation of the western U.S., including the Cenozoic extension of the Basin and Range Province.
Published by Eunseo Choi and Michael Gurnis in Geophysical Researh Letter, 2003, 30(2), 1076, doi: 10.1029/2002GL016129.
1
8
2.2
Introduction
Plate boundary and buoyancy forces could be significant in the extensional deformation of the western U.S., but since the deformation of this region is probably the result of a dynamic balance among multiple forces, attempts to determine which force dominates have been inconclusive. Since continental crust is pervasively faulted, England and Jackson (1989) and subsequent workers have assumed that at sufficiently large scales, the crust can be approximated as a non-linear viscous fluid. Making this assumption, with a single power law rheology approximating the vertically averaged lithosphere, Sonder et al. (1986) showed that boundary velocity decayed rapidly with distance away from the margin and that only a narrow region near either a transform or divergent boundary could be deformed. Sonder and Jones (1999) reached the conclusion that buoyancy forces dominate the western U.S. force balance by first showing that buoyancy forces are sufficiently large to cause the deformation and then by arguing that the influence of plate boundary forces is negligible within the continental interior. Using a thin viscous sheet model with velocity boundary conditions representing Pacific-North American transform motion, they rejected plate boundary forces as a viable mechanism for Basin and Range Province extension and concluded that plate boundary forces could not have generated the inferred strain rates throughout the western U.S. However, these models ignored lateral heterogeneities in lithospheric strength. Since the mantle lid is the strongest part of the lithosphere, its thickness is closely related to the strength of the lithosphere. Seismic studies within the Basin and Range Province have shown marked differences in the thickness of the mantle lid (Fig.2.1) that should presumably be reflected in the strength of the lithosphere. For example, by compiling seismic refraction profiles, Braile et al. (1989) showed that there are lower seismic velocities within the upper mantle beneath the northern Basin and Range Province while seismic velocities under the Sierra Nevada and Colorado plateau are higher than average. York and Helmberger (1973) detected a band of low seismic velocities across the northern Basin and Range Province from the middle of
9
Figure 2.1: Shaded relief map of the northern Basin and Range Province, roughly showing the locations of low seismic velocity zones (York and Helmberger, 1973; Braile et al., 1989; Humphreys and Dueker, 1994b) and of the high strain rate region found with the BARGEN data (Niemi et al., 2004). It also shows a probable division of the northern Basin and Range Province into regions of relatively high strain rate (with brighter shades) and low strain rate (with darker shades).
10 Utah and traversing southern Nevada and southern California. This upper mantle structure is consistent in a more recent tomographic study by Humphreys and Dueker (1994a). This low velocity band may be more aptly described as a thin upper mantle lid beneath the Basin and Range-Colorado Plateau boundary, which was revealed by teleseismic P wave travel times (Zandt et al., 1995). Melbourne and Helmberger (2001) showed that lithospheric lid thickness varies significantly from 55 km along the coast of California to near zero thickness along eastern California. We reappraise the argument against plate boundary forces and show that previous studies may not have adequately considered an important aspect of regional structure. New models presented here include lateral variations in viscosity to reflect the observed seismic velocity variations as a proxy for lithospheric strength. With the new formulation, we argue that plate boundary forces could generate substantial deformation within continental interiors in both strike-slip and extensional environments and that plate boundary forces should still be considered as a viable component in the force balance of the western U.S. during Cenozoic time.
2.3
Model setup
A three-dimensional (3D) spherical finite element method (Zhong et al., 2000) is used to model the influence of a discontinuous viscosity distribution on lithospheric deformation. Using the Boussinesq approximation, the equations of mass and momentum conservation are: u=0 (2.1)
-P + 2 u + 0 g(T - T0 )e3 = 0,
(2.2)
where e 3 is a unit vector in the radial direction; u is the velocity; T and T0 are the temperature and the reference temperature, respectively; P is the pressure; is the viscosity; is the coefficient of thermal expansion; 0 is the mantle density; and g is
11 Table 2.1: Model parameters Symbol Name 0 0 c1 c2 T0 T1 Coefficient of thermal expansion Reference density Thermal diffusivity Reference viscosity Constant for viscosity Constant for viscosity Reference mantle temperature Surface temperature Value 210-5 K-1 3300 kg/m3 10-6 m2 /s 1021 Pas 9.50614 1.02126 1300 C 0 C
the acceleration of gravity. Viscosity is both temperature- and position-dependent: (T, , r) = C(, r) 0 exp c1 c1 - c2 + T /T0 c2 + 1.0 , (2.3)
where 0 is the viscosity when T =T0 (viscosity within the asthenosphere); and r are the polar and the radial coordinates, respectively; c1 and c2 are constants chosen so that the viscosity at T =0 C is higher than 0 by two orders of magnitude. The role of C(,r ) is the spatially varying part of the viscosity. Model parameters are given in Table 2.1. The spatial extent of the model domain is 60 in longitude, 40 in latitude, and 637 km in depth (Fig.2.2a). Only the inner region, 30 20 637 km, is of interest; the outer region used to reduce the influence of the vertical walls. Velocity boundary conditions enforce plate motion on the top within the inner region (within region A, Fig.2.2a). By creating a "transform boundary" along the equator with the length of the boundary being L, we can directly compare our result to earlier thin viscous sheet models with homogeneous rheologies through parameter conversion (Sonder et al., 1986). The temperature distribution is from a half-space cooling model defining a 100 km-thick thermal lithosphere (Fig.2.2b) and is horizontally uniform. The coefficient, C(,r ), is dependent on the polar angle, (0-180 ), allowing weak zones through the lithosphere at different distances from the plate boundary (Fig.2.2c). The weak zone in spherical shell models is defined as a narrow region of lower
12 viscosity compared to adjacent regions. The ratio of viscosity of the weak zone to that of the adjacent normal region at a certain depth is denoted as weak /normal . The zone of reduced viscosity extends to the base of the lithosphere, is 210 km in width, and is of the same length as the longitudinal dimension of the outer region (~ 6672 km, Fig.2.2a). There are two controlling parameters for the weak zone: the viscosity ratio and its distance from the plate margin. weak /normal varies from 10-1 to 10-4 . Three locations for the weak zone are considered as in the channel flow models: on the transform boundary (Fig.2.3a,b), far inboard of the margins (about 0.6L), meant to represent the miogeoclinal hinge of the North American Plate (Fig.2.3c,d). Noting that a natural weak zone is developed at the plate boundary in thin viscous sheet models, we also consider models with two weak zones: one at the plate boundary and the other at 0.6L (Fig.2.3e,f).
2.4
Results
The results from the spherical shell models show a consistent dependence of deformation on viscosity reduction and weak zone position. We first consider the influence of the weak zones in spherical shell models in which plate motion along the boundary is strike-slip. We show the component (v) of velocity parallel to the strike-slip margin (see Fig.2.2a for a coordinate layout), normalized by the boundary velocity (v0 ), versus normalized distance (x/L) from the plate boundary along the 30 meridian. The gradient of velocity is particularly steep within the weak zones while it is mild outside of the weak zones, especially for models with the weak zone on the boundary between the two plates. In other words, there are mildly deforming blocks separated by shear zones, a conclusion which is not particularly sensitive to the magnitude of the viscosity reduction once the velocity has been reduced by one order of magnitude (Fig.2.3a,c). There has also been a normal component to the motion of the Pacific plate outboard of the northern Basin and Range Province (e.g., Atwater and Stock, 1998) and so a second boundary condition is considered, which can also be compared with earlier
13
Figure 2.2: (a) Coordinates of the finite element mesh and configuration for boundary conditions. The inner region is the volume outlined by thick solid and dashed lines. x-axis is along the 30 meridian and y-axis along the equator. u and v are velocity tangential to x- and y-axis, respectively. Velocity boundary conditions are applied on the surface in the regions shaded with dark gray: u=0.0, v=v0 =5.0 cm/yr in region A, while u = v = 0 in B. region v=0.0 and u/x=0.0 on the vertical boundaries. The mesh has 32 elements in longitude and 64 in latitude (with regular spacing) and 16 in depth (with higher spatial resolution through the upper 100 km) (b) Variations of temperature with depth. Temperature is horizontally uniform. (c) Viscosity variation with depth. Viscosity for the blocks is drawn in a solid line. Viscosity for a weak zone with viscosity ratio weak /block of 10-2 is drawn in a dashed line.
14 models (Figure 3b,d). With purely normal velocity boundary conditions (u 0 = -5 cm/yr and v = 0.0) applied on the region A (Fig.2.2a), the normal component of velocities for each configuration of weak zone shows a similar pattern compared with the transform boundary condition: strain rate concentration in a weak zone between mildly deforming blocks. It is noteworthy that the behavior of these numerical models is in sharp contrast to the analytical solutions for a semi-infinite, thin viscous sheet model for a lithosphere with a non-linear rheology. For example, in the weak-zone model with the weak zone placed far inboard of the margin (at x/L=0.6 Fig.2.3c), the velocity decreases to 50 % of the boundary value at x/L=0.3, while in the thin sheet model with n > 1 the corresponding velocity is essentially zero. The characteristic length scale of velocity decay is thus larger in these numerical models then in the thin viscous sheet models. It is even larger for the normal velocity boundary condition. The decrease of velocity is only about 10 % at x/L =0.3 for the model with the weak zone far inboard of the boundary (Fig.2.3d). In thin viscous sheet models, as the non-linearity of the viscosity increases (i.e. as the exponent, n, increases), the characteristic length scale of the decay of velocity decreases (England et al., 1985; Sonder et al., 1986). However, lowering the viscosity within the weak zone leads to an increase in characteristic length scale over which the deformation is distributed. Models with two weak zones (Fig.2.3e,f) show the importance of incorporating pre-existing strength variations in the lithosphere even though a non-Newtonian rheology may result in a weak plate boundary. A weak plate boundary zone, simulating the effect of a power-law rheology, does take up a large amount of deformation, but the interior weak zone still deforms faster than the strong interior region. This deformation pattern becomes obvious in models with a relatively strong strike-slip boundary (Fig.2.3e), or with imposed normal velocity boundary conditions (Fig.2.3f). Thus, our results are relevant to the deformation of the western U.S., because this region has a significant component of extension.
15
Figure 2.3: Velocity profiles for transform boundary conditions and for normal velocity boundary conditions. (left) y-velocity, v, plotted as a function of distance from the transform boundary, x. v is normalized with respect to v0 =5 cm/yr, the boundary velocity, and x with respect to L. (right) Normalized x velocity, u/u0 , versus normalized distance from the plate boundary, x/L, for the normal velocity boundary conditions: u = u0 =-5 cm/yr and v=0.0. Different symbols are used for numerical solutions with different viscosity ratio of the weak zone to the normal value (w /b ). For comparison, analytical solutions for thin viscous sheet model (TVS) with the corresponding velocity boundary condition (England et al., 1985) are plotted with different types of lines. (a,b) The weak zone is on the boundary between the two plates, (c,d) at x/L =0.6, far inboard of the plate. (e,f) Models with two weak zones: the outer one at the plate boundary and the interior one at x/a =0.6, of which viscosities are o and i , respectively. While i /b fixed at 10-2 , o /i changed from 1 to 102 .
16
2.5
Discussion
In order to verify the importance of either plate boundary forces or buoyancy forces, heterogeneities in the lithospheric strength need to be properly incorporated into deformation models. Variations in lithospheric strength are generally critical for understanding the dynamics of lithospheric deformation, but studies such as Neil and Houseman (1997) and Ellis et al. (1998) have only applied such complications to compressional settings. As for the Basin and Range Province, Sonder and Jones (1999) treated the Sierra-Nevada and Great Valley block as a relatively rigid area in a thin-sheet model, showed that deformation remained concentrated near the transform boundary, and argued that the inclusion of strong regions did not alter the decay of the plate boundary forces into continental interiors. However, their specific lithospheric strength structures did not incorporated all of the crustal and upper mantle features suggested by seismic studies. The eastern part of the northern Basin and Range Province appears to be different in terms of its seismically determined crust and upper mantle compared to adjacent regions. The eastern Basin and Range Province (at 39 -40 N) with a width of 300-km and located west of the Colorado Plateau coincides with the northern end of a low velocity zone (York and Helmberger, 1973; Humphreys and Dueker, 1994b) (Fig.2.1). The association of this seismic low velocity region with high heat flow (120 mW/m2 ) and recent active volcanism (along the St. George volcanic trend) led Humphreys and Dueker (1994a) to suggest the existence of partial melts in the upper mantle. They thought about 1-3 % of partial melt would be necessary to explain the 4 % of P wave velocity reduction. Since partial melts generally reduce viscosity, the heat flow and seismic studies support the presence of a weak zone beneath the eastern part of the northern Basin and Range Province. This putative weak zone is geometrically similar to our models with the weak zone far inboard of plate boundary. Consequently, based on our models (Fig.2.3) we would expect that this region to be straining much faster compared to the Basin and Range Province further to the west or the Colorado Plateau further to the east.
17 This interpretation is consistent with recent geodetically determined strain rate distributions. Wit...