Rheological dependence of extension in wedge models of convergent orogens

ELSEVIER

Tectonophysics 305 (1999) 419–435

Rheological dependence of extension in wedge models of convergent orogens Sean D. Willett * Department of Geosciences, Pennsylvania State University, University Park, PA, 16802, USA Received 14 October 1997; accepted 19 January 1999

Abstract Although several mechanical models of extension in contractional orogens have been proposed, in many cases these models have been based on viscous constitutive laws. Lack of more complicated rheological models has primarily been due to the technical difficulty of solving for deformation with complex rheological behavior and boundary conditions. A finite-element model is presented that allows for solution of the wedge mechanical problem with viscous, non-linear viscous and Coulomb plastic rheologies. This model is used to investigate extension in a doubly-vergent wedge model in which deformation is primarily contractional, driven by a reversal of basal traction. Proposed models showing extension in the crust at shallow crustal levels, in spite of overall contraction and contraction at depth, are confirmed, but only for high Argand number (low crustal viscosity). At low Argand number little or no extension is observed. Non-linear stress dependence (power-law flow) also suppresses extension, but not entirely; some extension is observed at high effective Argand number. However with a Coulomb plastic rheology, no extension is observed for a wide range of parameter values. No extension is observed for models with viscous deformation at depth and plastic deformation near the surface. This suggests that the extension is a function of the rheologic model, not the effective strength of the crust. Coeval extension and contraction in convergent orogens with steady external forcing therefore requires viscous deformation to very shallow depths.  1999 Elsevier Science B.V. All rights reserved. Keywords: orogen; extension; critical wedge; finite element

1. Introduction Observations of structures, such as normal-sense shear-zones, that indicate near-horizontal extension of upper crustal rocks within otherwise contractional orogenic belts are common (Platt, 1986; Dewey, 1988; Ratschbacher et al., 1989; Hodges et al., 1992; Malavielle, 1993; Crespi et al., 1996). In many cases these structures indicate extension parallel to the Ł Present

address: Department of Geological Sciences, University of Washington, Seattle, WA 98195.

strike of the orogenic belt and as such, could represent gravitationally-driven extension away from a region of thickened crust, or may be kinematic features of lateral extrusion oblique or normal to the direction of convergence. However, in many instances, extension is perpendicular to strike indicating extension in the direction of relative convergence. Timing control when available, often indicates that this extension and contraction are simultaneous. For example, in the Betic Cordillera of Spain, Platt and Vissers (1989) and Vissers et al. (1995) argue that extensional and contractional structures are synchronously

0040-1951/99/$ – see front matter  1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 0 - 1 9 5 1 ( 9 9 ) 0 0 0 3 4 - 7

420

S.D. Willett / Tectonophysics 305 (1999) 419–435

active and are kinematically linked. In the Himalaya, motion on the Main Central thrust and normal-sense motion on the South Tibetan detachment appear to have been concurrent (Burchfiel and Royden, 1985; Hodges et al., 1992). A potential modern example of this process may be provided in the arc–continent collision of Taiwan where geodetic data (Yu et al., 1995) and mapped normal faults (Crespi et al., 1996) indicate extension at high elevation in the interior of the otherwise convergent orogen. Several mechanisms have been proposed to explain these extensional features. Burchfiel and Royden (1985) proposed that the normal-sense displacement on the South Tibetan shear zone was the result of extrusion of the Himalayan crystalline rocks accompanied by thrusting along the Main Central thrust with a geometry similar to that shown in Fig. 1a. Analog models (Chemenda et al., 1995, 1996) have supported this theory for the Himalaya as well as providing a possible explanation for exhumation of high pressure rocks in subduction zones. Mechanical explanations of these models include strong buoyancy forces due to low density crust subducted into the mantle, or pre-existing weaknesses in the crust to allow normal-sense displacement in a horizontally compressive stress regime. In other cases where extension is spatially associated with subduction, it has been attributed to retreat or rollback of the subducting slab away from the interior of the overriding plate (Royden, 1993; Doglioni, 1994; Waschbusch and Beaumont, 1996) (Fig. 1b). Proposed examples of this process include the Apennines (Royden, 1993) and the Aegean (Dewey and Sengor, 1979; Angelier et al., 1982). Post-orogenic collapse of crust thickened by contraction is frequently cited as the process leading to extension of convergent orogens (Dewey, 1988;

Fig. 1. Models for extension in convergent settings. Double arrow indicates region of crustal extension. (a) Vertical block escape model of Chemenda et al. (1995). (b) Slab retreat. (c) Post-orogenic collapse. (d) Underplating driven extension in a viscous wedge (Platt, 1986). (e) Lower crustal extrusion from beneath an orogenic plateau (Royden, 1996). (f) Traction reversal model of Buck and Sokoutis (1994). Note that (d) through (f) assume viscous flow of the crust. See text for detailed descriptions of models. This paper investigates model (f).

S.D. Willett / Tectonophysics 305 (1999) 419–435

Brown and Carr, 1990; Malavielle, 1993) (Fig. 1c). This mechanism is poorly defined and the term has likely been overused. This mechanism should be distinguishable from the other mechanisms discussed here in that extensional deformation would post-date convergence, rather than being coeval. However, in practice the distinction is difficult as extensional collapse may drive contraction on the periphery of an orogen after the cessation of relative plate convergence. Also, if gravitational collapse is triggered by an external force, it might precede the cessation of contraction. An example of external forcing is collapse driven by uplift due to convective removal of the thickened lithospheric mantle (Houseman et al., 1981; Molnar et al., 1993). This mechanism has been proposed as being responsible for much of the east– west extension of the Tibetan Plateau (Mercier et al., 1987). There are also models for extension which do not rely on external forcing from mantle processes, but rather achieve extension through crustal mechanical processes. Platt (1986) proposed that extension at shallow levels in an accretionary wedge is driven by tectonic underplating, or transfer of accreted mass across a deep decollement at the base of the deforming wedge (Fig. 1d). Platt (1986) demonstrated that the stability of a wedge with a linear viscous constitutive law could be driven unstable by such underplating. A model recently proposed by Royden (1996) investigated viscous flow of continental crust thickened by orogeny to the scale of the Tibetan Plateau. With deformation through the entire crust and crustal thickness up to twice typical values, viscous flow takes on a complex, two- or three-dimensional pattern. For most thermal conditions and crustal viscosities, Royden found that crustal flow became localized within a lower crustal channel. Flow in this channel was always outward away from the highest topography. Under certain conditions, this outward flow became upwardly directed leading to localized surface extension near the edges of the orogenic plateau (Royden, 1996; Royden et al., 1997) (Fig. 1e). Buck and Sokoutis (1994) demonstrated through the use of both an analogue and a numerical model that a linear-viscous orogenic wedge can exhibit surface extension simply as the result of the convergent basal boundary motions. Convergent mountain belts

421

typically exhibit a large-scale structure of outwardly vergent structures. Insofar as structural vergence is the result of underthrusting and basal drag, this vergence reversal implies a reversal in the sense of basal traction across the orogenic belt (Willett et al., 1993). Buck and Sokoutis (1994) demonstrated that this basal traction reversal induces surface extension in the vicinity of the boundary velocity change (Fig. 1f). These last three models have in common the fact that extension is occurring during convergence. That is, net convergence is continuing between the two ‘plates’ involved in collision, but extension is still occurring in some spatially restricted region of the orogen. In addition, each of these models (Platt (1986), Buck and Sokoutis (1994), Royden (1996)) includes a crustal rheology which is linear viscous. The argument can be made that the models are more general and would hold for more complex rheologies, but it is not clear that extensional characteristics would persist. For example, Buck and Sokoutis (1994) noted that nonlinear viscosity decreased extension in their numerical model so that, insofar as a power-law viscosity with a stress exponent approaching infinity approximates a plastic material, no extension would be present with a plastic constitutive law. Rheological complexities thought to be characteristic of the crust include (1) temperature dependent viscosity (included in the model of Royden (1996)), (2) nonlinearity in the stress dependence of the flow such as exhibited by power-law flow, and (3) perhaps most importantly, low-temperature frictional plastic or brittle deformation of the shallow crust. The mechanical behavior of plastic, frictional wedges has been studied in detail (Chapple, 1978; Davis et al., 1983; Dahlen, 1984; Dahlen, 1990) and the resulting critical Coulomb wedge theory implies a mechanical response quite different from that of viscous wedges. In particular, it is much more difficult for steady or non-steady tectonic mechanisms to drive Coulomb wedges into a state of horizontal extension. According to critical wedge theory, a slight increase in taper angle increases the normal stresses within the wedge causing it to become stable, not extensional. In order for a wedge to extend, the taper angle must increase sufficiently to drive the wedge beyond this stable field and into collapse (Dahlen, 1984; Willett, 1992). An underplating mechanism

422

S.D. Willett / Tectonophysics 305 (1999) 419–435

(Fig. 1d) would thus be required to be efficient enough to increase the wedge taper angle far beyond that required for extension of a viscous wedge. The effects of a basal traction reversal are more difficult to ascertain, but sandbox experiments (Malavieille, 1984; Wang and Davis, 1996) have not documented extension for similar boundary conditions. These rheological complexities can also be addressed through the use of numerical models. In this paper I present a finite element model of crustal orogenic wedge formation with a simple mechanical parameterization. I use this model to test the rheological dependence of the traction-reversal model (Fig. 1f). The analysis is limited to this one model in order to present a suite of models exploring rheological variations, but many of the conclusions may apply to other models such as the Platt (1986) and Royden (1996) model. Models are also limited to plane-strain wedge conditions and so do not address the class of models which include out of plane motion.

2. Numerical model In order to test models with the range of rheological properties and boundary conditions proposed above, the versatility afforded by finite element techniques proves useful. The large deformation inherent to geologic problems precludes the straightforward use of more common elastic–plastic formulations in which an elastic solution is modified by spatially limited plastic deformation. Methods that work with a spatial grid and an Eulerian statement of equilibrium avoid problems of large mesh deformation and have proven useful in modeling orogenic deformation (Villote et al., 1982; Willett, 1992; Beaumont and Quinlan, 1994; Beaumont et al., 1995; Willett and Beaumont, 1994; Fullsack, 1995). Two dimensional force equilibrium can be expressed as a quasi-static equation of motion: @¦i j C ²g j ; j D 1; 2 @ xi

(1)

where summation is implied over repeated indices, ¦ is the 2-d stress tensor in spatial dimensions x1 and x2 , ² is the density and g j is the acceleration of gravity in the jth direction. If we assume a viscous

material, the constitutive law relating stress to rate of deformation is: ¦i j D

¦ Ži j C 2¼"i j

(2)

where "i j is the rate-of-deformation tensor, defined by the spatial velocity derivatives:   1 @vi @v j C "i j D 2 @x j @ xi ¼ is the viscosity, ¦ is the mean stress or pressure defined as: ¦kk ¦ D 3 and Ži j is the Kronecker delta, defined to be 1 when i D j and 0 otherwise. Rate-of-deformation differs from more common definitions of strain rate, but as the appropriate measure of the Eulerian rate of strain or deformation, it will be used interchangeably with the term strain rate in this paper. Combining Eq. 1 and Eq. 2 gives the equation of motion for a slowly-deforming (quasi-static) viscous material:   @ @vi @v j @¦ C¼ C C ²g j D 0 j D 1; 2 Žx j @ xi Žx j Žxi (3) with the 2 components of velocity, vi , i D 1; 2 and pressure, ¦ as dependent variables. The third equation to solve this system comes from conservation of mass, which for an incompressible material can be written in terms of the divergence of the velocity field: @vk D0 (4) @ xk Eqs. 3 and 4 make up a set of 3 equations and 3 unknowns which can be solved given appropriate boundary conditions. However, in applications to solid earth deformation, the appropriate constitutive law is frequently not linear viscous as suggested by the constant, scalar value of ¼ in Eq. 2. Many high temperature deformational mechanisms determined experimentally suggest that strain rate depends on stress to the power n giving a constitutive law of the form:   Q n (5) "P D A¦ exp RT

S.D. Willett / Tectonophysics 305 (1999) 419–435

where "P is the strain rate, A and n are constants, ¦ represents a differential stress and the temperature, T , dependence is determined through an activation energy, Q. R is the gas constant. Material constants, A, n, and Q can be determined experimentally (Kohlstedt et al., 1995). If n D 1 in Eq. 5 the viscosity can be calculated directly from the temperature and Eq. 5, but for the general case of n 6D 1, the constitutive law becomes nonlinear in stress. For the tensor version of the constitutive law, the stress, or strain rate, dependence can be written in terms of the second invariant of the tensor quantities. Using the square root of the second invariant of strain rate, the power-law viscosity, ¼pl , is defined as:   Q 1 0.1 n/ (6) ¼pl D I2 2n A n exp n RT where I20 is the second invariant of the deviatoric rate of deformation tensor defined as: 1 I20 D "i0 j "i0 j 2 Although the deviatoric rate-of-deformation is general, for the case of incompressible flow, the deviatoric and full rate-of-deformation are equivalent. The other deformational behavior proposed for crustal rocks is plastic, in particular plastic behavior governed by a Coulomb yield stress. Plastic deformation is fundamentally different from viscous deformation in that there is no deformation at stress levels less than the yield stress. For stress equal to the yield stress, strain rate is undefined, although many theories have been proposed for post-yield deformation. Since the problem here has been formulated as an Eulerian, velocity-based equilibrium problem, a natural formulation of the plastic deformation is provided by the Levy-Mises theory (Malvern, 1968). The Mises yield stress is defined as a limit on the second invariant of deviatoric stress, such that a scaler yield-stress, ¦ Y defines the material strength and the yield criterion is given by: q J20 D ¦ Y (7) where J20 is the second invariant of the deviatoric stress: 1 J20 D ¦i0j ¦i0j 2

423

The yield criterion of Eq. 7 is consistent with a viscous constitutive law (Eq. 2) if a stress-dependent viscosity is defined as: ¦Y ¼p D p 0 2 I2

(8)

The subscript on ¼ indicates that this is the effective viscosity for a plastic material. For geologic materials which exhibit frictional behavior at low temperature, a Coulomb yield criterion is more appropriate than the Mises criterion described above. This can be included by writing this criterion in terms of J20 . The Coulomb criterion can be written in terms of the three stress invariants, but for the case of 2-D, plane-strain, incompressible deformation, the third invariant is zero, so that the criterion can be written in terms of the first and second invariants of the deviatoric stress, J1 and J20 , respectively. q 1 J20 D c cos./ C J1 sin./ (9) 3 The cohesion, c, and the angle of friction,  are the two physical constants defining the Coulomb yield criterion. Eq. 9 is of the same form as Eq. 7 so the right side of Eq. 9 is the Coulomb equivalent of the Mises yield stress, ¦ Y , a scaler quantity, and Eq. 8 can be used to define an effective viscosity for a Coulomb plastic material. The equations of motion (Eqs. 3 and 4) with a choice of constant viscosity, a nonlinear power-law viscosity (Eq. 6), or a nonlinear, plastic viscosity (Eqs. 7–9) define a general statement of deformational behavior for a wide range of rheological behavior. Eqs. 3 and 4 are solved using standard finiteelement techniques for nonlinear viscous materials (Bathe, 1996). Details on the particular method can be found in Fullsack (1995).

3. Viscous wedge model With a linear strain rate response to stress and no spatial variation in material properties, a linear viscous rheology is the simplest material behavior to be considered in this study. The boundary and initial conditions for the formation of a doubly-vergent wedge are shown in Fig. 2. With this domain and boundary conditions, the equation of motion (Eq. 3)

424

S.D. Willett / Tectonophysics 305 (1999) 419–435

in terms of Ar and the other non-dimensional parameters. Note that the value of t Ł , the time parameter, corresponds to the total convergence in units of the layer thickness, so that t Ł of 5 corresponds to a total convergence of 5 times the layer thickness. Fig. 2. Boundary and initial conditions for traction reversal models. Domain is an infinite layer of thickness h with a free upper surface and specified velocity on lower surface. Horizontal velocity on base is constant at Vp for x < 0 and equal to 0 for x > 0. Vertical velocity is 0 on entire base.

can be non-dimensionalized using a characteristic length (h) and velocity (Vp ). With this scaling, nondimensional parameters include length: x Ł D x= h, time: t Ł D vp t= h, strain rate: "P Ł D "P h=vp , and a single parameter that includes the material properties and domain characteristics, the Argand number (England and McKenzie, 1982). The Argand number in this case is defined as (Buck and Sokoutis, 1994): Ar D

²gh 2 ¼vp

(10)

The Argand number in most cases reflects primarily the viscosity, as this is the parameter with the greatest potential range in geologic processes. Viscous models discussed below will be parameterized

3.1. Linear viscous models Using an Argand number of 10, a result similar to that of Buck and Sokoutis (1994) is obtained (Fig. 3). Fig. 3 shows the strain and strain rate fields in a format to be used throughout this paper. The lower figure shows the deformed domain along with the horizontal component of the instantaneous strain rate (rate of deformation) field. The upper two figures show the horizontal strain and strain rate along the upper surface defined to be the fractional change in length in the x (horizontal) direction and the change in the x component of velocity with respect to x, respectively. Fig. 3 shows deformation at a nondimensional time of 5, which corresponds to total shortening of 5 times the layer thickness. With this degree of convergence there is significant shortening and thickening with the formation of back to back wedges verging in opposite directions. Most of each wedge is experiencing horizontal shortening (light

Fig. 3. Traction reversal model for deformation of a crustal layer with a linear viscous rheology and an Argand number of 10. Results are shown for a non-dimensional time, t Ł of 5. Boundary and initial conditions are as shown in Fig. 2. Vertical ticks in strain rate and strain plots indicate point of basal velocity change. The lower figure shows the Lagrangian mesh with the horizontal component of the strain rate contoured with extension negative. Horizontal strain and horizontal strain rate of the upper surface are plotted above. Strain rate is scaled by the layer thickness, h, and the convergence velocity, Vp . Significant horizontal extension is present in the center of the model in both the strain rate and the total strain.

S.D. Willett / Tectonophysics 305 (1999) 419–435

425

Fig. 4. Traction reversal model with a linear viscous rheology and an Argand number of 100. Results are shown for a non-dimensional time, t Ł of 2. Horizontal extension is very large (fractional strain near 1.0) and evident in both strain and strain rate plots. See caption to Fig. 3 for other model details.

colors in lower Fig. 3), but at shallow depths and high elevations, the wedge is extending horizontally. The region of instantaneous extension is centered at the point of maximum elevation above the point of detachment of the deeper plate. The region of surface extensional strain (Fig. 3 center) is offset to the right from the instantaneous extension (Fig. 3 top) because of the left to right velocity of the deforming media which advects structure to the right. Extension in the near surface is quite large, with maximum strain reaching nearly 0.5, although the strain is depth dependent and extension is confined to shallow depths. Comparison of Fig. 3 with the numerical solution of Buck and Sokoutis (1994) (their fig. 3) indicates close agreement, consistent with the expected validity of the shallow channel flow assumption used in their paper, as well as providing a test of the finite element code used here. Both numerical solutions are consistent with the physical model presented in Buck and Sokoutis (1994). This result is parameter dependent with the magnitude of the surface extension being a strong function of the Argand number, Ar. The extension becomes more pronounced at larger Ar. With an Ar of 100 (Fig. 4) both contractional and extensional deformation are spread out over a much broader region. Horizontal extension is again present in the shallow crust at high elevation; in this case fractional strain values reach nearly 1.0 after a non-dimensional time of only 2.0 (Fig. 4). Correspondingly, a decrease in the Argand number decreases the extension. A model with an Argand number of 1 (Fig. 5) has a limited

region of instantaneous surface extension (Fig. 5, upper) and has no region of surface extensional strain. The lack of extensional strain reflects in part the large contraction experienced by material prior to extension. Since the deformed wedge is narrow with steep taper angles at low Ar, accreted material is highly contracted near the toe of either wedge. If this material is subsequently extended near the elevation high, the contraction must be fully reversed before net extension is obtained. In an orogenic belt characterized by low Ar, extensional structures would presumably be observed superimposed on earlier contractional structures, but this extension would be insufficient to result in any net crustal thinning. An even lower Argand number of 0.1 (Fig. 6) implies a very viscous crust and results in steep surface slopes. As with the model of Fig. 5, there is a slight amount of surface and near-surface extension, but no net extension of the upper surface of the model. 3.2. Power-law viscous wedges Rock deformational mechanisms show a constitutive behavior including linear viscous deformation by diffusional and pressure solution mechanisms. However, high-temperature deformation is frequently characterized by a dislocation creep mechanism which is described by a power dependence on deviatoric stress (Eq. 5). This nonlinear stress dependence is investigated in this section. To isolate the effect of the nonlinear stress dependence from the temperature dependence of dislocation creep, only isothermal models are presented. The Argand number is not

426

S.D. Willett / Tectonophysics 305 (1999) 419–435

Fig. 5. Traction reversal model with a linear viscous rheology and an Argand number of 1. Results are shown for a non-dimensional time, t Ł of 10. Horizontal extension is observed in the strain rate, but not strain. See caption to Fig. 3 for other model details.

constant for a nonlinear model, but by using Vc = h as an average, characteristic strain rate, an Argand number of 12 is obtained for the model of Fig. 7 which has a stress exponent, n, of 3. Use of a characteristic strain rate gives an Argand number which is correct only in average, such that this model exhibits

an average strength similar to the linear model with Ar D 10 (Fig. 3). In contrast to the linear viscous model, extension is greatly reduced, both in rate and in the net amount. Both are reduced by a factor of about 3 (note different scales on strain and strain rate plots). There is no net surface extension with

Fig. 6. Traction reversal model with a linear viscous rheology and an Argand number of 0.1. Results are shown for a non-dimensional time, t Ł of 5. Very slight extension of the upper surface is observed in strain rate. See caption to Fig. 3 for other model details.

S.D. Willett / Tectonophysics 305 (1999) 419–435

427

Fig. 7. Traction reversal model with a power-law flow rheology. Stress exponent is 3. Effective Argand number is 12. Results shown for a non-dimensional time, t Ł of 10. Extension of near surface is observed in strain rate and as a decrease in contraction, but with no net extension. See caption to Fig. 3 for other model details.

the power-law rheology, although the total surface shortening is clearly reduced in the region affected by extension. The decrease in extension is a characteristic of the stress dependence of the rheology, not the bulk properties. Even with a lower effective viscosity, corresponding to an effective Ar of 100. (Fig. 8),

extension is reduced relative to the comparable linear model (Fig. 4). With the lower effective viscosity (higher effective Ar) surface extension is at a rate sufficient to reverse the initial contraction and give net extension (Fig. 8, center panel), but, again, it is reduced relative to the model of Fig. 4, presumably reflecting the higher order of the stress dependence.

Fig. 8. Traction reversal model with a power-law flow rheology. Stress exponent is 3. Effective Argand number is 100. Results shown for a non-dimensional time, t Ł of 6. Extension is observed in both strain rate and strain with a maximum fractional strain of 0.25. See caption to Fig. 3 for other model details.

428

S.D. Willett / Tectonophysics 305 (1999) 419–435

4. Coulomb plastic wedge model Rocks at low-temperatures and low confining pressures deform in a brittle, plastic manner that is frequently described by a perfectly plastic constitutive law with a Coulomb yield stress (Eqs. 7 and 9). Deformation at the crustal scale has also be modeled with this material behavior. In particular, wedge mechanical models based on a Coulomb yield criterion have seen extensive development and application (Davis et al., 1983; Dahlen, 1984 and references therein). As discussed above, surface extension in critical Coulomb wedges has not been observed under conditions of continuous convergence, although doubly-vergent boundary conditions have not commonly been used (for exceptions, see discussions in Willett et al., 1993 and Wang and Davis, 1996). For this reason, a fairly extensive analysis of the mechanics of doubly-vergent Coulomb wedges is presented here. As in the previous models, the boundary and initial conditions for the Coulomb wedge models are as given in Fig. 2. However, in contrast to the viscous models, the material within the domain does not have constant properties. Critical Coulomb wedge theory is based on a model in which a wedge forms over a weaker basal detachment or decollement. In order to compare model results to critical wedge theory, similar conditions are imposed by including a layer along the base of the model with a lower coefficient of friction. There are thus two parameters that characterize the strength of the material, the internal angle of friction, , and the basal angle of friction, b . The plastic models do not scale in the same way as viscous models. There is no time or velocity dependence to the plastic deformation, so the problem is independent of convergence velocity and the appropriate time scaling is given by the total convergence, t Ł D Vc t= h. This also is evident by considering the Argand number for a plastic material. By substituting the definition of effective viscosity from Eq. 8 and the characteristic strain rate, Vc = h in the expression for Ar, (Eq. 9) we obtain: 2²gh Ar D (11) ¦Y Note that this expression is independent of the convergence velocity. It represents the ratio of two stresses, a characteristic pressure, ²gh and the yield

stress, ¦ Y . With a non-cohesive Coulomb yield stress and assuming that ²gh represents a characteristic normal stress, Eq. 11 reduces to Ar D 2= tan '. This implies that the coefficient of friction of the material completely characterizes the parameters in the models. The addition of a contrast in material property at the basal detachment introduces a second coefficient of friction and these two parameters characterize each of the models that follow. With a weak base and the basal boundary conditions with reversed traction, the development of the doubly-vergent Coulomb wedge follows a specific evolution as described in Willett et al. (1993) (Fig. 9). Initial deformation is focused at the detachment point. The plastic rheology results in strain localization into two conjugate shear-zones, the proand retro-shear zones that ‘step-up’ deformation to the surface (Fig. 9a). In an ideal plastic material these would be discrete faults, but in this case the localization is limited by the finite element size. With progressive deformation, the triangular plug between these shear zones is uplifted as more material is accreted to the left of the deformed region and the pro-shear zone propagates outward. Eventually, the initial plug moves to the right and a new retro-shear zone forms to the right (Fig. 9b). This process continues as more material is accreted into the deforming region; both shear-zones marking the edge of deformation propagate outward in discrete steps (Fig. 9a–d). Between these limiting shear zones deformation is also localized onto conjugate shear zones although some distributed deformation occurs in the blocks between them. The resulting structure of back-to-back wedges has a characteristic asymmetry with a relatively long, thin pro-wedge verging to the left and a shorter, steeper retro-wedge to the right. This asymmetry is reduced by a low-friction basal detachment as in this model; later models will illustrate a stronger asymmetry. The strain and strain rate at a mature stage of development of this model are shown in more detail in Fig. 10. The upper figures demonstrate that there is no horizontal extension of the upper surface in either the instantaneous or net strain. Strain and strain rate are highly variable across the upper surface in response to movement on the discrete shear zones (Fig. 10, lower), but are everywhere positive, indicating contraction.

S.D. Willett / Tectonophysics 305 (1999) 419–435

429

Fig. 9. Doubly-vergent Coulomb wedge model at non-dimensional times, t Ł of (a) 0, (b) 1, (c) 2, (d) 5. Internal angle of friction is 25º; basal angle of friction is 5º. Pro- and Retro- prefixes refer to features above the underthrusting (left) and overriding (right) plate respectively. Double-couple symbols show sense of shear in the deformation-limiting shear-zones. Boundary and initial conditions are shown in Fig. 2.

The geometry of either wedge in this system is a function of the internal angle of friction of the material and of the basal friction (Eq. 11). However, the taper angle is also a function of the kinematics of wedge formation as illustrated by the asymmetry of the two sides of the model. Coulomb wedges are stable during movement on a basal detachment for a range of taper angles. A critical taper is attained at the minimum taper angle only in the event of continuous accretion of new material into the toe of the wedge (Davis et al., 1983). A second critical wedge with a larger taper exists for any set of physical properties,  and b (Dahlen, 1984) and can be found in nature under appropriate kinematic conditions. The pro-wedge in the model of Fig. 11 is created by frontal accretion and therefore retains the minimum critical taper angle. The analytic theory predicts a minimum taper

angle of 2.1º and the pro-wedge in the model is in good agreement with that value (Fig. 11a). In contrast the retro-wedge is formed only in part by frontal accretion. Much of the material in the retro-wedge has passed over the detachment point from the prowedge. This material forms a steeper wedge just to the right of the highest point of elevation and reaches a taper close to the predicted maximum critical taper angle. Further to the right of this steep section, the wedge is again dominated by frontal accretion and the taper angle decreases to the minimum taper angle (Fig. 11). This close agreement with the analytical theory is maintained for other parameter values as demonstrated by the second example in Fig. 11 for a larger b of 15º. This pattern of a minimum critical taper pro-wedge and maximum and minimum tapers in the retro-wedge was also found in sandbox exper-

430

S.D. Willett / Tectonophysics 305 (1999) 419–435

Fig. 10. Strain and strain rate for the final configuration (t Ł D 5) of model in Fig. 9. Lower figure shows location and sense of shear on internal shear zones. Other figures as described in caption of Fig. 3. Strain rate is highly variable due to strain localization into discrete shear zones separating regions of low strain rate.

iments by Wang and Davis (1996). This consistency between sandbox experiments, analytical theory and numerical models provides important support for the use of this numerical approach in modelling plastic materials. The lack of horizontal extension in Coulomb wedge models is consistent across a range of parameter values. A low b as in the model of Fig. 10 encourages frontal accretion and hence overall contraction, but a higher b does not reverse the contraction. With a b of 15º (Fig. 12) the taper angles are all larger (Fig. 11), but the overall deformation pattern is similar to the case with b equal to 5º (Fig. 10). The largest taper angles that can be obtained with a material with an internal angle of friction of 25º is for b D  D 25º (Fig. 13). In this case, the theoretical minimum and maximum taper angles are equal and, consistent with theory, there is no detachment along the base of the model. There are numerical difficulties with this case as the retro-wedge is now at an angle approaching the angle of repose of the material and there appears to be small-scale, nearsurface failure of the model material which appears as extension in the strain rate plot (Fig. 13). This is even more apparent in a model with a lower angle of friction and no basal detachment, i.e. b D  D 15º (Fig. 14). There appears to be very large extension of

the surface of the model (Fig. 14), but this extension is confined to the very near-surface and is a form of down slope flow or slope failure. It is uncertain whether this is a numerical artifact, but it is unlikely that this type of deformation would be interpreted in terms of crustal scale extension.

5. Mixed rheology models The models of the preceding sections provide a useful documentation of the mechanical behavior of wedges with end-member rheologies, but geological materials in the crust inevitably show more complex behavior (Handy, 1989; Kohlstedt et al., 1995). For example, it is likely that deformation at low pressures and temperatures is by a Coulomb, plastic mechanism and at higher temperatures is by a thermally-activated viscous mechanism. As a viscous rheology predicts significant extension during convergence and a plastic rheology does not, it is important to determine which mechanism controls the crustal deformation pattern. A model with a mixed rheology is shown in Fig. 15. In this model, a viscous-plastic rheology has been used. Material deforms according to a nonlinear viscous constitutive law of the form of Eq. 5 un-

S.D. Willett / Tectonophysics 305 (1999) 419–435

431

Fig. 11. Elevation sections of Coulomb models with internal angle of friction of 25º and basal angle of friction of (a) 5º and (b) 15º. Models are shown in Figs. 10 and 12, respectively. Elevation is shown for non-dimensional times, t Ł of 1 through 5. Analytical solutions from Dahlen (1984) are shown as dashed lines. Both models show a minimum taper pro-wedge to left and a variable taper retro-wedge to the right with maximum and minimum taper sections. Agreement with analytical solution is good at non-dimensional times greater than about 3.

Fig. 12. Coulomb wedge model with internal angle of friction of 25º and a basal angle of friction of 15º. Wedge slopes are steeper than the model of Figs. 9 and 10, consistent with predictions of critical wedge theory. Retro-wedge shows both maximum and minimum tapers. No horizontal extension is observed in model.

432

S.D. Willett / Tectonophysics 305 (1999) 419–435

Fig. 13. Coulomb wedge model with basal angle of friction equal to the internal angle of friction which is 25º. Deformation is primarily localized on the limiting shear-zones. A small zone of surface extension is predicted on the steep retro-wedge slope due to near-surface slope failure.

less the Coulomb yield stress is exceeded, at which point the plastic, effective viscosity (Eq. 8) is used. Power-law parameters appropriate for a quartz rhe-

ology were used (Kohlstedt et al., 1995) including a stress exponent of 3 and an activation energy of 135 kJ=mol. A constant geothermal gradient of 20ºC=km

Fig. 14. Coulomb wedge model with basal angle of friction equal to the internal angle of friction which is 15º. Deformation is primarily localized on the limiting shear-zones. Surface extension on the retro-wedge is pronounced, but is confined to shallow depths and represents near-surface slope failure. Note absence of extension in strain rate plot indicating that the observed extension formed during a transient event.

S.D. Willett / Tectonophysics 305 (1999) 419–435

433

Fig. 15. Viscous-plastic rheology model. Material deforms according to a temperature-dependent power-law flow rule at stresses less than the Coulomb yield stress. Changes in crustal thickness are isostatically compensated. Note that plastic deformation in the shallow crust precludes any surface extension that otherwise would be present with a viscous rheology.

was assumed through the 40 km thick layer. An internal angle of friction of 20º was used for the Coulomb failure criterion, consistent with Byerlee’s Law with fluid pressure above hydrostatic (Kohlstedt et al., 1995). In addition the lower boundary condition was changed to allow for isostatic compensation of crustal thickening. Models presented earlier did not include isostatic compensation in order to keep models simple and to allow direct comparison to the analytical solutions of Buck and Sokoutis (1994) and Dahlen (1984). Crustal thickening in these last models is compensated regionally using the analytic solution for continuous plate flexure (Turcotte and Schubert, 1982) with a small flexural rigidity was used so that compensation is nearly local. The resulting orogen model (Fig. 15) shows characteristics of both deformational mechanisms. Strain is localized into bounding shearzones, but between these, deformation is more distributed than in the completely plastic case. Note that the ‘wave’ structure in the surface strain plot is due to structure developed during initiation of deformation. No horizontal extension is observed in the model implying that the plastic deformation at shallow depths precludes surface extension. This observation holds even if the crust is weaker on average. Increasing the effective Argand number of the crust by an order of magnitude gives a much thinner, broader model orogen (Fig. 16), but does not lead to extension. Again the Coulomb, plastic defor-

mation of the uppermost crust prevents any extension potentially driven by the viscous flow at depth.

6. Discussion and conclusions The suite of models presented here demonstrates that, although extension can occur in the context of orogenic wedge mechanics with homogeneous, isotropic materials, it can only occur under restricted rheological conditions. For the conditions studied here, traction reversal and formation of a doublyvergent orogenic wedge, extension occurs at shallow depths beneath the highest surface elevation most readily for weak, viscous materials. Specifically, extension occurs in linear-viscous wedges characterized by an Argand number greater than 1.0. Powerlaw, viscous materials also show extension, but at reduced rates and with lower net extensional strain. In contrast, Coulomb plastic materials exhibit no extensional deformation during wedge formation. This observation appears to hold for any specific yield parameters and reflects a fundamental rheological behavior difference in response to these applied boundary conditions, rather than reflecting an average strength of the material. Although only a single set of boundary conditions have been studied here, the conclusions likely apply to other mechanical models as well. The existence

434

S.D. Willett / Tectonophysics 305 (1999) 419–435

Fig. 16. Viscous-plastic rheology model. Material deforms according to a temperature-dependent power-law flow rule at stresses less than the Coulomb yield stress. Flow-law parameters are such that the average viscosity is an order of magnitude smaller than the model in Fig. 15. Plastic deformation in the shallow crust precludes any surface extension.

of a significantly thick Coulomb plastic layer above a viscous lower crust is likely to inhibit extension in the models of Royden (1996). As demonstrated by analytical critical wedge theory, Coulomb wedges have a broad region of stability that requires dramatic changes in strength or geometry in order to extend. Numerical modeling using the boundary conditions of Royden (1996), although not presented here, supports this view (Schneider and Willett, 1996). These models present implications for extension in convergent orogens. The significance of the implications depends in part on how well the bulk rheology of orogenic belts at a crustal scale can be determined. If a significant thickness of the upper crust deforms by a Coulomb plastic mechanism, it is difficult to envision steady convergence and the formation of an orogenic wedge leading to significant extensional strain. On the other hand, if viscous mechanisms dominate deformation even at shallow depths, the orogenic wedge could exhibit extension during convergence. Viscous rheologies need not be limited to high temperature mechanisms such as diffusion or dislocation creep, although the high erosion rates and high geothermal gradients common to orogenic belts do encourage high temperature mechanisms at shallow depths. It seems more likely that viscous deformation at very shallow depths (within a few kilometers of the Earth’s surface) would be by a pressure solution mechanism. In examples of limited

extension in orogenic belts such as the modern Taiwan collision (Crespi et al., 1996; Yu et al., 1995), the wedge traction-reversal model is attractive and the importance of pressure solution or other viscous deformational mechanisms should be explored.

Acknowledgements This research was supported by the National Science Foundation through grant EAR 94-17766. The manuscript was improved by helpful reviews by Roger Buck and Robert Holdsworth.

References Angelier, J., Liberis, N., Le Pichon, X., Barrier, E., Huchon, P., 1982. The tectonic development of the Hellenic Arc and the sea of Crete: A synthesis. Tectonophysics 86, 159–196. Bathe, K.-J., 1996. Finite Element Procedures. Prentice Hall, New Jersey, 1037 pp. Beaumont, C., Fullsack, P., Hamilton, J., 1995. Styles of crustal deformation caused by subduction of the underlying mantle. Tectonophysics 232, 119–132. Beaumont, C., Quinlan, G., 1994. A geodynamic framework for interpreting crustal-scale seismic-reflectivity patterns in compressional orogens. Geophys. J. Int. 116, 754–783. Brown, R.L., Carr, S.D., 1990. Lithospheric thickening and orogenic collapse within the Canadian cordillera. Pacific Rim 90 Congress, Proc., V. Ii: 1–10.

S.D. Willett / Tectonophysics 305 (1999) 419–435 Buck, W.R., Sokoutis, D., 1994. Analogue model of gravitational collapse and surface extension during continental convergence. Nature 369, 737–740. Burchfiel, B.C., Royden, B.C., 1985. North–south extension within the convergent Himalayan region. Geology 13, 679– 682. Chapple, W.M., 1978. Mechanics of thin-skinned fold-and-thrust belts. Geol. Soc. Am. Bull. 89, 1189–1198. Chemenda, A., Mattauer, M., Bokun, A., 1996. Continental subduction and a mechanism for exhumation of high-pressure metamorphic rocks; new modelling and field data from Oman. Earth Planet. Sci. Lett. 143, 173–182. Chemenda, A., Mattauer, M., Malavieille, J., Bokun, A., 1995. A mechanism for syn-collisional rock exhumation and associated normal faulting; results from physical modelling. Earth Planet. Sci. Lett. 132, 225–232. Crespi, J.M., Chan, Y.-C., Swaim, M., 1996. Synorogenic extension and exhumation of the Taiwan hinterland. Geology 24, 247–250. Dahlen, F.A., 1984. Noncohesive critical Coulomb wedges: An exact solution. J. Geophys. Res. 89, 10125–10133. Dahlen, F.A., 1990. Critical taper model of fold-and-thrust belts and accretionary wedges. Annu. Rev. Earth Planet. Sci. 18, 55–99. Davis, D., Suppe, J., Dahlen, F.A., 1983. Mechanics of fold-andthrust belts and accretionary wedges. J. Geophys. Res. 88, 1153–1172. Dewey, J.F., 1988. Extensional collapse of orogens. Tectonics 7, 1123–1139. Dewey, J.F., Sengor, A., 1979. Aegean and surrounding regions; complex multiplate and continuum tectonics in a convergent zone. Geol. Soc. Am. Bull. 90, I84–192. Doglioni, C., 1994. Foredeeps versus subduction zones. Geology 22, 271–274. England, P., McKenzie, D.P., 1982. A thin viscous sheet model for continental deformation. Geophys. J. R. Astron. Soc. 70, 295–321. Handy, M.R., 1989. Deformation regimes and the rheological evolution of fault zones in the lithosphere: the effects of pressure, temperature, grain size and time. Tectonophysics 163, 119–152. Fullsack, P., 1995. An arbitrary Lagrangian-Eulerian formulation for creeping flows and its applications in tectonic models. Geophys. J. R. Astron. Soc. 120, 1–23. Hodges, K.V., Parrish, R.R., Housh, T.B., Lux, D.R., Burchfiel, B.C., Royden, L.H., Chen, Z., 1992. Simultaneous Miocene extension and shortening in the Himalayan orogen. Science 258, 1466–1470. Houseman, G.A., McKenzie, D.P., Molnar, P., 1981. Convective instability of a thickened boundary layer and its relevance for the thermal evolution of continental convergent belts. J. Geophys. Res. 86, 6115–6132. Kohlstedt, D.L., Evans, B., Mackwell, S.J., 1995. Strength of the lithosphere: Constraints imposed by laboratory experiments. J. Geophys. Res. 100, 17587–17602. Malavieille, J., 1984. ModKelisation expe´rimentale des chevauchements imbrique´s: Application aux chaıˆnes des montagnes. Soc. Geol. Fr. Bull. 26, 129–138.

435

Malavielle, J., 1993. Late orogenic extension in mountain belts: insights from the basin and range and the late Paleozoic Variscan belt. Tectonics 12, 1115–1130. Malvern, 1968. An Introduction to the Mechanics of a Continuous Medium. Prentice-Hall, New Jersey. Mercier, J.-L., Armijo, R., Tapponnier, P., Carey-Gailhardis, E., Tonglin, H., 1987. Change from Tertiary compression to Quaternary extension in southern Tibet during the India–Asia collision. Tectonics 6, 275–304. Molnar, P., England, P., Martinod, J., 1993. Mantle dynamics, uplift of the Tibetan Plateau, and the Indian Monsoon. Rev. Geophys. 31, 357–396. Platt, J.P., 1986. Dynamics of orogenic wedges and the uplift of high-pressure metamorphic rocks. Geol. Soc. Am. Bull. 79, 1037–1053. Platt, J.P., Vissers, R.L.M., 1989. Extensional collapse of thickened continental lithosphere; a working hypothesis for the Alboran Sea and Gibraltar Arc. Geology 17, 540–543. Ratschbacher, L., Frisch, W., Neubauer, F., Schmid, S.M., Neugebauer, J., 1989. Extension in compressional orogenic belts: The eastern Alps. Geology 17, 404–407. Royden, L., 1993. Evolution of retreating subduction boundaries formed during continental collision. Tectonics 12, 629–638. Royden, L., 1996. Coupling and decoupling of crust and mantle in convergent orogens: Implications for strain partitioning in the crust. J. Geophys. Res. 101, 17679–17705. Royden, L., Burchfiel, B.C., King, R.W., Wang, E., Chen, Z., Shen, F., Liu, Y., 1997. Surface deformation and lower crustal flow in eastern Tibet. Science 276, 788–790. Schneider, C., Willett, S.D., 1996. Rheological controls on extension in convergent orogens. Trans. Am. Geophys. Union 77, S268. Turcotte, D.L., Schubert, G., 1982. Geodynamics: Applications of Continuum Physics to Geological Problems. Wiley, New York. Villote, J.P., Daignieres, M., Madariaga, R., 1982. Numerical modeling of intraplate deformation: simple mechanical models of continental collision. J. Geophys. Res. 87, 10709–10728. Vissers, R.L.M., Platt, J.P., van der Wal, D., 1995. Late orogenic extension of the Betic Cordillera and the Alboran Domain: A lithospheric view. Tectonics 14, 786–803. Wang, W.-H., Davis, D.M., 1996. Sandbox model simulation of forearc evolution and non-critical wedges. J. Geophys. Res. 11, 11329–11339. Waschbusch, P., Beaumont, C., 1996. Effect of a retreating subduction zone on deformation in simple regions of plaste convergence. J. Geophys. Res. 101, 28133–28148. Willett, S.D., 1992. Dynamic and kinematic growth and change of a Coulomb wedge. In: Thrust Tectonics. Chapman and Hall, London, pp. 19–31. Willett, S.D., Beaumont, C., 1994. Subduction of Asian lithospheric mantle beneath Tibet inferred from models of continental collision. Nature 369, 642–645. Willett, S.D., Beaumont, C., Fullsack, P., 1993. Mechanical model for the tectonics of doubly vergent compressional orogens. Geology 21, 371–374. Yu, S.-B., Chen, H.-Y., Kuo, L-C., 1995. Velocity field of GPS stations in the Taiwan area. ACT Symp., pp. 317–327.