Stresses and velocities in orogenic wedges with power-law rheology and linearly varying longitudinal strain rate

PII:

Journal of Structural Geology, Vol. 20, No. 12, pp. 1611 to 1623, 1998 # 1998 Elsevier Science Ltd. All rights reserved 0191-8141/98/$ - see front matter S0191-8141(98)00036-4

Stresses and velocities in orogenic wedges with power-law rheology and linearly varying longitudinal strain rate JIYANG LIU and GIORGIO RANALLI Department of Earth Sciences and Ottawa±Carleton Geoscience Centre, Carleton University, Ottawa K1S 5B6, Canada E-mail: [email protected] (Received 7 November 1996; accepted in revised form 30 March 1998) AbstractÐA two-dimensional analytical solution for stress, strain rate, and velocity is obtained for parallelsided and wedge-shaped blocks with generalized viscous rheology (linearly viscous and power-law) deforming in plane strain. The main assumptions used in the derivation of the solution are that the material is incompressible, the longitudinal gradient in shear stress is much less than the vertical gradient of vertical normal stress, and the longitudinal strain rate varies linearly in the horizontal direction. Velocity boundary conditions are speci®ed at the top of the block, and shear stress boundary conditions at the base of the block. In the onedimensional case (where stress and strain rate do not vary in the longitudinal direction), the solution reduces to a well-known solution for the deformation of parallel-sided ice sheets [Nye, J. F. (1957) The distribution of stress and velocity in glaciers and ice sheets. Proceedings of the Royal Society of London A-239, 113±133]. The stress equilibrium for tapered wedges [Platt, J. P. (1986) Dynamics of orogenic wedges and the uplift of highpressure metamorphic rocks. Geological Society of America Bulletin 97, 1037±1053] is a special case of the present stress solution. Implementation of the solution requires the subdivision of the wedge into vertical segments, and yields the tectonic normal and shear stresses that must be applied to the rear of a block with speci®ed rheology in order to maintain a given longitudinal strain rate. The solution makes it possible to model deformation patterns analytically with longitudinally varying strain rate (including coeval compression and extension) and with vertical components of velocity re¯ecting the e€ects of underplating. # 1998 Elsevier Science Ltd. All rights reserved

INTRODUCTION

(1957) for incompressible plane strain ¯ow. Because of the similarities between the movement of glaciers and Most models for orogenic wedges, notably the critical ice sheets and that of nappes and thrust sheets (Elliott, taper model, have focused on the dynamic equilibrium 1976), this solution has been applied to the analysis of of foreland fold-and-thrust belts and accretionary shortening and elongating thrust sheets (Wojtal, prisms, i.e. supracrustal thin-skinned orogenic wedges 1992a,b). Since it does not specify a priori the stress which are usually less than 5±10 km thick at the rear exponent n for the material, it is suciently general to (Elliott, 1976; Chapple, 1978; Davis et al., 1983; accommodate a wide range of deformation modes Dahlen, 1984, 1990; Dahlen et al., 1984; Dahlen and which might be applicable to thick orogenic wedges: Suppe, 1988; Zhao et al., 1986; Yin, 1989; Liu and linear Newtonian ¯ow (n = 1), power-law creep Ranalli, 1992). These thin-skinned wedge models, how- (n = 20 4), and perfect plasticity (n 4 1). However, ever, do not account quantitatively for some features Nye's solution is applicable only to parallel-sided often observed in thick-skinned orogenic wedges invol- sheets with uniform shortening or extending ¯ows, ving basement rocks (Pavlis and Bruhn, 1983; Platt, since it is based on the assumptions that the stress ten1986, 1987, 1993), namely: (a) coeval development of sor depends on depth only and that the vertical velcompressional structures at the thin-skinned foreland ocity component does not vary with longitudinal `front' and extensional structures at the thick-skinned position. The constancy of longitudinal strain rate `rear' towards the hinterland; (b) exhumation of high makes this solution unsuitable to the modelling of pressure±low temperature metamorphic rocks near the thick-skinned orogenic wedges where the deformation rear of the wedge; and (c) nonlinear rheology of ma- pattern varies not only in the vertical but also in the terials (which is Coulomb frictional only in the upper longitudinal direction. parts, and most likely nonlinearly viscous in the bulk In this paper we present a two-dimensional plane of the wedge; Platt, 1986, 1987, 1993; Liu and Ranalli, strain analytical solution for stresses and velocities in 1994; Ranalli, 1995; Liu, 1996). blocks and wedges with nonlinear rheology. This solThus, the need exists for models incorporating the ution removes some of the conditions limiting Nye's nonlinear bulk rheology of thick-skinned wedges. (1957) solution, and is therefore applicable to cases Analytical solutions for nonlinear (power-law) viscous where the strain rate varies in the longitudinal direcrheology exist in glaciology (Nye, 1951, 1957, 1969; tion. It reduces to Nye's solution in the one-dimenBudd, 1970; Shoemaker and Morland, 1984). One of sional case. The focus of the argument is on the these is the classic ice sheet solution obtained by Nye procedure for obtaining a closed-form solution, and 1611

1612

J. LIU and G. RANALLI e_ xx ˆ ls0xx

not on computational details. However, the results of four model tests are given, to show the potential applicability of the solution to thick-skinned orogenic wedges. GOVERNING EQUATIONS The geometry of the model, coordinate system, and sign convention for stresses are given in Fig. 1. The equilibrium condition for two-dimensional incompressible power-law ¯ow is (cf. e.g. Ranalli, 1995) @sxx @sxy ‡ ‡ rgx ˆ 0 @x @y @sxy @syy ‡ ‡ rgy ˆ 0 @x @y

e_ yy ˆ ls0yy e_ xy ˆ ls0xy

…2†

where s'ij is the stress deviator, Çeij is the strain rate tensor, and l is a stress-dependent parameter having dimensions of reciprocal viscosity. The strain rate tensor is related to velocity components u, v in the x, y direction as e_ xx ˆ

@u @x

@v @u ˆÿ @y @x   1 @u @v ‡ ˆ 2 @y @x

e_ yy ˆ …1†

where sij is the stress tensor, r the density, and gi the acceleration of gravity. The rheological equations are

e_ xy

and the stress deviator is given by

Fig. 1. (a) Simpli®ed two-dimensional orogenic wedge model: a is the surface slope, b the basal slope, u, v the velocity components, and g the acceleration of gravity; positive stress components are shown. (b) Parallel-sided block and (c) wedge-shaped block. tb is basal shear stress, h is local wedge thickness.

…3†

Stresses and velocities in orogenic wedges
1ÿ sxx ÿ syy 2
1ÿ ˆ syy ÿ sxx ˆ ÿs0xx 2

s0xx ˆ s0yy

s0xy ˆ sxy

…4†

The explicit form of the parameter l in equation (2) for power-law creep is (cf. e.g. Ranalli, 1995)    Q ‡ PV 1 0 0 …n ÿ1†=2 s s …5† l ˆ A0 exp ÿ RT 2 ij ij where T is absolute temperature, P pressure, R the universal gas constant, and Q, V, A0 and n are material parameters [Q and V are the creep activation energy and activation volume, respectively; the term PV is negligible when P<
A TWO-DIMENSIONAL LINEAR LONGITUDINAL STRAIN RATE SOLUTION

1613

usually small (aR 68, cf. Boyer and Elliot, 1982; Suppe, 1985), this has no major consequences. The ®rst assumption [equation (6)] greatly simpli®es the problem while still allowing the simulation of synconvergence extension, which is a ®rst-order feature of many orogenic wedges (Pavlis and Bruhn, 1983; Platt, 1986, 1987, 1993). The introduction of this assumption implies that the problem becomes that of ®nding force boundary conditions at the rear of the wedge resulting in a linearly varying longitudinal strain rate and realistic fault patterns. Consequently, the validity of the assumption is tested a posteriori from the success of the solution to model realistic tectonic forces and deformation. The second assumption [equation (7)] is generally held to be correct to the ®rst order and has been adopted in numerous studies (Elliott, 1976; Boyer and Elliott, 1982; Davis et al., 1983; Dahlen, 1984; Suppe, 1985; Platt, 1986; Liu and Ranalli, 1992, 1994; Yin, 1988, 1989, 1993). The third assumption [equation (8)] should be a good approximation whenever the material is well lithi®ed (cf. Dahlen and Barr, 1989 for details); however, it is not valid when there is large volume reduction due to loss of ¯uids. The fourth assumption [equation (9)] is simply a ®rst approximation to the usual density variations in the upper crust of the Earth.

Basic assumptions

Derivation of solution

Four assumptions are adopted in the derivation of the present solution:

A general stress solution to equation (1) under the second assumption [equation (7)] is

1. The longitudinal strain rate Çexx varies linearly in the longitudinal direction at any point within the wedge

sxx ˆ ÿ rgy y ‡ F…y† ÿ xG0 …y†

@2 e_ xx ˆ0 @x2

rgy y syy ˆ ÿ

…6†

2. The gradient of shear stress in the x-direction is negligible compared with the gradient of vertical normal stress in the y-direction @sxy << @syy …7† @x @y 3. The wedge material is incompressible and in plane strain @u @v ‡ ˆ0 …8† @x @y

rgx y ‡ G…y† sxy ˆ ÿ

where r is the depth-averaged density and F(y), G(y) are two arbitrary functions of depth. Physically, F(y) and G(y) are the non-gravitational (i.e. tectonic) normal and shear stress components applied on the rear of the block (Fig. 2). Using the stress solution [equation (10)] in the expressions for strain rate [equation (2)] and taking into account equation (4), the strain rate components are

e_ yy e_ xy

…9†

It should be noted that the terms `horizontal' and `vertical' are used in this paper in the sense of `x-direction' and `y-direction'. Since the surface slope angle a is

 l F…y† ÿ xG0 …y† 2  l ˆ ls0yy ˆ ÿ F…y† ÿ xG0 …y† 2   0  gx y ˆ lsxy ˆ l G…y† ÿ r

e_ xx ˆ ls0xx ˆ

4. The density varies with depth only @r ˆ 0: @x

…10†

…11a†

where the parameter l as de®ned in equation (5), can be expressed in terms of the stress solution as  l ˆ B0

2  2 1  gx y F…y† ÿ xG0 …y† ‡ G…y† ÿ r 4

…n ÿ1†=2

…11b†

1614

J. LIU and G. RANALLI x2 b1 a1 …a1 ‡b1 y†ÿx…a2 ‡b2 y†ÿ y3 ÿ y2 ÿb3 y 2 6 2  ) Z z Z y( a2 ‡b2 z a1 ‡b1 t F…t† dt dz rgx ‡ ÿ a3 ‡4 F…z† 0 a2 ‡b2 t 0

uˆÿ

vˆ ‡

Fig. 2. State of stress along the wedge rear boundary (x = 0). The stress axis refers to magnitude. According to the adopted sign convention, F(y) < 0 when compressive; G(y) < 0 when shearing downwards.

  where B0 ˆ A0 exp ÿ Q‡PV is a function of temperaRT ture, pressure and material parameters [see equation (5)]. The general expressions for the velocity components are obtained by integrating the ®rst two of equation (3) with respect to x and y, respectively, Z   1 x l…x; y† F…y† ÿ xG0 …y† dx ‡ H…y† u…x; y† ˆ 2 0 Z   1 y l…x; t† F…t† ÿ xG0 …t† dt ‡ I…x† …12† v…x; y† ˆ ÿ 2 0 where I(x) and H(y) are two integration functions. Physically, H(y) is the vertical distribution of the longitudinal velocity component at the rear of the block, since H(y) = u(x = 0, y); and I(x) is the longitudinal distribution of vertical velocity component at the surface, since I(x) = v(x, y = 0). Combination of the four assumptions [equations (6)± (9)] and appropriate boundary conditions allow all the integration functions appearing in equations (10)±(12) to be expressed as functions of F(y) and seven arbitrary constants (a1, a2, a3, b1, b2, b3, b4). The details of the derivation are given in the Appendix. The ®nal expressions for stress, strain rate and velocity are   a1 ‡ b1 y x F…y† rgy y ‡ 1 ‡ sxx ˆ ÿ a2 ‡ b2 y syy ˆ ÿ rgy y  Z y a1 ‡ b1 t F…t† dt rgx ‡ sxy ˆ ÿ a2 ‡ b2 t 0

…13a†

e_ xx ˆ ÿ…a1 ‡ b1 y†x ÿ …a2 ‡ b2 y† e_ yy ˆ ÿ_exx e_ xy

2…a2 ‡ b2 y† ˆ F…y†

Z y 0

 a1 ‡ b1 t F…t† dt rgx ‡ a2 ‡ b2 t

…13b†

  b1 3 b2 2 b1 2 x ‡ x ‡x y ‡ a1 y ‡ b3 6 2 2 b2 2 y ‡ a2 y ‡ b4 2

…13c†

The parameter l in equation (2) can be written in terms of F(y) and arbitrary constants as l = ÿ 2(a2+b2y)/F(y), which, when combined with equations (5) and (11b), links the function F(y) with the rheology of the block through the following expression ( #2  2…a2 ‡ b2 y† F…y† a1 ‡ b1 y x ˆ 1‡ ÿ F…y†B0 2 a2 ‡ b2 y "Z   2 )…n ÿ1†=2 y a1 ‡ b1 t ‡ : F…t† dt rgx ‡ a2 ‡ b2 t 0

…13d†

Equations (13a±d) constitute a complete general solution for a two-dimensional block with power-law rheology and linear longitudinal strain rate. The general solution is equally applicable to rectangular blocks and to triangular wedges (refer to Fig. 1). Discussion of the solution Several points about equation (13) should be noted: 1. There are two kinds of unknowns which are involved in equation (13), and need to be determined from boundary conditions: the function F(y), and seven constants of integration. The function F(y) is related to some of the integration constants [cf. equation (13d)]. According to their kinematic implications, these seven constants can be classi®ed into two groups: those related to longitudinal strain rate (a1, a2, b1 and b2), and those related to rigidbody motion corresponding to velocity and rotation at the origin (a3, b3 and b4). 2. The state of stress at the rear of the block (x = 0), represented by the functions F(y) and G(y), is related to rheology and strain rate [cf. equation (13d) and (A6) in the Appendix]. The solution therefore requires F(y) to be computed for a given rheology and prescribed strain rate. Physically, therefore, F(y) is the tectonic normal stress required at the rear of a block with given material rheology and temperature distribution to maintain a given longitudinal strain rate. Another possible way of applying

Stresses and velocities in orogenic wedges equation (13d) is to compute the strain rate for a given rheology and independently speci®ed F(y). However, this approach leads to severe diculties, since it is not feasible to ensure the compatibility of the speci®ed F(y) with the given rheology and temperature distribution. 3. The assumption of linear longitudinal strain rate [equation (6)] requires the parameter l to be a function of y only [see equation (A5) in the Appendix]. This results in the left side of equation (13d) being a function of y only while the right side is a function of both x and y. This apparent paradox is the direct consequence of the linearity assumption for the longitudinal strain rate. To demonstrate the impact of this assumption on the relation between strain rate and deviatoric stress [i.e. equation (2)], we have derived other closed-form solutions based on di€erent assumptions for the longitudinal strain rate. For example, assuming that the longitudinal strain rate at any depth y is a cubic polynomial function of position x (i.e. @4 e_ xx =@x4 ˆ 0), we obtain another set of solutions for stress, strain rate and velocity. Both sides of the corresponding equation (13d) for this new set of solutions are functions of x and y. Although it is much more complicated, this new set of solutions can be reduced to the simpler solution reported in this paper if we degrade the assumption on the longitudinal strain rate from a cubic polynomial function of x into a linear function of x. Obviously, the complexity of the solution increases with the complexity of the assumed x-dependence of longitudinal strain rate. On the other hand, the signi®cant bene®t of the linearity assumption is the relative simplicity of the derived analytical solution. Numerically, this diculty is solved by dividing the wedge into vertical segments of width Dx and applying equation (13) to each segment, within which the longitudinal variation of deviatoric normal stress is negligible. Physically, the interdependence of tectonic stress and longitudinal strain rate is a consequence of the rheological properties of the wedge. Mathematically, the determination of F(y) from equation (13d) is a two-point boundary value problem, which is solved here by the so-called shooting method (for detailed discussion see Press et al., 1992). For a typical wedge segment with a longitudinal width less than 10 km, the longitudinal variation of the computed function F(y) is less than 5% (Liu, 1996). A smooth variation of the solution for stress, strain rate and velocity from one segment to the next supports its validity. Continuity for adjacent Dx is achieved by linear interpolation. Two sets of boundary conditions are speci®ed for each segment: three points on the top boundary with assigned velocity and one point on the basement with assigned basal shear stress. The state of stress along

1615

the basement, together with the smoothness of the solution in the whole wedge, is another test of the validity of the present approach. The required basal shear stress tb should be smooth, either constant or increasing from the front to the rear, and of reasonable magnitude (0±100 MPa).

COMPARISON BETWEEN TWO-DIMENSIONAL AND ONE-DIMENSIONAL SOLUTIONS The linear longitudinal strain rate solution given by equation (13) is a two-dimensional extension of the classic one-dimensional solution for glacier ¯ow (Nye, 1957). Table 1 compares the two solutions in detail, and shows that the present solution reduces to Nye's for a1=b1=b2=b3=0, which is equivalent to the assumption @v/@x = 0 used in Nye's analysis. In the one-dimensional solution, both the second stress invariant and the shear stress are assumed to be independent of longitudinal position. In the two-dimensional solution, the less restrictive assumption that the longitudinal gradient in shear stress is negligible compared with the vertical gradient in vertical normal stress results in a horizontal normal deviatoric stress which is a linear function of longitudinal position. This allows a more complex velocity ®eld to be accommodated.

COMPARISON WITH EQUILIBRIUM CONDITION FOR TAPERED WEDGES In this section, we compare the present solution with Platt's dynamic stability criterion for orogenic wedges (Platt, 1986), and show that this criterion is a particular case of the present stress solution. Assuming an unspeci®ed viscous rheology, Platt (1986) derived a stability criterion for two-dimensional plane strain thick orogenic wedges, which relates the basal shear stress (tb) to the topographic slope (a), the depth-averaged longitudinal deviatoric stress (txx) and its longitudinal gradient acting in a tapered wedge (cf. Fig. 3) tb ˆ rgha ÿ 2txx y ÿ 2

@txx h @x

…14†

where y(=a + b) is the taper of the wedge, h the local thickness of the wedge, r density, and g gravity. It should be pointed out that an additional implicit assumption is needed in Platt's model in order to derive equation (14). This is that the longitudinal deviatoric stress txx be equal to zero on the upper surface of the wedge. This assumption is used when integration with respect to y is performed on the stress equilibrium equation in the x-direction [see equations (11) and (12) in Platt, 1986]. This additional condition means that the near-surface part of the

1616

J. LIU and G. RANALLI Table 1. Comparison between one-dimensional solution (Nye, 1957) and the present two-dimensional solution.

*Minor modi®cations to Nye's original expressions have been made to facilitate comparison.

wedge is always under a nearly lithostatic state of stress. To recast the present stress solution [equation (13a)] in a format that is directly comparable with Platt's model, the longitudinal deviatoric stress (s'xx) and the depth-averaged longitudinal deviatoric stress (txx) are computed from equation (13a)

F…y†ÿxG0 …y† 2 " Z # Z @txx @ 1 h 0 1 h 0 G…h† ˆ sxx dy ˆ ÿ G …y†dy ˆ ÿ @x @x h 0 2h 0 2h s0xx ˆ

…15†

where h is the thickness of the wedge at the point

Stresses and velocities in orogenic wedges

di€erence in coordinate systems (compare Figs 1c & 3). The second di€erence is caused by the aforementioned assumption about the stress state along the wedge upper surface. Therefore, Platt's (1986) stability criterion is compatible withÐindeed, it is a particular case ofÐthe stress solution derived in this paper, subject to the assumption that the longitudinal deviatoric stress does not vary signi®cantly along any depth pro®le. This condition most likely holds in thin-skinned wedges but may not be satis®ed in thick-skinned wedges, where deformation patterns may vary signi®cantly in the vertical direction.

Fig. 3. Platt's accretionary wedge model and its coordinates (Platt, 1986): a and b are the top and basal slopes, respectively; y(=a + b) is the taper of the wedge, h thickness, and tb basal shear stress.

under consideration. The basal shear stress tb is related to the state of stress within the wedge as (cf. e.g. Ranalli, 1995)
ÿ
1ÿ syy ÿ sxx yˆh sin 2y ÿ sxy yˆh cos 2y 2   1 rgha ÿ y F…y† ÿ xG0 …y† yˆh ÿG…h†

EXAMPLES Four examples of the two-dimensional solution are discussed in this section. The ®rst two are rectangular blocks with Newtonian and power-law rheologies; the remaining two are wedges with Newtonian and powerlaw rheologies. The parameter values used in the computations are listed in Table 2. All cases are parameterized in the following way:

tb ˆ

…16†

Combining equations (15) and (16) leads to ÿ
@txx gha ÿ 2y s0xx yˆh ‡2h tb ˆ r @x

1617

…17†

(a) The velocity boundary conditions along the top surface are of the general form

where small-angle approximations for a and y are used. Comparison of equations (14) and (17) shows two di€erences: (1) the signs of the third term on the right side are di€erent; and (2) the stress component in the second term on the right side of equation (17) is the longitudinal deviatoric stress on the basement (i.e. y = h) while the corresponding term in equation (14) is the depth-averaged longitudinal deviatoric stress. The sign di€erence is a trivial consequence of the

utop ˆ u0 ‡ u1 x ‡ u2 x2 vtop ˆ v0 ‡ v1 x

…18a†

where u0, u1, u2, v0 and v1 are constants. Assuming that displacements of the top and the base of the block are related, this boundary condition can represent indirectly basal underplating and basal erosion.

Table 2. Model parameters for rectangular and wedge-shaped blocks Rectangular block (L = 250 km, a = b = 38, r = 2800 Wedge block (L = 200 km, a = 38, b = 68, r = 2700 kg mÿ3) kg mÿ3) Parameter ÿ10

ÿ1

u0 (10 ms ) u1 (10ÿ15 sÿ1) u2 (10ÿ20 mÿ1 sÿ1) v0 (10ÿ10 m sÿ1) v1 (10ÿ15 sÿ1) t0 (MPa) t1 (MPa kmÿ1) T0 (K) T1 (K kmÿ1) A*0 (10ÿ23 Paÿ1 sÿ1) A*1 (10ÿ23 Paÿ1 sÿ1 kmÿ1) power-law rheology

A

B

31.71$ ÿ1 0 ÿ0.4 0 ÿ57.44 0 % % 5 0 %

31.71 ÿ1 0 ÿ0.4 0 ÿ57.44 0 750 2.5 % % anorthosite}

Case*

*Cases A, B, C, and D refer to Figs 4±7, respectively. A:Newtonian, constant Z(1022 Pa s), constant Çexx and constant v on surface, constant tb. B:Power-law, depth-dependent T, constant exx and constant v on surface, constant tb. C:Newtonian, constant Z(1022 Pa s), variable Çexx and variable v on surface, variable tb. D:Power-law, constant T, variable Çexx and variable v on surface, variable tb. $100 mm aÿ1131.71  10ÿ10 m sÿ1. %not applicable (constant Newtonian viscosity or computed power-law viscosity). }rheological parameters: A0=2.06  10ÿ23 Paÿ3.2 sÿ1, n = 3.2, Q = 238 kJ molÿ1.

C

D

0 ÿ1 ÿ2.25 ÿ1.5855 0.7927 ÿ80 0.4 % % 5 0 %

0 1 ÿ1 ÿ1.5855 0.7927 ÿ80 0.4 750 0 % % anorthosite}

1618

J. LIU and G. RANALLI

(b) The basal shear stress along the basement is a linear function of position tb ˆ t 0 ‡ t 1 x

…18b†

where t0 and t1 are constants. (c) The temperature within the block is either constant or varies linearly with depth T ˆ T0 ‡ T1 y

…18c†

where T0 and T1 are constants. (d) The rheology of the block can be nonlinear (stress exponent n>1) or linear (n = 1). In case of linear rheology, Newtonian viscosity is assumed to be of the form

1
Zˆ ÿ  2 A0 ‡ A1 y

…18d†

where and A*0 are A*1 constants. This allows the viscosity to be either constant or decreasing with depth. The ®rst example (Fig. 4) serves as a reference case for a rectangular block with constant Newtonian viscosity (Z = 1022 Pa s) and constant basal shear stress. The top velocity boundary condition is such that the u component varies linearly with x (resulting in a constant longitudinal strain rate of Çexx= ÿ 10ÿ15 sÿ1 on the upper surface of the block) and the v component is

Fig. 4. Stress and deformation within a rectangular block with Newtonian rheology (case A of Table 2). (a) Velocity ®eld; (b) maximum shear stress trajectories and type of potential faults; (c) tectonic (non-gravitational) normal and shear stresses along the rear of the block.

Stresses and velocities in orogenic wedges

1619

Fig. 5. Stress and deformation within a rectangular block with power-law rheology (case B of Table 2). (a), (b) & (c) as in Fig. 4.

constant. This velocity boundary condition is compatible with Nye's one-dimensional solution (Nye, 1957). Consequently, results for this case yield simply a twodimensional representation of Nye's one-dimensional solution in the case of Newtonian viscosity. Figure 4 shows (a) the velocity ®eld within the block; (b) the maximum shear stress trajectories, which coincide with potential faulting surfaces in a frictionless (plastic) material; (c) the depth distribution of the `non-gravitational' normal and shear stress components [i.e. the functions and F(y) and G(y)] along the rear vertical edge of the block (note that a tectonic compression is required to maintain the deformation regime). The second example (Fig. 5) considers a rectangular block of nonlinear rheology with the power-law creep parameters of anorthosite (cf. Ranalli, 1995) and the same velocity and basal shear stress boundary con-

ditions as in the previous case. The temperature increases from 750 K on the top to 850 K at the bottom of the block (the relatively large temperature at the top is necessary to yield realistic strain rates, given the rheology of the material). Even this low temperature gradient increases considerably the ductility of the material: predominant laminar ¯ow is found near the base of the block. F(y) decreases in magnitude non linearly from about 76 MPa at the top to about 0 MPa at the bottom of the block. The third example (Fig. 6) consists of a wedge with constant Newtonian viscosity (Z = 1022 Pa s) and a longitudinal component of the top velocity boundary condition resulting in a longitudinal strain rate Çexx varying from ÿ10ÿ15 sÿ1 at the rear to ÿ10ÿ14 sÿ1 at the front of the wedge. The velocity boundary conditions result in a dominant movement pattern of

1620

J. LIU and G. RANALLI

Fig. 6. Stress and deformation within a wedge-shaped block with Newtonian rheology (case C of Table 2). (a), (b) & (c) as in Fig. 4.

wedge material which is longitudinal at the front, and mainly vertical towards the rear (for instance, as resulting from increasing underplating). The maximum shear stress trajectories show thrust deformation throughout the wedge. F(y) increases in magnitude from 40 MPa at the top to 50 MPa at the bottom of the wedge; G(y) varies roughly in a parabolic way with depth. The fourth example (Fig. 7) considers a wedge with the power-law rheology of anorthosite and with a constant temperature of 750 K. The longitudinal strain rate Çexx on the top surface of the wedge varies linearly from ÿ10ÿ15 sÿ1 at the rear to ÿ3  10ÿ15 sÿ1 at the front. There is some underplating near the rear of the wedge as in the previous case (note the di€erent velocity scale in the ®gures). The maximum shear stress trajectories show coeval development of an extensional

regime (potential normal faulting) at the rear and a compressional regime (thrusting) at the front of the wedge. The non-gravitational normal stress at the rear of the wedge necessary to maintain this deformation regime is tensional: F(y) decreases nonlinearly from about 80 MPa on the top to about 10 MPa at the bottom of the wedge. The depth distribution of G(y) is similar to the previous case. DISCUSSION AND CONCLUSIONS We have presented an analytical linear longitudinal strain rate solution for blocks and wedges with nonlinear rheology (nr 1), which is a two-dimensional extension of the classic one-dimensional solution derived by Nye (1957) for ice sheets. This solution allows

Stresses and velocities in orogenic wedges

1621

Fig. 7. Stress and deformation within a wedge-shaped block with power-law rheology (case D of Table 2). (a), (b) & (c) As in Fig. 4.

the modelling of stresses and deformation in twodimensional orogenic wedges with deformation styles varying in the longitudinal direction. The basic kinematic assumption is that the longitudinal strain rate, at any depth, is a linear function of the distance from the rear of the block. Together with wedge segmentation, this assumption allows the simulation of varying deformation patterns, including coeval compression and extension, in di€erent parts of the wedge. The speci®cation of velocity boundary condition on the surface of the block signi®cantly simpli®es the implementation of the solution. A more complicated algorithm would be needed to specify the velocity boundary conditions on the basement rather than on the top surface. However, kinematic conditions on the top of an incompressible block may be taken as a ®rst approximation to conditions at the bot-

tom, for instance, with an upward component of velocity simulating underplating. In the solution, the nongravitational stresses F(y) and G(y) along the rear vertical boundary of the block are derived quantities, i.e. they are the tectonic stresses necessary to maintain the speci®ed linear longitudinal strain rate for a given rheology under speci®ed top (velocity) and bottom (shear stress) boundary conditions (this approach is the same as the one taken in the one-dimensional solution; Nye, 1957). Alternatively, they could be assigned a priori as part of the boundary conditions, but this is not possible in the framework of the present solution. It remains to be discussed, in any particular case, whether the required normal and shear tectonic stresses are realistic. In the four simple examples considered, their magnitudes are of the right order, being comprised between 280 MPa.

1622

J. LIU and G. RANALLI

Their distribution varies with depth, but is generally rather simple. If a priori knowledge of F(y) and G(y) is assumed, it is still possible, within the framework of the present solution, to modify iteratively the boundary conditions until the required tectonic stress distribution is satisfactorily approximated. The linear strain rate solution can also be applied to investigate the evolution of wedges with varying taper, by solving iteratively for di€erent surface slopes. In this context, it is also relevant to the problem of the exhumation of high pressure±low temperature metamorphic rocks near the rear of an unstable orogenic wedge (Rubie, 1984; Platt and Lister, 1985; Platt, 1993). Combined with a given erosion rate at the top surface of the wedge, the velocity expression of the solution [equation (13c)] provides an estimation of the exhumation rate and, therefore, of the time variation of the pressure±temperature conditions of rocks throughout the wedge. Assuming that the surface erosion rate is the same as the rate at which material is uplifted at the basement of the wedge, the exhumation rate at the rear of the wedge of the fourth example would be about 6 mm per y (Fig. 7). This is, of course, an upper boundary since it assumes that no additional topography is generated. Naturally, the rheology of orogenic wedges is more complicated than power-law throughout. As the temperature decreases approaching the surface, plastic± brittle deformation takes over. In principle, this change can be accommodated by the present solution, by letting the stress exponent increase (pure plasticity requires n 4 1), and ensuring that the deformation patterns in the power-law and plastic domains are compatible (Liu, 1996). Despite the limitations imposed by its analytical character, the solution presented in this paper has the advantage of generality, and provides an additional approach to the study of the tectonic evolution of orogenic wedges. AcknowledgementsÐThis paper represents part of the doctoral dissertation by one of the authors (J. L.). The work was supported by NSERC (Natural Sciences and Engineering Research Council of Canada) grants to G. R. and by a Carleton University Scholarship to J. L. We thank P.-Y. Robin and K. Rohr for constructive reviews.

REFERENCES Boyer, S. and Elliott, D. (1982) Thrust systems. American Association of Petroleum Geologists Bulletin 66, 1196±1230. Budd, W. F. (1970) The longitudinal stress and strain-rate gradients in ice masses. Journal of Glaciology 9, 19±47. Chapple, W. M. (1978) Mechanics of thin-skinned fold-and-thrust belts. Geological Society of America Bulletin 89, 1189±1198. Dahlen, F. A. (1984) Noncohesive critical Coulomb wedges: An exact solution. Journal of Geophysical Research 89, 10,125±10,133. Dahlen, F. A. (1990) Critical taper model of fold-and-thrust belts and accretionary wedges. Annual Review of Earth and Planetary Sciences 18, 55±99. Dahlen, F. A. and Barr, T. D. (1989) Brittle frictional mountain building 1, Deformation and mechanical energy budget. Journal of Geophysical Research 94, 3906±3922.

Dahlen, F. A. and Suppe, J. (1988) Mechanics, growth, and erosion of mountain belts. In Processes in Continental Lithospheric Deformation. eds S. P. Clark Jr, B. C. Burch®el and J. Suppe. pp. 161±208. Geological Society of America Special Paper, 218. Dahlen, F. A., Suppe, J. and Davis, D. M. (1984) Mechanics of fold-and-thrust belts and accretionary wedges: cohesive Coulomb theory. Journal of Geophysical Research 89, 10,087±10,101. Davis, D. M., Suppe, J. and Dahlen, F. A. (1983) Mechanics of fold-and-thrust belts and accretionary wedges. Journal of Geophysical Research 89, 1153±1172. Elliott, D. (1976) The motion of thrust sheets. Journal of Geophysical Research 81, 949±963. Liu, J. Y. (1996) Stress and Deformation in Orogenic Wedges with Power-Law Rheology. Ph.D. thesis. Carleton University, Ottawa. Liu, J. Y. and Ranalli, G. (1992) Stress in an overthrust sheet and propagation of thrusting: An Airy stress function solution. Tectonics 11, 549±559. Liu, J. Y. and Ranalli, G. (1994) A two-dimensional power-law rheology model for subduction-accretionary orogenic complexes (Abstract). EOS 75, 598. Nye, J. F. (1951) The ¯ow of glaciers and ice-sheets as a problem in plasticity. Proceedings Royal Society of London A±207, 554±572. Nye, J. F. (1957) The distribution of stress and velocity in glaciers and ice sheets. Proceedings Royal Society of London AA±239, 113± 133. Nye, J. F. (1969) The e€ect of longitudinal stress on the shear stress at the base of an ice sheet. Journal of Glaciology 8, 207±213. Pavlis, T. L. and Bruhn, R. L. (1983) Deep-seated ¯ow as a mechanism for the uplift of broad forearc ridges and its role in the exposure of high P/T metamorphic terranes. Tectonics 2, 473±497. Platt, J. P. (1986) Dynamics of orogenic wedges and the uplift of high-pressure metamorphic rocks. Geological Society of America Bulletin 97, 1037±1053. Platt, J. P. (1987) The uplift of high-pressure±low-temperature metamorphic rocks. Philosophical Transactions Royal Society of London, A±321, 87±103. Platt, J. P. (1993) Exhumation of high-pressure rocks. Terra Nova 5, 119±133. Platt, J. P. and Lister, G. S. (1985) Structural history of high-pressure metamorphic rocks in the Vanoise massif, French Alps, and their relation to Alpine tectonic events. Journal of Structural Geology 7, 19±36. Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P. (1992) Numerical recipes. 2nd edn. Cambridge University Press, Cambridge. Ranalli, G. (1995) Rheology of the Earth, 2nd edn. Chapman & Hall, London. Rubie, D. C. (1984) A thermal-tectonic model for high-pressure metamorphism and deformation in the Sesia zone, western Alps. Journal of Geology 92, 21±36. Shoemaker, E. M. and Morland, L. W. (1984) A glacier ¯ow model incorporating longitudinal deviatoric stress. Journal of Glaciology 30, 334±339. Spiegel, M. R. (1971) Schaum's Outline of Theory and Problems of Advanced Mathematics for Engineers and Scientists. McGraw-Hill, New York. Suppe, J. (1985) Principles of Structural Geology. Prentice-Hall, Englewood Cli€s. Wojtal, S. (1992a) One-dimensional models for plane and non-plane power-law ¯ow in shortening and elongating thrust zones. In Thrust tectonics, ed. K. R. McClay, pp. 41±52. Chapman & Hall, London. Wojtal, S. (1992b) Shortening and elongation of thrust zones within the Appalachian foreland fold-thrust belt. In Structural Geology of Fold and Thrust Belts, eds S. Mitra and G. W. Fisher, pp. 93±103. Johns Hopkins University Press, Baltimore. Yin, A. (1988) Geometry, kinematics, and a mechanical analysis of a strip of the Lewis allochthon from Peril Peak to Bison Mountain, Glacier National Park, Montana. Ph.D. thesis. University of Southern California, Los Angeles. Yin, A. (1989) Origin of regional, rooted low-angle normal faults: a mechanical model and its implications. Tectonics 8, 469±482. Yin, A. (1993) Mechanics of wedge-shaped fault blocks 1, an elastic solution for compressional wedges. Journal of Geophysical Research 98, 14,245±14,256.

Stresses and velocities in orogenic wedges Zhao, W. L., Davis, D. M., Dahlen, F. A. and Suppe, J. (1986) Origin of convex accretionary wedges: evidence from Barbados. Journal of Geophysical Research 91, 10,246±10,258.

APPENDIX In this Appendix, we derive the general solution given in the text [equation (13)]. Di€erentiating the ®rst of equation (11a) twice with respect to x and combining it with equation (6) leads to F…y† ÿ xG0 …y† @2 l @l ÿ G0 …y† ˆ 0 2 @x2 @x

…A1†

which is a second-order partial di€erential equation which can be transformed into a ®rst-order equation and then solved by standard procedures (cf. Spiegel, 1971). A general solution to equation (A1) is lˆ

k…y† ‡ f0 …y† F…y† ÿ xG0 …y†

2_

2_

8 gx y G…y† ˆ r > > < f …y†G0 …y† ˆ 2…a ‡ b y† 0 1 1 …c†  1 00 1 d2  > > f0 …y†F…y† ˆ 0 : k …y† ‡ 2 2 dy2

G0 …y† ˆ ÿ

…y G…y† ˆ ÿ

2_

0

 x d2   1 d2  f0 …y†F…y† ÿ f0 …y†G0 …y† 2 2 dy 2 dy2

 d      k…y†G0 …y† G0 …y†k…y† G0 …y†ÿrgx ‡ G…y†ÿ rgx y dy ‡2  2 F…y†ÿxG0 …y†    G0 …y†k…y† F0 …y†ÿxG00 …y† G…y†ÿ rgx y ÿ4 ˆ0  3 F…y†ÿxG0 …y†

…A3†

where the double prime denotes the second derivative with respect to y. Equation (A3) contains several terms that are functions of x of di€erent orders. In order that equation (A3) be satis®ed at all points, one of following three conditions must apply 8 < k…y† ˆ 0 …a† f0 …y†G0 …y† ˆ 2…a1 ‡ b1 y† …A4a† : f0 …y†F…y† ˆ ÿ2…a2 ‡ b2 y† 8 0 < G …y† ˆ 0  1 d2  …b† 1 00 f0 …y†F…y† ˆ 0 : 2 k …y† ‡ 2 dy2

…A4b†

a1 ‡ b1 y F…y† a2 ‡ b2 y

…A5†

Integrating the second of equations (A5) with respect to y gives

@ exy 1 00 @ exx @ eyy ˆ k …y† ‡ ÿ2 @y2 @x2 @x@y 2 ‡

…A4c†

where a1, a2, b1 and b2 are arbitrary constants. The second condition [equation (A4b)] represents a stress ®eld corresponding to ¯ow of a material compressed between two parallel plates (Nye, 1957), and is not applicable to an orogenic wedge with variable longitudinal deformation pattern. The third condition [equation (A4c)] is not applicable either, because it is a stress ®eld with zero shear stress everywhere in the wedge, which is not able to produce any shear deformation. Therefore, only the ®rst case [equation (A4a)] is relevant. Equation (A4a) gives ÿ
2 a2 ‡ b2 y l ˆ f0 …y† ˆ ÿ F…y†

…A2†

where k(y) and f0(y) are two integration functions. Combination of equations (A2) and (11a) with the equation of compatibility for two-dimensional plane strain (cf. e.g. Ranalli, 1995) gives

1623

a1 ‡ b1 t F…t†dt a2 ‡ b2 t

…A6†

where the traction-free boundary condition on the surface of the wedge [i.e. G(0) = 0] has been used to determine the integration constant. Equation (A6) implies that the normal and shear stress components of non-gravitational stress along the rear boundary of the wedge are not independent of each other. The velocity integration functions I(x) and H(y) in equation (12) can be determined by combining equations (A5), (A6) and (11a) with (2). After some algebra, we have b1 3 b2 2 x ‡ x ‡ b4 6 2 9 8 y  > …> = > a2 ‡ b2 t ; : F…z† I…x† ˆ

0

0

b a1 ÿ y3 ÿ y2 ÿ b3 y ÿ a3 6 2

…A7†

where a3, b3 and b4 are arbitrary constants. Thus, all the integration functions for stress, strain rate, and velocity have been expressed as functions of seven arbitrary constants and the function F(y), the horizontal nongravitational stress on the rear of the block.